octave_packages/signal-1.1.3/bilinear.m

   1 ## Copyright (C) 1999 Paul Kienzle <pkienzle@users.sf.net>
   2 ##
   3 ## This program is free software; you can redistribute it and/or modify it under
   4 ## the terms of the GNU General Public License as published by the Free Software
   5 ## Foundation; either version 3 of the License, or (at your option) any later
   6 ## version.
   7 ##
   8 ## This program is distributed in the hope that it will be useful, but WITHOUT
   9 ## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  10 ## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
  11 ## details.
  12 ##
  13 ## You should have received a copy of the GNU General Public License along with
  14 ## this program; if not, see <http://www.gnu.org/licenses/>.
  15
  16 ## usage: [Zz, Zp, Zg] = bilinear(Sz, Sp, Sg, T)
  17 ##        [Zb, Za] = bilinear(Sb, Sa, T)
  18 ##
  19 ## Transform a s-plane filter specification into a z-plane
  20 ## specification. Filters can be specified in either zero-pole-gain or
  21 ## transfer function form. The input form does not have to match the
  22 ## output form. 1/T is the sampling frequency represented in the z plane.
  23 ##
  24 ## Note: this differs from the bilinear function in the signal processing
  25 ## toolbox, which uses 1/T rather than T.
  26 ##
  27 ## Theory: Given a piecewise flat filter design, you can transform it
  28 ## from the s-plane to the z-plane while maintaining the band edges by
  29 ## means of the bilinear transform.  This maps the left hand side of the
  30 ## s-plane into the interior of the unit circle.  The mapping is highly
  31 ## non-linear, so you must design your filter with band edges in the
  32 ## s-plane positioned at 2/T tan(w*T/2) so that they will be positioned
  33 ## at w after the bilinear transform is complete.
  34 ##
  35 ## The following table summarizes the transformation:
  36 ##
  37 ## +---------------+-----------------------+----------------------+
  38 ## | Transform     | Zero at x             | Pole at x            |
  39 ## |    H(S)       |   H(S) = S-x          |    H(S)=1/(S-x)      |
  40 ## +---------------+-----------------------+----------------------+
  41 ## |       2 z-1   | zero: (2+xT)/(2-xT)   | zero: -1             |
  42 ## |  S -> - ---   | pole: -1              | pole: (2+xT)/(2-xT)  |
  43 ## |       T z+1   | gain: (2-xT)/T        | gain: (2-xT)/T       |
  44 ## +---------------+-----------------------+----------------------+
  45 ##
  46 ## With tedious algebra, you can derive the above formulae yourself by
  47 ## substituting the transform for S into H(S)=S-x for a zero at x or
  48 ## H(S)=1/(S-x) for a pole at x, and converting the result into the
  49 ## form:
  50 ##
  51 ##    H(Z)=g prod(Z-Xi)/prod(Z-Xj)
  52 ##
  53 ## Please note that a pole and a zero at the same place exactly cancel.
  54 ## This is significant since the bilinear transform creates numerous
  55 ## extra poles and zeros, most of which cancel. Those which do not
  56 ## cancel have a "fill-in" effect, extending the shorter of the sets to
  57 ## have the same number of as the longer of the sets of poles and zeros
  58 ## (or at least split the difference in the case of the band pass
  59 ## filter). There may be other opportunistic cancellations but I will
  60 ## not check for them.
  61 ##
  62 ## Also note that any pole on the unit circle or beyond will result in
  63 ## an unstable filter.  Because of cancellation, this will only happen
  64 ## if the number of poles is smaller than the number of zeros.  The
  65 ## analytic design methods all yield more poles than zeros, so this will
  66 ## not be a problem.
  67 ##
  68 ## References:
  69 ##
  70 ## Proakis & Manolakis (1992). Digital Signal Processing. New York:
  71 ## Macmillan Publishing Company.
  72
  73 function [Zz, Zp, Zg] = bilinear(Sz, Sp, Sg, T)
  74
  75   if nargin==3
  76     T = Sg;
  77     [Sz, Sp, Sg] = tf2zp(Sz, Sp);
  78   elseif nargin!=4
  79     print_usage;
  80   end
  81
  82   p = length(Sp);
  83   z = length(Sz);
  84   if z > p || p==0
  85     error("bilinear: must have at least as many poles as zeros in s-plane");
  86   end
  87
  88 ## ----------------  -------------------------  ------------------------
  89 ## Bilinear          zero: (2+xT)/(2-xT)        pole: (2+xT)/(2-xT)
  90 ##      2 z-1        pole: -1                   zero: -1
  91 ## S -> - ---        gain: (2-xT)/T             gain: (2-xT)/T
  92 ##      T z+1
  93 ## ----------------  -------------------------  ------------------------
  94   Zg = real(Sg * prod((2-Sz*T)/T) / prod((2-Sp*T)/T));
  95   Zp = (2+Sp*T)./(2-Sp*T);
  96   if isempty(Sz)
  97     Zz = -ones(size(Zp));
  98   else
  99     Zz = [(2+Sz*T)./(2-Sz*T)];
 100     Zz = postpad(Zz, p, -1);
 101   end
 102
 103   if nargout==2, [Zz, Zp] = zp2tf(Zz, Zp, Zg); endif
 104 endfunction