% M05-ztsurf.diary
% 3-D surface visualization of |G(z)|
% 2007-10-02

% Plot of 3-D |G(z)| surface
% Define a grid of regularly spaced complex values on the zplane
[x,y]=meshgrid(-1.8:0.05:1.8);
% The actual complex values
zzz = x + j*y;
% Define the z-transform by its poles and zeros...
% Two real zeros
zz = [0.8 -1.5];
% A pair of complex-conjugate poles
pp = [.9*exp(j*.6), .9*exp(j*-.6)];
% Convert these into numerator and denominator polynomial coefficients
[nn,dd]=zp2tf(zz',pp',1);
% Evaluate the numerator and denominator at every grid point
numv = polyval(nn,zzz);
denv = polyval(dd,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([-2 2 -2 2 -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(nn,ww)./polyval(dd,ww);
% Plot a red line just above the surface around the unit circle
plot3(real(ww), imag(ww), 1+20*log10(abs(gw)),'r')
% Add similar lines just above the surface for the Re and Im axes
uu = x(1,:);
gx = polyval(nn,uu)./polyval(dd,uu);
gy = polyval(nn,j*uu)./polyval(dd,j*uu);
plot3(uu, 0*uu, 1+20*log10(abs(gx)), 'b');
plot3(uu*0, uu, 1+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
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