1 ## Copyright (C) 1999-2001 Paul Kienzle <pkienzle@users.sf.net>
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
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
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/>.
16 ## usage: [Sz, Sp, Sg] = sftrans(Sz, Sp, Sg, W, stop)
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).
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:
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
53 ## ---------------- ------------------------- ------------------------
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:
62 ## H(S)=g prod(S-Xi)/prod(S-Xj)
64 ## The transforms are from the references. The actual pole-zero-gain
65 ## changes I derived myself.
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.
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.
84 ## Proakis & Manolakis (1992). Digital Signal Processing. New York:
85 ## Macmillan Publishing Company.
87 function [Sz, Sp, Sg] = sftrans(Sz, Sp, Sg, W, stop)
97 error("sftrans: must have at least as many poles as zeros in s-plane");
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 ## ---------------- ------------------------- ------------------------
111 Sg = Sg * real (1./ prod(-Sp));
113 Sg = Sg * real(prod(-Sz));
115 Sg = Sg * real(prod(-Sz)/prod(-Sp));
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)];
121 Sz = [extend(1+rem([1:2*p],2))];
123 b = (C*(Fh-Fl)/2)./Sz;
124 Sz = [b+sqrt(b.^2-Fh*Fl), b-sqrt(b.^2-Fh*Fl)];
126 Sz = [Sz, extend(1+rem([1:2*(p-z)],2))];
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)];
142 b = Sz*((Fh-Fl)/(2*C));
143 Sz = [b+sqrt(b.^2-Fh*Fl), b-sqrt(b.^2-Fh*Fl)];
145 Sz = [Sz, zeros(1, (p-z))];
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 ## ---------------- ------------------------- ------------------------
158 Sg = Sg * real (1./ prod(-Sp));
160 Sg = Sg * real(prod(-Sz));
162 Sg = Sg * real(prod(-Sz)/prod(-Sp));
170 Sz = [Sz, zeros(1,p-z)];
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);