octave_packages/signal-1.1.3/sftrans.m

   1 ## Copyright (C) 1999-2001 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: [Sz, Sp, Sg] = sftrans(Sz, Sp, Sg, W, stop)
  17 ##
  18 ## Transform band edges of a generic lowpass filter (cutoff at W=1)
  19 ## represented in splane zero-pole-gain form.  W is the edge of the
  20 ## target filter (or edges if band pass or band stop). Stop is true for
  21 ## high pass and band stop filters or false for low pass and band pass
  22 ## filters. Filter edges are specified in radians, from 0 to pi (the
  23 ## nyquist frequency).
  24 ##
  25 ## Theory: Given a low pass filter represented by poles and zeros in the
  26 ## splane, you can convert it to a low pass, high pass, band pass or
  27 ## band stop by transforming each of the poles and zeros individually.
  28 ## The following table summarizes the transformation:
  29 ##
  30 ## Transform         Zero at x                  Pole at x
  31 ## ----------------  -------------------------  ------------------------
  32 ## Low Pass          zero: Fc x/C               pole: Fc x/C
  33 ## S -> C S/Fc       gain: C/Fc                 gain: Fc/C
  34 ## ----------------  -------------------------  ------------------------
  35 ## High Pass         zero: Fc C/x               pole: Fc C/x
  36 ## S -> C Fc/S       pole: 0                    zero: 0
  37 ##                   gain: -x                   gain: -1/x
  38 ## ----------------  -------------------------  ------------------------
  39 ## Band Pass         zero: b ± sqrt(b^2-FhFl)   pole: b ± sqrt(b^2-FhFl)
  40 ##        S^2+FhFl   pole: 0                    zero: 0
  41 ## S -> C --------   gain: C/(Fh-Fl)            gain: (Fh-Fl)/C
  42 ##        S(Fh-Fl)   b=x/C (Fh-Fl)/2            b=x/C (Fh-Fl)/2
  43 ## ----------------  -------------------------  ------------------------
  44 ## Band Stop         zero: b ± sqrt(b^2-FhFl)   pole: b ± sqrt(b^2-FhFl)
  45 ##        S(Fh-Fl)   pole: ±sqrt(-FhFl)         zero: ±sqrt(-FhFl)
  46 ## S -> C --------   gain: -x                   gain: -1/x
  47 ##        S^2+FhFl   b=C/x (Fh-Fl)/2            b=C/x (Fh-Fl)/2
  48 ## ----------------  -------------------------  ------------------------
  49 ## Bilinear          zero: (2+xT)/(2-xT)        pole: (2+xT)/(2-xT)
  50 ##      2 z-1        pole: -1                   zero: -1
  51 ## S -> - ---        gain: (2-xT)/T             gain: (2-xT)/T
  52 ##      T z+1
  53 ## ----------------  -------------------------  ------------------------
  54 ##
  55 ## where C is the cutoff frequency of the initial lowpass filter, Fc is
  56 ## the edge of the target low/high pass filter and [Fl,Fh] are the edges
  57 ## of the target band pass/stop filter.  With abundant tedious algebra,
  58 ## you can derive the above formulae yourself by substituting the
  59 ## transform for S into H(S)=S-x for a zero at x or H(S)=1/(S-x) for a
  60 ## pole at x, and converting the result into the form:
  61 ##
  62 ##    H(S)=g prod(S-Xi)/prod(S-Xj)
  63 ##
  64 ## The transforms are from the references.  The actual pole-zero-gain
  65 ## changes I derived myself.
  66 ##
  67 ## Please note that a pole and a zero at the same place exactly cancel.
  68 ## This is significant for High Pass, Band Pass and Band Stop filters
  69 ## which create numerous extra poles and zeros, most of which cancel.
  70 ## Those which do not cancel have a "fill-in" effect, extending the
  71 ## shorter of the sets to have the same number of as the longer of the
  72 ## sets of poles and zeros (or at least split the difference in the case
  73 ## of the band pass filter).  There may be other opportunistic
  74 ## cancellations but I will not check for them.
  75 ##
  76 ## Also note that any pole on the unit circle or beyond will result in
  77 ## an unstable filter.  Because of cancellation, this will only happen
  78 ## if the number of poles is smaller than the number of zeros and the
  79 ## filter is high pass or band pass.  The analytic design methods all
  80 ## yield more poles than zeros, so this will not be a problem.
  81 ##
  82 ## References:
  83 ##
  84 ## Proakis & Manolakis (1992). Digital Signal Processing. New York:
  85 ## Macmillan Publishing Company.
  86
  87 function [Sz, Sp, Sg] = sftrans(Sz, Sp, Sg, W, stop)
  88
  89   if (nargin != 5)
  90     print_usage;
  91   end
  92
  93   C = 1;
  94   p = length(Sp);
  95   z = length(Sz);
  96   if z > p || p == 0
  97     error("sftrans: must have at least as many poles as zeros in s-plane");
  98   end
  99
 100   if length(W)==2
 101     Fl = W(1);
 102     Fh = W(2);
 103     if stop
 104 ## ----------------  -------------------------  ------------------------
 105 ## Band Stop         zero: b ± sqrt(b^2-FhFl)   pole: b ± sqrt(b^2-FhFl)
 106 ##        S(Fh-Fl)   pole: ±sqrt(-FhFl)         zero: ±sqrt(-FhFl)
 107 ## S -> C --------   gain: -x                   gain: -1/x
 108 ##        S^2+FhFl   b=C/x (Fh-Fl)/2            b=C/x (Fh-Fl)/2
 109 ## ----------------  -------------------------  ------------------------
 110       if (isempty(Sz))
 111         Sg = Sg * real (1./ prod(-Sp));
 112       elseif (isempty(Sp))
 113         Sg = Sg * real(prod(-Sz));
 114       else
 115         Sg = Sg * real(prod(-Sz)/prod(-Sp));
 116       endif
 117       b = (C*(Fh-Fl)/2)./Sp;
 118       Sp = [b+sqrt(b.^2-Fh*Fl), b-sqrt(b.^2-Fh*Fl)];
 119       extend = [sqrt(-Fh*Fl), -sqrt(-Fh*Fl)];
 120       if isempty(Sz)
 121         Sz = [extend(1+rem([1:2*p],2))];
 122       else
 123         b = (C*(Fh-Fl)/2)./Sz;
 124         Sz = [b+sqrt(b.^2-Fh*Fl), b-sqrt(b.^2-Fh*Fl)];
 125         if (p > z)
 126           Sz = [Sz, extend(1+rem([1:2*(p-z)],2))];
 127         end
 128       end
 129     else
 130 ## ----------------  -------------------------  ------------------------
 131 ## Band Pass         zero: b ± sqrt(b^2-FhFl)   pole: b ± sqrt(b^2-FhFl)
 132 ##        S^2+FhFl   pole: 0                    zero: 0
 133 ## S -> C --------   gain: C/(Fh-Fl)            gain: (Fh-Fl)/C
 134 ##        S(Fh-Fl)   b=x/C (Fh-Fl)/2            b=x/C (Fh-Fl)/2
 135 ## ----------------  -------------------------  ------------------------
 136       Sg = Sg * (C/(Fh-Fl))^(z-p);
 137       b = Sp*((Fh-Fl)/(2*C));
 138       Sp = [b+sqrt(b.^2-Fh*Fl), b-sqrt(b.^2-Fh*Fl)];
 139       if isempty(Sz)
 140         Sz = zeros(1,p);
 141       else
 142         b = Sz*((Fh-Fl)/(2*C));
 143         Sz = [b+sqrt(b.^2-Fh*Fl), b-sqrt(b.^2-Fh*Fl)];
 144         if (p>z)
 145           Sz = [Sz, zeros(1, (p-z))];
 146         end
 147       end
 148     end
 149   else
 150     Fc = W;
 151     if stop
 152 ## ----------------  -------------------------  ------------------------
 153 ## High Pass         zero: Fc C/x               pole: Fc C/x
 154 ## S -> C Fc/S       pole: 0                    zero: 0
 155 ##                   gain: -x                   gain: -1/x
 156 ## ----------------  -------------------------  ------------------------
 157       if (isempty(Sz))
 158         Sg = Sg * real (1./ prod(-Sp));
 159       elseif (isempty(Sp))
 160         Sg = Sg * real(prod(-Sz));
 161       else
 162         Sg = Sg * real(prod(-Sz)/prod(-Sp));
 163       endif
 164       Sp = C * Fc ./ Sp;
 165       if isempty(Sz)
 166         Sz = zeros(1,p);
 167       else
 168         Sz = [C * Fc ./ Sz];
 169         if (p > z)
 170           Sz = [Sz, zeros(1,p-z)];
 171         end
 172       end
 173     else
 174 ## ----------------  -------------------------  ------------------------
 175 ## Low Pass          zero: Fc x/C               pole: Fc x/C
 176 ## S -> C S/Fc       gain: C/Fc                 gain: Fc/C
 177 ## ----------------  -------------------------  ------------------------
 178       Sg = Sg * (C/Fc)^(z-p);
 179       Sp = Fc * Sp / C;
 180       Sz = Fc * Sz / C;
 181     end
 182   end
 183 endfunction