2 # Copyright 2004,2005 Free Software Foundation, Inc.
4 # This file is part of GNU Radio
6 # GNU Radio is free software; you can redistribute it and/or modify
7 # it under the terms of the GNU General Public License as published by
8 # the Free Software Foundation; either version 3, or (at your option)
11 # GNU Radio is distributed in the hope that it will be useful,
12 # but WITHOUT ANY WARRANTY; without even the implied warranty of
13 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 # GNU General Public License for more details.
16 # You should have received a copy of the GNU General Public License
17 # along with GNU Radio; see the file COPYING. If not, write to
18 # the Free Software Foundation, Inc., 51 Franklin Street,
19 # Boston, MA 02110-1301, USA.
23 Routines for designing optimal FIR filters.
25 For a great intro to how all this stuff works, see section 6.6 of
26 "Digital Signal Processing: A Practical Approach", Emmanuael C. Ifeachor
27 and Barrie W. Jervis, Adison-Wesley, 1993. ISBN 0-201-54413-X.
31 from gnuradio import gr
35 # ----------------------------------------------------------------
37 def low_pass (gain, Fs, freq1, freq2, passband_ripple_db, stopband_atten_db,
39 passband_dev = passband_ripple_to_dev (passband_ripple_db)
40 stopband_dev = stopband_atten_to_dev (stopband_atten_db)
41 desired_ampls = (gain, 0)
42 (n, fo, ao, w) = remezord ([freq1, freq2], desired_ampls,
43 [passband_dev, stopband_dev], Fs)
44 taps = gr.remez (n + nextra_taps, fo, ao, w, "bandpass")
47 # FIXME high_passs is broken...
48 def high_pass (Fs, freq1, freq2, stopband_atten_db, passband_ripple_db,
51 passband_dev = passband_ripple_to_dev (passband_ripple_db)
52 stopband_dev = stopband_atten_to_dev (stopband_atten_db)
53 desired_ampls = (0, 1)
54 (n, fo, ao, w) = remezord ([freq1, freq2], desired_ampls,
55 [stopband_dev, passband_dev], Fs)
56 taps = gr.remez (n + nextra_taps, fo, ao, w, "bandpass")
59 # ----------------------------------------------------------------
61 def stopband_atten_to_dev (atten_db):
62 """Convert a stopband attenuation in dB to an absolute value"""
63 return 10**(-atten_db/20)
65 def passband_ripple_to_dev (ripple_db):
66 """Convert passband ripple spec expressed in dB to an absolute value"""
67 return (10**(ripple_db/20)-1)/(10**(ripple_db/20)+1)
69 # ----------------------------------------------------------------
71 def remezord (fcuts, mags, devs, fsamp = 2):
73 FIR order estimator (lowpass, highpass, bandpass, mulitiband).
75 (n, fo, ao, w) = remezord (f, a, dev)
76 (n, fo, ao, w) = remezord (f, a, dev, fs)
78 (n, fo, ao, w) = remezord (f, a, dev) finds the approximate order,
79 normalized frequency band edges, frequency band amplitudes, and
80 weights that meet input specifications f, a, and dev, to use with
83 * f is a sequence of frequency band edges (between 0 and Fs/2, where
84 Fs is the sampling frequency), and a is a sequence specifying the
85 desired amplitude on the bands defined by f. The length of f is
86 twice the length of a, minus 2. The desired function is
89 * dev is a sequence the same size as a that specifies the maximum
90 allowable deviation or ripples between the frequency response
91 and the desired amplitude of the output filter, for each band.
93 Use remez with the resulting order n, frequency sequence fo,
94 amplitude response sequence ao, and weights w to design the filter b
95 which approximately meets the specifications given by remezord
96 input parameters f, a, and dev:
98 b = remez (n, fo, ao, w)
100 (n, fo, ao, w) = remezord (f, a, dev, Fs) specifies a sampling frequency Fs.
102 Fs defaults to 2 Hz, implying a Nyquist frequency of 1 Hz. You can
103 therefore specify band edges scaled to a particular applications
106 In some cases remezord underestimates the order n. If the filter
107 does not meet the specifications, try a higher order such as n+1
115 for i in range (len (fcuts)):
116 fcuts[i] = float (fcuts[i]) / fsamp
124 raise ValueError, "Length of mags and devs must be equal"
126 if nf != 2 * (nbands - 1):
127 raise ValueError, "Length of f must be 2 * len (mags) - 2"
129 for i in range (len (mags)):
130 if mags[i] != 0: # if not stopband, get relative deviation
131 devs[i] = devs[i] / mags[i]
133 # separate the passband and stopband edges
139 for i in range (len (f1)):
140 if f2[i] - f1[i] < min_delta:
142 min_delta = f2[i] - f1[i]
145 # lowpass or highpass case (use formula)
146 l = lporder (f1[n], f2[n], devs[0], devs[1])
148 # bandpass or multipass case
149 # try different lowpasses and take the worst one that
150 # goes through the BP specs
152 for i in range (1, nbands-1):
153 l1 = lporder (f1[i-1], f2[i-1], devs[i], devs[i-1])
154 l2 = lporder (f1[i], f2[i], devs[i], devs[i+1])
157 n = int (math.ceil (l)) - 1 # need order, not length for remez
159 # cook up remez compatible result
160 ff = [0] + fcuts + [1]
161 for i in range (1, len (ff) - 1):
169 wts = [1] * len(devs)
170 for i in range (len (wts)):
171 wts[i] = max_dev / devs[i]
173 return (n, ff, aa, wts)
175 # ----------------------------------------------------------------
177 def lporder (freq1, freq2, delta_p, delta_s):
179 FIR lowpass filter length estimator. freq1 and freq2 are
180 normalized to the sampling frequency. delta_p is the passband
181 deviation (ripple), delta_s is the stopband deviation (ripple).
183 Note, this works for high pass filters too (freq1 > freq2), but
184 doesnt work well if the transition is near f == 0 or f == fs/2
186 From Herrmann et al (1973), Practical design rules for optimum
187 finite impulse response filters. Bell System Technical J., 52, 769-99
189 df = abs (freq2 - freq1)
190 ddp = math.log10 (delta_p)
191 dds = math.log10 (delta_s)
208 dinf=((t1 + t2 + a3) * dds) + (t3 + t4 + a6)
209 ff = b1 + b2 * (ddp - dds)
210 n = dinf / df - ff * df + 1
214 def bporder (freq1, freq2, delta_p, delta_s):
216 FIR bandpass filter length estimator. freq1 and freq2 are
217 normalized to the sampling frequency. delta_p is the passband
218 deviation (ripple), delta_s is the stopband deviation (ripple).
220 From Mintzer and Liu (1979)
222 df = abs (freq2 - freq1)
223 ddp = math.log10 (delta_p)
224 dds = math.log10 (delta_s)
238 cinf = dds * (t1 + t2 + a3) + t3 + t4 + a6
239 ginf = -14.6 * math.log10 (delta_p / delta_s) - 16.9
240 n = cinf / df + ginf * df + 1