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 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,
 118                  nextra_taps=2):
 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,
 125                                desired_ripple, Fs)
 126     # Make sure we use an odd number of taps
 127     if((n+nextra_taps)%2 == 1):
 128         n += 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")
 131     return taps
 132
 133
 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,
 143                nextra_taps=2):
 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):
 152         n += 1
 153
 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")
 156     return taps
 157
 158 # ----------------------------------------------------------------
 159
 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)
 163
 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)
 167
 168 # ----------------------------------------------------------------
 169
 170 def remezord (fcuts, mags, devs, fsamp = 2):
 171     '''
 172     FIR order estimator (lowpass, highpass, bandpass, mulitiband).
 173
 174     (n, fo, ao, w) = remezord (f, a, dev)
 175     (n, fo, ao, w) = remezord (f, a, dev, fs)
 176
 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
 180     the remez command.
 181
 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
 186       piecewise constant.
 187
 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.
 191
 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:
 196
 197     b = remez (n, fo, ao, w)
 198
 199     (n, fo, ao, w) = remezord (f, a, dev, Fs) specifies a sampling frequency Fs.
 200
 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
 203     sampling frequency.
 204
 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
 207     or n+2.
 208     '''
 209     # get local copies
 210     fcuts = fcuts[:]
 211     mags = mags[:]
 212     devs = devs[:]
 213
 214     for i in range (len (fcuts)):
 215         fcuts[i] = float (fcuts[i]) / fsamp
 216
 217     nf = len (fcuts)
 218     nm = len (mags)
 219     nd = len (devs)
 220     nbands = nm
 221
 222     if nm != nd:
 223         raise ValueError, "Length of mags and devs must be equal"
 224
 225     if nf != 2 * (nbands - 1):
 226         raise ValueError, "Length of f must be 2 * len (mags) - 2"
 227
 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]
 231
 232     # separate the passband and stopband edges
 233     f1 = fcuts[0::2]
 234     f2 = fcuts[1::2]
 235
 236     n = 0
 237     min_delta = 2
 238     for i in range (len (f1)):
 239         if f2[i] - f1[i] < min_delta:
 240             n = i
 241             min_delta = f2[i] - f1[i]
 242
 243     if nbands == 2:
 244         # lowpass or highpass case (use formula)
 245         l = lporder (f1[n], f2[n], devs[0], devs[1])
 246     else:
 247         # bandpass or multipass case
 248         # try different lowpasses and take the worst one that
 249         #  goes through the BP specs
 250         l = 0
 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])
 254             l = max (l, l1, l2)
 255
 256     n = int (math.ceil (l)) - 1               # need order, not length for remez
 257
 258     # cook up remez compatible result
 259     ff = [0] + fcuts + [1]
 260     for i in range (1, len (ff) - 1):
 261         ff[i] *= 2
 262
 263     aa = []
 264     for a in mags:
 265         aa = aa + [a, a]
 266
 267     max_dev = max (devs)
 268     wts = [1] * len(devs)
 269     for i in range (len (wts)):
 270         wts[i] = max_dev / devs[i]
 271
 272     return (n, ff, aa, wts)
 273
 274 # ----------------------------------------------------------------
 275
 276 def lporder (freq1, freq2, delta_p, delta_s):
 277     '''
 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).
 281
 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
 284
 285     From Herrmann et al (1973), Practical design rules for optimum
 286     finite impulse response filters.  Bell System Technical J., 52, 769-99
 287     '''
 288     df = abs (freq2 - freq1)
 289     ddp = math.log10 (delta_p)
 290     dds = math.log10 (delta_s)
 291
 292     a1 = 5.309e-3
 293     a2 = 7.114e-2
 294     a3 = -4.761e-1
 295     a4 = -2.66e-3
 296     a5 = -5.941e-1
 297     a6 = -4.278e-1
 298
 299     b1 = 11.01217
 300     b2 = 0.5124401
 301
 302     t1 = a1 * ddp * ddp
 303     t2 = a2 * ddp
 304     t3 = a4 * ddp * ddp
 305     t4 = a5 * ddp
 306
 307     dinf=((t1 + t2 + a3) * dds) + (t3 + t4 + a6)
 308     ff = b1 + b2 * (ddp - dds)
 309     n = dinf / df - ff * df + 1
 310     return n
 311
 312
 313 def bporder (freq1, freq2, delta_p, delta_s):
 314     '''
 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).
 318
 319     From Mintzer and Liu (1979)
 320     '''
 321     df = abs (freq2 - freq1)
 322     ddp = math.log10 (delta_p)
 323     dds = math.log10 (delta_s)
 324
 325     a1 = 0.01201
 326     a2 = 0.09664
 327     a3 = -0.51325
 328     a4 = 0.00203
 329     a5 = -0.57054
 330     a6 = -0.44314
 331
 332     t1 = a1 * ddp * ddp
 333     t2 = a2 * ddp
 334     t3 = a4 * ddp * ddp
 335     t4 = a5 * ddp
 336
 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
 340     return n
 341