function [H,x,y] = vizZTsurf(b, a, maxx, nsteps)
% [H,x,y] = vizZTsurf(b, a, maxx, nsteps)
%      Visualize the ZT surface defined by b (zeros) and a (poles).
%      H returns a 2D complex array of the surface samples.
%      x and y are the corresponding sample points, which range
%      from -maxx to +maxx (default 1.5) in nsteps steps (default
%      300).
% 2013-10-30 Dan Ellis dpwe@ee.columbia.edu

if nargin < 3; maxx = 1.5; end
if nargin < 4; nsteps = 300; end

% Plot of 3-D |G(z)| surface
% Define a grid of regularly spaced complex values on the zplane
[x,y]=meshgrid(-maxx:(2*maxx/nsteps):maxx);
% The actual complex values
zzz = x + j*y;
% Define the z-transform by its poles and zeros...

% Evaluate the numerator and denominator at every grid point
numv = polyval(b,zzz);
denv = polyval(a,zzz);
% .. so the Z-transform is the ratio of the two
gg = numv ./ denv;
% Plot the surface
subplot(111)
% Note: vertical axis in dB (log scale) so zeros go to -Inf
surfl(x, y, 20*log10(abs(gg)))
shading interp
colormap(gray)
axis([-maxx maxx -maxx maxx -30 30])
hold on
% Add the unit circle
% Define the values of z around the unit circle
ww = exp(j*[0:100]*2*pi/100);
% Evaluate the ZT function on the unit circle
gw = polyval(b,ww)./polyval(a,ww);
% Plot a red line just above the surface around the unit circle
plot3(real(ww), imag(ww), 0.2+20*log10(abs(gw)),'r')
% Add similar lines just above the surface for the Re and Im axes
uu = x(1,:);
gx = polyval(b,uu)./polyval(a,uu);
gy = polyval(b,j*uu)./polyval(a,j*uu);
plot3(uu, 0*uu, 0.2+20*log10(abs(gx)), 'b');
plot3(uu*0, uu, 0.2+20*log10(abs(gy)), 'b');
% Plot a 'flat' z-plane at height -30: axes and unit circle
plot3(uu*0, uu, -30*ones(1,length(uu)), 'k');
plot3(uu, 0*uu, -30*ones(1,length(uu)), 'k');
plot3(real(ww), imag(ww), -30*ones(1,length(ww)),'k')
% Add the actual poles and zeros on the 2D zplane
zz = roots(b);
pp = roots(a);
plot3(real(zz), imag(zz), [-30 -30], 'ob')
plot3(real(pp), imag(pp), [-30 -30], 'xr')
% Now try spinning it in 3D to get a sense of the surface...
axis vis3d
hold off