function [Hz,zgrid] = zplot(b,a,N,zmax)
%  zplot --> Plot of rational z-transforms.
%
%  <Synopsis>
%    [Hz,zgrid] = zplot(b,a)
%    [Hz,zgrid] = zplot(b,a,N)
%    [Hz,zgrid] = zplot(b,a,N,zmax)
%
%  <Description>
%    The function evaluates the z-transform for the rational function
%
%                                   -1                -m 
%               B(z)    b(1) + b(2)z  + .... + b(m+1)z
%        H(z) = ---- = ----------------------------------
%                                   -1                -n
%               A(z)    a(1) + a(2)z  + .... + a(n+1)z
%
%    given numerator and denominator coefficients in vectors b and a.
%    The argument zmax is the maximum magnitude of z for which the
%    z-transform H(z) is evaluated, and N is the number of grid points
%    equally spaced from z = 0 to zmax. If zmax and N isn't specified,
%    the defaults are N = 20 and zmax = 1.5.

%  <Revision>
%    Peter S.K. Hansen, IMM, Technical University of Denmark
%
%    Last revised: March 20, 2000
%-----------------------------------------------------------------------

% Check the required input arguments and set defaults.
if (nargin < 2)
  error('Not enough input arguments.')
end
if (nargin < 3)
  N = 20;
end
if (nargin < 4)
  zmax = 1.5;
end

% Make sure that parameter vectors a and b are row vectors.
a = a(:)'; b = b(:)';

% Grid of points in the complex z-plane.
zreal = ones(2*N+1,1)*(-zmax:(zmax/N):zmax);
zgrid = zreal + j*zreal';

% Evaluate the z-transform H(z)=B(z)/A(z) for all z in zgrid.
Hz = zgrid.^(length(a)-length(b)).*polyval(b,zgrid)./polyval(a,zgrid);

% Evaluate the Fourier-transform H(w)=B(w)/A(w).
[H,w] = freqz(b,a,2048,'whole');

% Suppress warnings for taking log to zero, i.e., in dB plots.
warning off;

% Graph of magnitude surface |H(z)|.
figure(1),mesh(zreal,zreal',20*log10(abs(Hz)));
view(45,45);
xlabel('Real Part, Re({\it z} )');
ylabel('Imaginary Part, Im({\it z} )');
zlabel('Magnitude, |{\it H(z)} |  [dB]');

% Add graph for |H(z)| evaluated on the unit-circle, i.e., |H(w)|.
hold on;
plot3(cos(w),sin(w),20*log10(abs(H)),'LineWidth',1);
hold off;

% Graph of magnitude response |H(w)|.
figure(2),subplot(2,1,1),plot(w/2/pi,20*log10(abs(H)));
axis([0 1 -inf inf]),grid;
xlabel('Normalized Frequency, {\it f}  [cycles/sample]');
ylabel('Magnitude, |{\it H(\omega)} |  [dB]');

% Zero-pole plot of H(z).
subplot('Position',[0.35 0.11 0.30 0.34]),zplane(b,a);
title('Zero-Pole plot')

warning on;

%-----------------------------------------------------------------------
% End of function zplot
%-----------------------------------------------------------------------