octave_packages/signal-1.1.3/zplane.m

   1 ## Copyright (C) 1999, 2001 Paul Kienzle <pkienzle@users.sf.net>
   2 ## Copyright (C) 2004 Stefan van der Walt <stefan@sun.ac.za>
   3 ##
   4 ## This program is free software; you can redistribute it and/or modify it under
   5 ## the terms of the GNU General Public License as published by the Free Software
   6 ## Foundation; either version 3 of the License, or (at your option) any later
   7 ## version.
   8 ##
   9 ## This program is distributed in the hope that it will be useful, but WITHOUT
  10 ## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  11 ## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
  12 ## details.
  13 ##
  14 ## You should have received a copy of the GNU General Public License along with
  15 ## this program; if not, see <http://www.gnu.org/licenses/>.
  16
  17 ## usage: zplane(b [, a]) or zplane(z [, p])
  18 ##
  19 ## Plot the poles and zeros.  If the arguments are row vectors then they
  20 ## represent filter coefficients (numerator polynomial b and denominator
  21 ## polynomial a), but if they are column vectors or matrices then they
  22 ## represent poles and zeros.
  23 ##
  24 ## This is a horrid interface, but I didn't choose it; better would be
  25 ## to accept b,a or z,p,g like other functions.  The saving grace is
  26 ## that poly(x) always returns a row vector and roots(x) always returns
  27 ## a column vector, so it is usually right.  You must only be careful
  28 ## when you are creating filters by hand.
  29 ##
  30 ## Note that due to the nature of the roots() function, poles and zeros
  31 ## may be displayed as occurring around a circle rather than at a single
  32 ## point.
  33 ##
  34 ## The transfer function is
  35 ##
  36 ##        B(z)   b0 + b1 z^(-1) + b2 z^(-2) + ... + bM z^(-M)
  37 ## H(z) = ---- = --------------------------------------------
  38 ##        A(z)   a0 + a1 z^(-1) + a2 z^(-2) + ... + aN z^(-N)
  39 ##
  40 ##               b0          (z - z1) (z - z2) ... (z - zM)
  41 ##             = -- z^(-M+N) ------------------------------
  42 ##               a0          (z - p1) (z - p2) ... (z - pN)
  43 ##
  44 ## The denominator a defaults to 1, and the poles p defaults to [].
  45
  46 ## TODO: Consider a plot-like interface:
  47 ## TODO:       zplane(x1,y1,fmt1,x2,y2,fmt2,...)
  48 ## TODO:    with y_i or fmt_i optional as usual.  This would allow
  49 ## TODO:    legends and control over point colour and filters of
  50 ## TODO:    different orders.
  51 function zplane(z, p = [])
  52
  53   if (nargin < 1 || nargin > 2)
  54     print_usage;
  55   end
  56   if columns(z)>1 || columns(p)>1
  57     if rows(z)>1 || rows(p)>1
  58       ## matrix form: columns are already zeros/poles
  59     else
  60       ## z -> b
  61       ## p -> a
  62       if isempty(z), z=1; endif
  63       if isempty(p), p=1; endif
  64
  65       M = length(z) - 1;
  66       N = length(p) - 1;
  67       z = [ roots(z); zeros(N - M, 1) ];
  68       p = [ roots(p); zeros(M - N, 1) ];
  69     endif
  70   endif
  71
  72
  73   xmin = min([-1; real(z(:)); real(p(:))]);
  74   xmax = max([ 1; real(z(:)); real(p(:))]);
  75   ymin = min([-1; imag(z(:)); imag(p(:))]);
  76   ymax = max([ 1; imag(z(:)); imag(p(:))]);
  77   xfluff = max([0.05*(xmax-xmin), (1.05*(ymax-ymin)-(xmax-xmin))/10]);
  78   yfluff = max([0.05*(ymax-ymin), (1.05*(xmax-xmin)-(ymax-ymin))/10]);
  79   xmin = xmin - xfluff;
  80   xmax = xmax + xfluff;
  81   ymin = ymin - yfluff;
  82   ymax = ymax + yfluff;
  83
  84   text();
  85   plot_with_labels(z, "o");
  86   plot_with_labels(p, "x");
  87   refresh;
  88
  89   r = exp(2i*pi*[0:100]/100);
  90   plot(real(r), imag(r),'k'); hold on;
  91   axis equal;
  92   grid on;
  93   axis(1.05*[xmin, xmax, ymin, ymax]);
  94   if (!isempty(p))
  95     h = plot(real(p), imag(p), "bx");
  96     set (h, 'MarkerSize', 7);
  97   endif
  98   if (!isempty(z))
  99     h = plot(real(z), imag(z), "bo");
 100     set (h, 'MarkerSize', 7);
 101   endif
 102   hold off;
 103 endfunction
 104
 105 function plot_with_labels(x, symbol)
 106   if ( !isempty(x) )
 107
 108     x_u = unique(x(:));
 109
 110     for i = 1:length(x_u)
 111       n = sum(x_u(i) == x(:));
 112       if (n > 1)
 113         text(real(x_u(i)), imag(x_u(i)), [" " num2str(n)]);
 114        endif
 115     endfor
 116
 117     col = "rgbcmy";
 118     for c = 1:columns(x)
 119       plot(real( x(:,c) ), imag( x(:,c) ), [col(mod(c,6)),symbol ";;"]);
 120     endfor
 121
 122   endif
 123 endfunction
 124
 125 %!demo
 126 %! ## construct target system:
 127 %! ##   symmetric zero-pole pairs at r*exp(iw),r*exp(-iw)
 128 %! ##   zero-pole singletons at s
 129 %! pw=[0.2, 0.4, 0.45, 0.95];   #pw = [0.4];
 130 %! pr=[0.98, 0.98, 0.98, 0.96]; #pr = [0.85];
 131 %! ps=[];
 132 %! zw=[0.3];  # zw=[];
 133 %! zr=[0.95]; # zr=[];
 134 %! zs=[];
 135 %!
 136 %! ## system function for target system
 137 %! p=[[pr, pr].*exp(1i*pi*[pw, -pw]), ps]';
 138 %! z=[[zr, zr].*exp(1i*pi*[zw, -zw]), zs]';
 139 %! M = length(z); N = length(p);
 140 %! sys_a = [ zeros(1, M-N), real(poly(p)) ];
 141 %! sys_b = [ zeros(1, N-M), real(poly(z)) ];
 142 %! disp("The first two graphs should be identical, with poles at (r,w)=");
 143 %! disp(sprintf(" (%.2f,%.2f)", [pr ; pw]));
 144 %! disp("and zeros at (r,w)=");
 145 %! disp(sprintf(" (%.2f,%.2f)", [zr ; zw]));
 146 %! disp("with reflection across the horizontal plane");
 147 %! subplot(231); title("transfer function form"); zplane(sys_b, sys_a);
 148 %! subplot(232); title("pole-zero form"); zplane(z,p);
 149 %! subplot(233); title("empty p"); zplane(z);
 150 %! subplot(234); title("empty a"); zplane(sys_b);
 151 %! disp("The matrix plot has 2 sets of points, one inside the other");
 152 %! subplot(235); title("matrix"); zplane([z, 0.7*z], [p, 0.7*p]);