% Make animation/movie of sampling the Z-transform of a single root around the 
% unit circle
% 2002-10-09 dpwe@ee.columbia.edu

% Do we want to write a movie? (otherwise, just display on screen)
makemovie = 0;

% Define the poles/zeros

% Single zero example
%zz = [0.8*exp(j*pi*0.3)].';
%pp = [].';
%moviename = 'onezero.mov';
%ymax = 2;

% Single pole example
%pp = [0.9*exp(j*pi*0.3)].';
%zz = [].';
%moviename = 'onepole.mov';
%ymax = 10;

% Two-pole, two-zero example
zz = [0.8*exp(j*pi*0.3) 0.8*exp(j*pi*-0.3)].';
pp = [0.9*exp(j*pi*0.3) 0.9*exp(j*pi*-0.3)].';
moviename = '2p2z.mov';
ymax = 2;

% zp2tf doesn't work for complex systems??
%[bb,aa] = zp2tf(zz,pp,1);
bb = poly(zz);
aa = poly(pp);

% Steps around the top half of the unit circle
ww = [0:200]/200*pi;

% Initialize the movie
if makemovie
    eval(['MakeQTMovie start ', moviename]);
end

% Lay out the display
subplot(121);
zplane(pp,zz);

% fixed axes for frq plot
subplot(222)
fax = [0 1 0 ymax];
axis(fax)
grid

% fixed axes for phs plot (pi-normalized)
subplot(224)
pax = [0 1 -1 1];
axis(pax);
grid

%HH = NaN*ww;
%PP = HH;
% Preset the plots to have the whole lines, 
% not just points so far 
GG = polyval(bb,exp(j*ww))./polyval(aa,exp(j*ww));
HH = abs(GG);
PP = angle(GG);

% plot the frames
for i = 1:length(ww);

  w = ww(i);
  z = exp(j*w);

  % Evaluate the z transform at this point
  Gz = polyval(bb,z)./polyval(aa,z);
  
  HH(i) = abs(Gz);
  PP(i) = angle(Gz);
  
  % Make the plots
  subplot(121)
  zplane([0],[]);
  hold on
  plot(real(z),imag(z),'sg');
  
  % Add omega parameters
  for www = -0.8:0.2:0.8
    ejw = exp(j*www*pi);
    ll = sprintf('%.1f\\pi',www);
    text(real(ejw), imag(ejw), ll);
  end
  
  % Plot and connect to all the poles and zeros
  for r = pp.'
	plot(real(r),imag(r), 'xr');
	plot([real(r), real(z)], [imag(r), imag(z)], '-r');
  end
  for r = zz.'
	plot(real(r),imag(r), 'ob');
	plot([real(r), real(z)], [imag(r), imag(z)], '-g');
  end
  hold off
  
  subplot(222)
  plot(ww/pi, HH, w/pi, HH(i), 'sg',[w/pi w/pi], [0 HH(i)], '-g');
  axis(fax)
  grid
  title('magnitude');

  subplot(224)
  plot(ww/pi, PP/pi, w/pi, PP(i)/pi, 'sg');
  axis(pax)
  grid
  title('phase');
  xlabel('\omega / \pi');
  set(gca, 'YTick', [-1 -0.5 0 0.5 1])
  set(gca, 'YTickLabel', [' -pi ';'-pi/2';'  0  ';' pi/2'; '  pi '])

  % Save this image
  if makemovie
    MakeQTMovie addfigure
  else
    drawnow
  end

%  if (w > 0.28*pi)
%    break
%  end

end

if makemovie
  MakeQTMovie finish
end