1 ## Copyright (C) 1999 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: [Zz, Zp, Zg] = bilinear(Sz, Sp, Sg, T)
17 ## [Zb, Za] = bilinear(Sb, Sa, T)
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.
24 ## Note: this differs from the bilinear function in the signal processing
25 ## toolbox, which uses 1/T rather than T.
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.
35 ## The following table summarizes the transformation:
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 ## +---------------+-----------------------+----------------------+
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
51 ## H(Z)=g prod(Z-Xi)/prod(Z-Xj)
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.
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
70 ## Proakis & Manolakis (1992). Digital Signal Processing. New York:
71 ## Macmillan Publishing Company.
73 function [Zz, Zp, Zg] = bilinear(Sz, Sp, Sg, T)
77 [Sz, Sp, Sg] = tf2zp(Sz, Sp);
85 error("bilinear: must have at least as many poles as zeros in s-plane");
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
93 ## ---------------- ------------------------- ------------------------
94 Zg = real(Sg * prod((2-Sz*T)/T) / prod((2-Sp*T)/T));
95 Zp = (2+Sp*T)./(2-Sp*T);
99 Zz = [(2+Sz*T)./(2-Sz*T)];
100 Zz = postpad(Zz, p, -1);
103 if nargout==2, [Zz, Zp] = zp2tf(Zz, Zp, Zg); endif