% 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