function [hwr,wwr,f] = dosurf(coefs);
%
% DOSURF  3-D view of complex FIR polynomial magnitude
%    dosurf(B) creates a figure showing the magnitude of the 
%    z-transform expression defined by the polynomial coefficients B.
%    B is assumed to be ascending in powers of z^-1
%
%    [hw,ww]=dosurf(B) creates the figure and returns the frequency
%    response magnitude vector hw at the frequencies ww.
%
%    [hw,ww,f]=dosurf(B) also returns the sequence of movie frames, f, that
%    can be used to see a helical animation of the 3-D plot using movie(f).
%
% RCM 5 March 2004
% (let me know if you find any bugs or errors!  rob.maher@montana.edu)
%

%set up z-plane grid
xx=-4:0.08:4;
xx2=0:0.08:4;
zz=ones(1,length(xx))'*(xx) +j*(xx)'*ones(1,length(xx));

%set up unit circle
ww= -pi:pi/400:pi;
uc=exp(j*ww);

%define z-transform (complex result)
h= zfunc(coefs,zz);

%define freq resp
hw= zfunc(coefs,uc);
topz=max(abs(hw));

%set up the figure:  response mag vs. normalized freq. in upper left corner
subplot('position',[0.1 0.7 0.2 0.25])
plot(ww,abs(hw));
xlabel('Normalized Radian Frequency');
ylabel('Magnitude');

%show pole/zero diagram in lower left corner
subplot('position',[0.1 0.2 0.2 0.4])
zplane(coefs);

%show 3-D plot of z-trans magnitude on right side
subplot('position',[0.3 0.05 0.7 0.9])
%show magnitude surface
surf(xx,xx,abs(h),'EdgeColor',[.8 .8 .82]);
hold on;

%indicate frequency response contour (z-trans on unit circle)
plot3(real(uc),imag(uc),abs(hw)+0.05,'r','LineWidth',2)

%show unit circle on imag vs. real face
plot3(real(uc),imag(uc),zeros(1,length(uc)),'b','LineWidth',1)

%define positive real axis contour
ax= zfunc(coefs,xx2);

plot3([0 max(xx2)],[0 0],[0 0],'g');
plot3(xx2,zeros(1,length(xx2)),abs(ax)+0.05,'g','LineWidth',2)
hold off;
colormap(gray);

%
%now do some scaling tweaks:
%
axis([-4 4 -4 4 0 1.25*topz]);
axis vis3d;
caxis([0 1.3*topz]);
xlabel('Real');
ylabel('Imag');

%
%finally, print the z-transform expression above the graph
%
str=sprintf('H(z)= %0.3g',coefs(1));

%only print the non-zero coefficients, and suppress printing unity
for k=1:(length(coefs)-1)
    if(coefs(k+1)<0)
        str=[str sprintf(' - ')];
    else if( coefs(k+1)~=0 )
            str=[str sprintf(' + ')];
        end
    end
    if(abs(coefs(k+1))==1)
        str=[str sprintf('z^{-%d}',k)];
    else if(abs(coefs(k+1))~=0)
            str=[str sprintf('%0.3gz^{-%d}',abs(coefs(k+1)),k)];
        end
    end

end
title(str);

view(20,12);

%
%if the user asked for output, make the assignments
if nargout > 0
    hwr=hw;
end
if nargout > 1
    wwr=ww;
end

%use getframe to make a movie sequence
if nargout > 2
    for k=1:25
        f(k)=getframe;
        view(10+(k*360/22),70-(k*60/22));
    end
end  % END of function dosurf


function h = zfunc(coefs,dat)
% ZFUNC
% h=zfunc(coefs,dat)
% function to compute h=coefs(1) + coefs(2)*dat^-1 + coefs(3)*dat^-2 + ...
%
% RCM 5 March 2004

h=coefs(1);
if( length(coefs)>1)
    for i=2:length(coefs)
        h = h + coefs(i)*dat.^-(i-1);
    end
end