function [hwr,wwr,f] = zview(b,a);
%
% ZVIEW  3-D view of complex FIR polynomial magnitude
%    zview(B,A) creates a figure showing the magnitude of the 
%    z-transform expression defined by the numerator polynomial coefficients B and
%    the denominator polynomial coefficients A.
%    A and B are assumed to be ascending in powers of z^-1
%
%    [hw,ww]=zview(B,A) creates the figure and returns the frequency
%    response magnitude vector hw at the frequencies ww.
%
%    [hw,ww,f]=zview(B,A) 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 6 March 2006  (original version started 11 March 2004)
% (let me know if you find any bugs or errors!  rob.maher@montana.edu)
%
if nargin < 1
    error('usage: zview(b,a)');
end
if nargin < 2
    a=1;
end

%clear existing figure or open new figure window
clf

%find pole or zero that is farthest from origin (for scaling)
maxz= max( [2; max( [abs(roots(b)); abs(roots(a))]) ]);

%set up z-plane grid
xx=-maxz:maxz/60:maxz;
%xx2=0:maxz/60:maxz;
xx2=maxz/60:maxz/60:maxz;
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(b,a,zz);

%define frequency response (evaluate z-transform along unit circle contour)
hw= zfunc(b,a,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(b,a);
hold on;
%show positive real axis in green
plot([0 1],[0 0],'g');
hold off;

%show 3-D plot of z-trans magnitude on right side
h3d=subplot('position',[0.3 0.05 0.7 0.8]);
%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.02*topz,'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, and show on surface in green
ax= zfunc(b,a,xx2);

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

%
%now do some scaling tweaks:
%
axis([-1.25*maxz 1.25*maxz -1.25*maxz 1.25*maxz 0 1.5*topz]);
axis vis3d;
caxis([0 1.3*topz]);
xlabel('Real');
ylabel('Imag');
view(20,12);

%
%finally, print the z-transform expression above the 3-D graph
%
subplot('position',[0.3 0.85 0.7 0.15])

str=sprintf('%0.3g',b(1));

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

end
text(.5,.5,str,'horizontalalignment','center','verticalalignment','bottom');  %,'backgroundcolor','w');
str=sprintf('%0.3g',a(1));

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

end
text(.5,.5,str,'horizontalalignment','center','verticalalignment','top'); %,'backgroundcolor','w');
line([.3 .7],[.5 .5],'color','k');
axis([0 1 0 1]);
axis off;

%
%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
    subplot(h3d);
    for k=1:25
        f(k)=getframe;
        view(10+(k*360/22),70-(k*60/22));
    end
end  % END of function zview


function h = zfunc(num,den,dat)
% ZFUNC
% h=zfunc(num,den,dat)
% function to compute h=coefs(1) + coefs(2)*dat^-1 + coefs(3)*dat^-2 + ...
%
% RCM 5 March 2004
hn=num(1);
if( length(num)>1)
    for i=2:length(num)
        hn = hn + num(i)*dat.^-(i-1);
    end
end
hd=den(1);
if( length(den)>1)
    for i=2:length(den)
        hd = hd + den(i)*dat.^-(i-1);
    end
end
% hd = not(isnan(hd) + isinf(hd)).*hd;
hd =  1e-5*(hd==0) + hd;
% hn = not(isnan(hn) + isinf(hn)).*hn;
h=hn./hd;
% end of function zfunc