## Copyright (C) 1999, 2001 Paul Kienzle ## Copyright (C) 2004 Stefan van der Walt ## ## This program is free software; you can redistribute it and/or modify it under ## the terms of the GNU General Public License as published by the Free Software ## Foundation; either version 3 of the License, or (at your option) any later ## version. ## ## This program is distributed in the hope that it will be useful, but WITHOUT ## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or ## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more ## details. ## ## You should have received a copy of the GNU General Public License along with ## this program; if not, see . ## -*- texinfo -*- ## @deftypefn {Function File} {} zplane (@var{z}, @var{p}) ## @deftypefnx {Function File} {} zplane (@var{b}, @var{a}) ## @deftypefnx {Function File} {} zplane (..., nods) ## ## Plot the poles and zeros. If the arguments are row vectors then they ## represent filter coefficients (numerator polynomial b and denominator ## polynomial a), but if they are column vectors or matrices then they ## represent poles and zeros. ## ## This is a horrid interface, but I didn't choose it; better would be ## to accept b,a or z,p,g like other functions. The saving grace is ## that poly(x) always returns a row vector and roots(x) always returns ## a column vector, so it is usually right. You must only be careful ## when you are creating filters by hand. ## ## Note that due to the nature of the roots() function, poles and zeros ## may be displayed as occurring around a circle rather than at a single ## point. ## ## The transfer function is ## ## @example ## @group ## B(z) b0 + b1 z^(-1) + b2 z^(-2) + ... + bM z^(-M) ## H(z) = ---- = -------------------------------------------- ## A(z) a0 + a1 z^(-1) + a2 z^(-2) + ... + aN z^(-N) ## ## b0 (z - z1) (z - z2) ... (z - zM) ## = -- z^(-M+N) ------------------------------ ## a0 (z - p1) (z - p2) ... (z - pN) ## @end group ## @end example ## ## The denominator a defaults to 1, and the poles p defaults to []. ## @end deftypefn ## 2/3/2018 modifications by Dr. Brian M. Reeder ## changes colour for each system plotted ## correctly marks multiple poles / zeros function zplane (z, p = []) if (nargin<1 || nargin>2) print_usage; endif if (columns (z)>1 || columns (p)>1) if (rows (z)<=1 || rows (p)<=1) if (isempty (z)) z = 1; endif # convert to poles / zeros if (isempty (p)) p = 1; endif M = length (z) - 1; N = length (p) - 1; z = [roots(z); zeros(N - M, 1)]; p = [roots(p); zeros(M - N, 1)]; endif endif a = axis; # plot range xmin = min ([a(1)/1.1; -1; real(z(:)); real(p(:))]); xmax = max ([a(2)/1.1; 1; real(z(:)); real(p(:))]); ymin = min ([a(3)/1.1; -1; imag(z(:)); imag(p(:))]); ymax = max ([a(4)/1.1; 1; imag(z(:)); imag(p(:))]); r = exp (2i*pi*[0:100]/100); # draw unit circle plot (real (r), imag (r), "k"); grid on; hold on; colororder = get (gca (), "colororder"); for c = 1:columns (z) # for every system ci = get (gca (), "colororderindex"); # get plot colour for zeros col = colororder (ci, :); plot_with_labels (c, z, col, "o"); set (gca (), "colororderindex", ci); # use same colour for poles plot_with_labels (c, p, col, "x"); endfor axis equal; axis (1.1*[xmin, xmax, ymin, ymax]); hold off endfunction function plot_with_labels (c, x, col, symbol) # plot poles or zeros if (! isempty (x)) x_u = unique (x(:,c)); plot (real (x_u), imag (x_u), symbol, "color", col); for i = 1:length (x_u) n = sum (x_u(i) == x(:,c)); if (n > 1) text (real (x_u(i)), imag (x_u(i)), [" " num2str(n)], "color", col); endif endfor endif endfunction %!demo %! ## construct target system: %! ## symmetric zero-pole pairs at r*exp(iw),r*exp(-iw) %! ## zero-pole singletons at s %! pw=[0.2, 0.4, 0.45, 0.95]; #pw = [0.4]; %! pr=[0.98, 0.98, 0.98, 0.96]; #pr = [0.85]; %! ps=[]; %! zw=[0.3]; # zw=[]; %! zr=[0.95]; # zr=[]; %! zs=[]; %! %! ## system function for target system %! p=[[pr, pr].*exp(1i*pi*[pw, -pw]), ps]'; %! z=[[zr, zr].*exp(1i*pi*[zw, -zw]), zs]'; %! M = length(z); N = length(p); %! sys_a = [ zeros(1, M-N), real(poly(p)) ]; %! sys_b = [ zeros(1, N-M), real(poly(z)) ]; %! disp("The first two graphs should be identical, with poles at (r,w)="); %! disp(sprintf(" (%.2f,%.2f)", [pr ; pw])); %! disp("and zeros at (r,w)="); %! disp(sprintf(" (%.2f,%.2f)", [zr ; zw])); %! disp("with reflection across the horizontal plane"); %! subplot(231); %! zplane(sys_b, sys_a); %! title("transfer function form"); %! subplot(232); %! zplane(z,p); %! title("pole-zero form"); %! subplot(233); %! zplane(z); %! title("empty p"); %! subplot(234); %! zplane(sys_b); %! title("empty a"); %! disp("The matrix plot has 2 sets of points, one inside the other"); %! subplot(235); %! zplane([z, 0.7*z], [p, 0.7*p]); %! title("matrix");