git.gag.com Git - debian/gnuradio/blob - gnuradio-core/src/python/gnuradio/optfir.py

   1 #
   2 # Copyright 2004,2005,2009 Free Software Foundation, Inc.
   3 #
   4 # This file is part of GNU Radio
   5 #
   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)
   9 # any later version.
  10 #
  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.
  15 #
  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.
  20 #
  21
  22 '''
  23 Routines for designing optimal FIR filters.
  24
  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.
  28 '''
  29
  30 import math, cmath
  31 from gnuradio import gr
  32
  33 remez = gr.remez
  34
  35 # ----------------------------------------------------------------
  36
  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,
  46               nextra_taps=2):
  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")
  54     return taps
  55
  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,
  68                nextra_taps=2):
  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,
  75                                desired_ripple, Fs)
  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")
  78     return taps
  79
  80
  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,
  94                        nextra_taps=2):
  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)]
 102     return taps
 103
 104
 105 ##  Builds a high pass filter.
 106 #   @param gain  Filter gain in the passband (linear)
 107 #   @param Fs    Sampling rate (sps)
 108 #   @param freq1 End of stop band (in Hz)
 109 #   @param freq2 Start of pass band (in Hz)
 110 #   @param passband_ripple_db Pass band ripple in dB (should be small, < 1)
 111 #   @param stopband_atten_db  Stop band attenuation in dB (should be large, >= 60)
 112 #   @param nextra_taps  Extra taps to use in the filter (default=2)
 113 def high_pass (gain, Fs, freq1, freq2, passband_ripple_db, stopband_atten_db,
 114                nextra_taps=2):
 115     passband_dev = passband_ripple_to_dev (passband_ripple_db)
 116     stopband_dev = stopband_atten_to_dev (stopband_atten_db)
 117     desired_ampls = (0, 1)
 118     (n, fo, ao, w) = remezord ([freq1, freq2], desired_ampls,
 119                                [stopband_dev, passband_dev], Fs)
 120     # For a HPF, we need to use an odd number of taps
 121     # In gr.remez, ntaps = n+1, so n must be even
 122     if((n+nextra_taps)%2 == 1):
 123         n += 1
 124
 125     # The remezord typically under-estimates the filter order, so add 2 taps by default
 126     taps = gr.remez (n + nextra_taps, fo, ao, w, "bandpass")
 127     return taps
 128
 129 # ----------------------------------------------------------------
 130
 131 def stopband_atten_to_dev (atten_db):
 132     """Convert a stopband attenuation in dB to an absolute value"""
 133     return 10**(-atten_db/20)
 134
 135 def passband_ripple_to_dev (ripple_db):
 136     """Convert passband ripple spec expressed in dB to an absolute value"""
 137     return (10**(ripple_db/20)-1)/(10**(ripple_db/20)+1)
 138
 139 # ----------------------------------------------------------------
 140
 141 def remezord (fcuts, mags, devs, fsamp = 2):
 142     '''
 143     FIR order estimator (lowpass, highpass, bandpass, mulitiband).
 144
 145     (n, fo, ao, w) = remezord (f, a, dev)
 146     (n, fo, ao, w) = remezord (f, a, dev, fs)
 147
 148     (n, fo, ao, w) = remezord (f, a, dev) finds the approximate order,
 149     normalized frequency band edges, frequency band amplitudes, and
 150     weights that meet input specifications f, a, and dev, to use with
 151     the remez command.
 152
 153     * f is a sequence of frequency band edges (between 0 and Fs/2, where
 154       Fs is the sampling frequency), and a is a sequence specifying the
 155       desired amplitude on the bands defined by f. The length of f is
 156       twice the length of a, minus 2. The desired function is
 157       piecewise constant.
 158
 159     * dev is a sequence the same size as a that specifies the maximum
 160       allowable deviation or ripples between the frequency response
 161       and the desired amplitude of the output filter, for each band.
 162
 163     Use remez with the resulting order n, frequency sequence fo,
 164     amplitude response sequence ao, and weights w to design the filter b
 165     which approximately meets the specifications given by remezord
 166     input parameters f, a, and dev:
 167
 168     b = remez (n, fo, ao, w)
 169
 170     (n, fo, ao, w) = remezord (f, a, dev, Fs) specifies a sampling frequency Fs.
 171
 172     Fs defaults to 2 Hz, implying a Nyquist frequency of 1 Hz. You can
 173     therefore specify band edges scaled to a particular applications
 174     sampling frequency.
 175
 176     In some cases remezord underestimates the order n. If the filter
 177     does not meet the specifications, try a higher order such as n+1
 178     or n+2.
 179     '''
 180     # get local copies
 181     fcuts = fcuts[:]
 182     mags = mags[:]
 183     devs = devs[:]
 184
 185     for i in range (len (fcuts)):
 186         fcuts[i] = float (fcuts[i]) / fsamp
 187
 188     nf = len (fcuts)
 189     nm = len (mags)
 190     nd = len (devs)
 191     nbands = nm
 192
 193     if nm != nd:
 194         raise ValueError, "Length of mags and devs must be equal"
 195
 196     if nf != 2 * (nbands - 1):
 197         raise ValueError, "Length of f must be 2 * len (mags) - 2"
 198
 199     for i in range (len (mags)):
 200         if mags[i] != 0:                        # if not stopband, get relative deviation
 201             devs[i] = devs[i] / mags[i]
 202
 203     # separate the passband and stopband edges
 204     f1 = fcuts[0::2]
 205     f2 = fcuts[1::2]
 206
 207     n = 0
 208     min_delta = 2
 209     for i in range (len (f1)):
 210         if f2[i] - f1[i] < min_delta:
 211             n = i
 212             min_delta = f2[i] - f1[i]
 213
 214     if nbands == 2:
 215         # lowpass or highpass case (use formula)
 216         l = lporder (f1[n], f2[n], devs[0], devs[1])
 217     else:
 218         # bandpass or multipass case
 219         # try different lowpasses and take the worst one that
 220         #  goes through the BP specs
 221         l = 0
 222         for i in range (1, nbands-1):
 223             l1 = lporder (f1[i-1], f2[i-1], devs[i], devs[i-1])
 224             l2 = lporder (f1[i], f2[i], devs[i], devs[i+1])
 225             l = max (l, l1, l2)
 226
 227     n = int (math.ceil (l)) - 1               # need order, not length for remez
 228
 229     # cook up remez compatible result
 230     ff = [0] + fcuts + [1]
 231     for i in range (1, len (ff) - 1):
 232         ff[i] *= 2
 233
 234     aa = []
 235     for a in mags:
 236         aa = aa + [a, a]
 237
 238     max_dev = max (devs)
 239     wts = [1] * len(devs)
 240     for i in range (len (wts)):
 241         wts[i] = max_dev / devs[i]
 242
 243     return (n, ff, aa, wts)
 244
 245 # ----------------------------------------------------------------
 246
 247 def lporder (freq1, freq2, delta_p, delta_s):
 248     '''
 249     FIR lowpass filter length estimator.  freq1 and freq2 are
 250     normalized to the sampling frequency.  delta_p is the passband
 251     deviation (ripple), delta_s is the stopband deviation (ripple).
 252
 253     Note, this works for high pass filters too (freq1 > freq2), but
 254     doesnt work well if the transition is near f == 0 or f == fs/2
 255
 256     From Herrmann et al (1973), Practical design rules for optimum
 257     finite impulse response filters.  Bell System Technical J., 52, 769-99
 258     '''
 259     df = abs (freq2 - freq1)
 260     ddp = math.log10 (delta_p)
 261     dds = math.log10 (delta_s)
 262
 263     a1 = 5.309e-3
 264     a2 = 7.114e-2
 265     a3 = -4.761e-1
 266     a4 = -2.66e-3
 267     a5 = -5.941e-1
 268     a6 = -4.278e-1
 269
 270     b1 = 11.01217
 271     b2 = 0.5124401
 272
 273     t1 = a1 * ddp * ddp
 274     t2 = a2 * ddp
 275     t3 = a4 * ddp * ddp
 276     t4 = a5 * ddp
 277
 278     dinf=((t1 + t2 + a3) * dds) + (t3 + t4 + a6)
 279     ff = b1 + b2 * (ddp - dds)
 280     n = dinf / df - ff * df + 1
 281     return n
 282
 283
 284 def bporder (freq1, freq2, delta_p, delta_s):
 285     '''
 286     FIR bandpass filter length estimator.  freq1 and freq2 are
 287     normalized to the sampling frequency.  delta_p is the passband
 288     deviation (ripple), delta_s is the stopband deviation (ripple).
 289
 290     From Mintzer and Liu (1979)
 291     '''
 292     df = abs (freq2 - freq1)
 293     ddp = math.log10 (delta_p)
 294     dds = math.log10 (delta_s)
 295
 296     a1 = 0.01201
 297     a2 = 0.09664
 298     a3 = -0.51325
 299     a4 = 0.00203
 300     a5 = -0.57054
 301     a6 = -0.44314
 302
 303     t1 = a1 * ddp * ddp
 304     t2 = a2 * ddp
 305     t3 = a4 * ddp * ddp
 306     t4 = a5 * ddp
 307
 308     cinf = dds * (t1 + t2 + a3) + t3 + t4 + a6
 309     ginf = -14.6 * math.log10 (delta_p / delta_s) - 16.9
 310     n = cinf / df + ginf * df + 1
 311     return n
 312