2 # Copyright 2004,2005,2009 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 ## Builds a low pass filter.
38 # @param gain Filter gain in the passband (linear)
39 # @param Fs Sampling rate (sps)
40 # @param freq1 End of pass band (in Hz)
41 # @param freq2 Start of stop band (in Hz)
42 # @param passband_ripple_db Pass band ripple in dB (should be small, < 1)
43 # @param stopband_atten_db Stop band attenuation in dB (should be large, >= 60)
44 # @param nextra_taps Extra taps to use in the filter (default=2)
45 def low_pass (gain, Fs, freq1, freq2, passband_ripple_db, stopband_atten_db,
47 passband_dev = passband_ripple_to_dev (passband_ripple_db)
48 stopband_dev = stopband_atten_to_dev (stopband_atten_db)
49 desired_ampls = (gain, 0)
50 (n, fo, ao, w) = remezord ([freq1, freq2], desired_ampls,
51 [passband_dev, stopband_dev], Fs)
52 # The remezord typically under-estimates the filter order, so add 2 taps by default
53 taps = gr.remez (n + nextra_taps, fo, ao, w, "bandpass")
56 ## Builds a band pass filter.
57 # @param gain Filter gain in the passband (linear)
58 # @param Fs Sampling rate (sps)
59 # @param freq_sb1 End of stop band (in Hz)
60 # @param freq_pb1 Start of pass band (in Hz)
61 # @param freq_pb2 End of pass band (in Hz)
62 # @param freq_sb2 Start of stop band (in Hz)
63 # @param passband_ripple_db Pass band ripple in dB (should be small, < 1)
64 # @param stopband_atten_db Stop band attenuation in dB (should be large, >= 60)
65 # @param nextra_taps Extra taps to use in the filter (default=2)
66 def band_pass (gain, Fs, freq_sb1, freq_pb1, freq_pb2, freq_sb2,
67 passband_ripple_db, stopband_atten_db,
69 passband_dev = passband_ripple_to_dev (passband_ripple_db)
70 stopband_dev = stopband_atten_to_dev (stopband_atten_db)
71 desired_ampls = (0, gain, 0)
72 desired_freqs = [freq_sb1, freq_pb1, freq_pb2, freq_sb2]
73 desired_ripple = [stopband_dev, passband_dev, stopband_dev]
74 (n, fo, ao, w) = remezord (desired_freqs, desired_ampls,
76 # The remezord typically under-estimates the filter order, so add 2 taps by default
77 taps = gr.remez (n + nextra_taps, fo, ao, w, "bandpass")
81 ## Builds a band pass filter with complex taps by making an LPF and
82 # spinning it up to the right center frequency
83 # @param gain Filter gain in the passband (linear)
84 # @param Fs Sampling rate (sps)
85 # @param freq_sb1 End of stop band (in Hz)
86 # @param freq_pb1 Start of pass band (in Hz)
87 # @param freq_pb2 End of pass band (in Hz)
88 # @param freq_sb2 Start of stop band (in Hz)
89 # @param passband_ripple_db Pass band ripple in dB (should be small, < 1)
90 # @param stopband_atten_db Stop band attenuation in dB (should be large, >= 60)
91 # @param nextra_taps Extra taps to use in the filter (default=2)
92 def complex_band_pass (gain, Fs, freq_sb1, freq_pb1, freq_pb2, freq_sb2,
93 passband_ripple_db, stopband_atten_db,
95 center_freq = (freq_pb2 + freq_pb1) / 2.0
96 lp_pb = (freq_pb2 - center_freq)/1.0
97 lp_sb = freq_sb2 - center_freq
98 lptaps = low_pass(gain, Fs, lp_pb, lp_sb, passband_ripple_db,
99 stopband_atten_db, nextra_taps)
100 spinner = [cmath.exp(2j*cmath.pi*center_freq/Fs*i) for i in xrange(len(lptaps))]
101 taps = [s*t for s,t in zip(spinner, lptaps)]
105 ## Builds a band reject filter
106 # spinning it up to the right center frequency
107 # @param gain Filter gain in the passband (linear)
108 # @param Fs Sampling rate (sps)
109 # @param freq_pb1 End of pass band (in Hz)
110 # @param freq_sb1 Start of stop band (in Hz)
111 # @param freq_sb2 End of stop band (in Hz)
112 # @param freq_pb2 Start of pass band (in Hz)
113 # @param passband_ripple_db Pass band ripple in dB (should be small, < 1)
114 # @param stopband_atten_db Stop band attenuation in dB (should be large, >= 60)
115 # @param nextra_taps Extra taps to use in the filter (default=2)
116 def band_reject (gain, Fs, freq_pb1, freq_sb1, freq_sb2, freq_pb2,
117 passband_ripple_db, stopband_atten_db,
119 passband_dev = passband_ripple_to_dev (passband_ripple_db)
120 stopband_dev = stopband_atten_to_dev (stopband_atten_db)
121 desired_ampls = (gain, 0, gain)
122 desired_freqs = [freq_pb1, freq_sb1, freq_sb2, freq_pb2]
123 desired_ripple = [passband_dev, stopband_dev, passband_dev]
124 (n, fo, ao, w) = remezord (desired_freqs, desired_ampls,
126 # Make sure we use an odd number of taps
127 if((n+nextra_taps)%2 == 1):
129 # The remezord typically under-estimates the filter order, so add 2 taps by default
130 taps = gr.remez (n + nextra_taps, fo, ao, w, "bandpass")
134 ## Builds a high pass filter.
135 # @param gain Filter gain in the passband (linear)
136 # @param Fs Sampling rate (sps)
137 # @param freq1 End of stop band (in Hz)
138 # @param freq2 Start of pass band (in Hz)
139 # @param passband_ripple_db Pass band ripple in dB (should be small, < 1)
140 # @param stopband_atten_db Stop band attenuation in dB (should be large, >= 60)
141 # @param nextra_taps Extra taps to use in the filter (default=2)
142 def high_pass (gain, Fs, freq1, freq2, passband_ripple_db, stopband_atten_db,
144 passband_dev = passband_ripple_to_dev (passband_ripple_db)
145 stopband_dev = stopband_atten_to_dev (stopband_atten_db)
146 desired_ampls = (0, 1)
147 (n, fo, ao, w) = remezord ([freq1, freq2], desired_ampls,
148 [stopband_dev, passband_dev], Fs)
149 # For a HPF, we need to use an odd number of taps
150 # In gr.remez, ntaps = n+1, so n must be even
151 if((n+nextra_taps)%2 == 1):
154 # The remezord typically under-estimates the filter order, so add 2 taps by default
155 taps = gr.remez (n + nextra_taps, fo, ao, w, "bandpass")
158 # ----------------------------------------------------------------
160 def stopband_atten_to_dev (atten_db):
161 """Convert a stopband attenuation in dB to an absolute value"""
162 return 10**(-atten_db/20)
164 def passband_ripple_to_dev (ripple_db):
165 """Convert passband ripple spec expressed in dB to an absolute value"""
166 return (10**(ripple_db/20)-1)/(10**(ripple_db/20)+1)
168 # ----------------------------------------------------------------
170 def remezord (fcuts, mags, devs, fsamp = 2):
172 FIR order estimator (lowpass, highpass, bandpass, mulitiband).
174 (n, fo, ao, w) = remezord (f, a, dev)
175 (n, fo, ao, w) = remezord (f, a, dev, fs)
177 (n, fo, ao, w) = remezord (f, a, dev) finds the approximate order,
178 normalized frequency band edges, frequency band amplitudes, and
179 weights that meet input specifications f, a, and dev, to use with
182 * f is a sequence of frequency band edges (between 0 and Fs/2, where
183 Fs is the sampling frequency), and a is a sequence specifying the
184 desired amplitude on the bands defined by f. The length of f is
185 twice the length of a, minus 2. The desired function is
188 * dev is a sequence the same size as a that specifies the maximum
189 allowable deviation or ripples between the frequency response
190 and the desired amplitude of the output filter, for each band.
192 Use remez with the resulting order n, frequency sequence fo,
193 amplitude response sequence ao, and weights w to design the filter b
194 which approximately meets the specifications given by remezord
195 input parameters f, a, and dev:
197 b = remez (n, fo, ao, w)
199 (n, fo, ao, w) = remezord (f, a, dev, Fs) specifies a sampling frequency Fs.
201 Fs defaults to 2 Hz, implying a Nyquist frequency of 1 Hz. You can
202 therefore specify band edges scaled to a particular applications
205 In some cases remezord underestimates the order n. If the filter
206 does not meet the specifications, try a higher order such as n+1
214 for i in range (len (fcuts)):
215 fcuts[i] = float (fcuts[i]) / fsamp
223 raise ValueError, "Length of mags and devs must be equal"
225 if nf != 2 * (nbands - 1):
226 raise ValueError, "Length of f must be 2 * len (mags) - 2"
228 for i in range (len (mags)):
229 if mags[i] != 0: # if not stopband, get relative deviation
230 devs[i] = devs[i] / mags[i]
232 # separate the passband and stopband edges
238 for i in range (len (f1)):
239 if f2[i] - f1[i] < min_delta:
241 min_delta = f2[i] - f1[i]
244 # lowpass or highpass case (use formula)
245 l = lporder (f1[n], f2[n], devs[0], devs[1])
247 # bandpass or multipass case
248 # try different lowpasses and take the worst one that
249 # goes through the BP specs
251 for i in range (1, nbands-1):
252 l1 = lporder (f1[i-1], f2[i-1], devs[i], devs[i-1])
253 l2 = lporder (f1[i], f2[i], devs[i], devs[i+1])
256 n = int (math.ceil (l)) - 1 # need order, not length for remez
258 # cook up remez compatible result
259 ff = [0] + fcuts + [1]
260 for i in range (1, len (ff) - 1):
268 wts = [1] * len(devs)
269 for i in range (len (wts)):
270 wts[i] = max_dev / devs[i]
272 return (n, ff, aa, wts)
274 # ----------------------------------------------------------------
276 def lporder (freq1, freq2, delta_p, delta_s):
278 FIR lowpass filter length estimator. freq1 and freq2 are
279 normalized to the sampling frequency. delta_p is the passband
280 deviation (ripple), delta_s is the stopband deviation (ripple).
282 Note, this works for high pass filters too (freq1 > freq2), but
283 doesnt work well if the transition is near f == 0 or f == fs/2
285 From Herrmann et al (1973), Practical design rules for optimum
286 finite impulse response filters. Bell System Technical J., 52, 769-99
288 df = abs (freq2 - freq1)
289 ddp = math.log10 (delta_p)
290 dds = math.log10 (delta_s)
307 dinf=((t1 + t2 + a3) * dds) + (t3 + t4 + a6)
308 ff = b1 + b2 * (ddp - dds)
309 n = dinf / df - ff * df + 1
313 def bporder (freq1, freq2, delta_p, delta_s):
315 FIR bandpass filter length estimator. freq1 and freq2 are
316 normalized to the sampling frequency. delta_p is the passband
317 deviation (ripple), delta_s is the stopband deviation (ripple).
319 From Mintzer and Liu (1979)
321 df = abs (freq2 - freq1)
322 ddp = math.log10 (delta_p)
323 dds = math.log10 (delta_s)
337 cinf = dds * (t1 + t2 + a3) + t3 + t4 + a6
338 ginf = -14.6 * math.log10 (delta_p / delta_s) - 16.9
339 n = cinf / df + ginf * df + 1