Digital Signal Processing-Second Edition-Matlab program

function [h1, h2, h3]=pzplot(z,p,ax); %PZPLOT Z-plane zero-pole plot. % PZPLOT(Z,P) plots the zeros Z and poles P (in column vectors) with the % unit circle for reference. Each zero is represented with a 'o' and % each pole with a 'x' on the plot. Multiple zeros and poles are % indicated by the multiplicity number shown to the upper right of the % zero or pole. PZPLOT(Z,P) where Z and/or P is a matrix plots the zeros % or poles in different columns with different colors. % % PZPLOT(B,A) where B and A are row vectors containing transfer function % polynomial coefficients plots the poles and zeros of B(z)/A(z). Note % that if B and A are both scalars they will be interpreted as Z and P. % % [H1,H2,H3]=PZPLOT(Z,P) returns vectors of handles to the lines and % text objects generated. H1 is a vector of handles to the zeros lines, % H2 is a vector of handles to the poles lines, and H3 is a vector of % handles to the axes / unit circle line and to text objects which are % present when there are multiple zeros or poles. In case there are no % zeros or no poles, H1 or H2 is set to the empty matrix []. % % PZPLOT(Z,P,AX) puts the plot into axes AX. % % See also FREQZ. % Author(s): T. Krauss, 3-19-93 % Copyright (c) 1988-98 by The MathWorks, Inc. % $Revision: 1.5 $ $Date: 1998/12/23 22:37:03 $ error(nargchk(1,2,nargin)) tf = 0; % By default assume inputs are poles and zeros prodsz = prod(size(z)); colsz = size(z,2); switch nargin, case 1, % Only Z was specified, check if it is indeed Z or if it is a numerator if (prodsz==colsz)&(~(prodsz == 1)), % Row vector, num specified p = 1; tf = 1; % Set TF flag else % Scalar or column vector, interpret as Z p = []; end case 2, prodsp = prod(size(p)); colsp = size(p,2); if xor(prodsz == 1,prodsp == 1), % One of the inputs is a scalar but not both, if (prodsz==colsz)&(prodsp==colsp), % The other input is a row, TF specified tf = 1; % Set TF flag end elseif ~((prodsz == 1)&(prodsp == 1)), % No scalars specified if (prodsz==colsz)&(prodsp==colsp), % num and den specified tf = 1; % Set TF flag elseif (prodsz==colsz)|(prodsp==colsp), % A row vector was specified with a column vector, this is not allowed error('When specifying polynomials, both vectors must be rows.'); end end end if tf == 1, % TF was specified, convert to z,p,k % Catch cases when the num or den are zero. if max(abs(p)) == 0, % Divide by zero not allowed error('Denominator cannot be zero.'); elseif max(abs(z)) == 0, % num is zero, there are no poles p = 1; end % Pad A or B with trailing zeros if B and A are of different length if length(z) < length(p) z = [z zeros(1,length(p)-length(z))]; elseif length(p) < length(z) p = [p zeros(1,length(z)-length(p))]; end % Remove trailing zeros if both num and den have them while z(end) == 0 & p(end) == 0, z(end) = []; p(end) = []; end % Find Poles and Zeros z = roots(z); p = roots(p); end if ~any(imag(z)), z = z + j*1e-50; end; if ~any(imag(p)), p = p + j*1e-50; end; % equivalent of 'hold on': ax = gca; set(ax,'nextplot','add') set(get(ax,'parent'),'nextplot','add') if ~isempty(z), zh = plot(z,'*','markersize',7); else zh = []; end if ~isempty(p), ph = plot(p,'+','markersize',6); else ph = []; end theta = linspace(0,2*pi,70); oh = plot(cos(theta),sin(theta),':'); %set(oh,'xdat',[get(oh,'xdat') NaN ... %xl(1)-diff(xl)*100 xl(2)+diff(xl)*100 NaN 0 0]); %set(oh,'ydat',[get(oh,'ydat') NaN 0 0 NaN ... %yl(1)-diff(yl)*100 yl(2)+diff(yl)*100]); handle_counter = 2; fuzz = diff(xl)/80; % horiz spacing between '*' or '+' and number fuzz=0; [r,c]=size(z); if (r>1)&(c>1), % multiple columns in z ZEE=z; else ZEE=z(:); c = min(r,c); end; for which_col = 1:c, % for each column of ZEE ... z = ZEE(:,which_col); [mz,z_ind]=mpoles(z); for i=2:max(mz), j=find(mz==i); for k=1:length(j), x = real(z(z_ind(j(k)))) + fuzz; y = imag(z(z_ind(j(k)))); if (j(k)~=length(z)), if (mz(j(k)+1)<mz(j(k))), oh(handle_counter) = text(x,y,num2str(i)); handle_counter = handle_counter + 1; end else oh(handle_counter) = text(x,y,num2str(i)); handle_counter = handle_counter + 1; end end end end [r,c]=size(p); if (r>1)&(c>1), % multiple columns in z PEE=p; else PEE=p(:); c = min(r,c); end; for which_col = 1:c, % for each column of PEE ... p = PEE(:,which_col); [mp,p_ind]=mpoles(p); for i=2:max(mp), j=find(mp==i); for k=1:length(j), x = real(p(p_ind(j(k)))) + fuzz; y = imag(p(p_ind(j(k)))); if (j(k)~=length(p)), if (mp(j(k)+1)<mp(j(k))), oh(handle_counter) = text(x,y,num2str(i)); handle_counter = handle_counter + 1; end else oh(handle_counter) = text(x,y,num2str(i)); handle_counter = handle_counter + 1; end end end end set(oh(2:length(oh)),'vertical','bottom'); if (nargout==1), h1 = zh; elseif (nargout==2), h1 = zh; h2 = ph; elseif (nargout==3), h1 = zh; h2 = ph; h3 = oh; end set(ax,'nextplot','replace') set(get(ax,'parent'),'nextplot','replace')