You must log in to use subscribtions.

Clear message

This cookbook example shows how to design and use a low-pass FIR filter using functions from scipy.signal.

The pylab module from matplotlib is used to create plots.

   1 
   2 from numpy import cos, sin, pi, absolute, arange
   3 from scipy.signal import kaiserord, lfilter, firwin, freqz
   4 from pylab import figure, clf, plot, xlabel, ylabel, xlim, ylim, title, grid, axes, show
   5 
   6 
   7 #------------------------------------------------
   8 # Create a signal for demonstration.
   9 #------------------------------------------------
  10 
  11 sample_rate = 100.0
  12 nsamples = 400
  13 t = arange(nsamples) / sample_rate
  14 x = cos(2*pi*0.5*t) + 0.2*sin(2*pi*2.5*t+0.1) + \
        0.2*sin(2*pi*15.3*t) + 0.1*sin(2*pi*16.7*t + 0.1) + \
            0.1*sin(2*pi*23.45*t+.8)
  15 
  16 
  17 #------------------------------------------------
  18 # Create a FIR filter and apply it to x.
  19 #------------------------------------------------
  20 
  21 # The Nyquist rate of the signal.
  22 nyq_rate = sample_rate / 2.0
  23 
  24 # The desired width of the transition from pass to stop,
  25 # relative to the Nyquist rate.  We'll design the filter
  26 # with a 5 Hz transition width.
  27 width = 5.0/nyq_rate
  28 
  29 # The desired attenuation in the stop band, in dB.
  30 ripple_db = 60.0
  31 
  32 # Compute the order and Kaiser parameter for the FIR filter.
  33 N, beta = kaiserord(ripple_db, width)
  34 
  35 # The cutoff frequency of the filter.
  36 cutoff_hz = 10.0
  37 
  38 # Use firwin with a Kaiser window to create a lowpass FIR filter.
  39 taps = firwin(N, cutoff_hz/nyq_rate, window=('kaiser', beta))
  40 
  41 # Use lfilter to filter x with the FIR filter.
  42 filtered_x = lfilter(taps, 1.0, x)
  43 
  44 #------------------------------------------------
  45 # Plot the FIR filter coefficients.
  46 #------------------------------------------------
  47 
  48 figure(1)
  49 plot(taps, 'bo-', linewidth=2)
  50 title('Filter Coefficients (%d taps)' % N)
  51 grid(True)
  52 
  53 #------------------------------------------------
  54 # Plot the magnitude response of the filter.
  55 #------------------------------------------------
  56 
  57 figure(2)
  58 clf()
  59 w, h = freqz(taps, worN=8000)
  60 plot((w/pi)*nyq_rate, absolute(h), linewidth=2)
  61 xlabel('Frequency (Hz)')
  62 ylabel('Gain')
  63 title('Frequency Response')
  64 ylim(-0.05, 1.05)
  65 grid(True)
  66 
  67 # Upper inset plot.
  68 ax1 = axes([0.42, 0.6, .45, .25])
  69 plot((w/pi)*nyq_rate, absolute(h), linewidth=2)
  70 xlim(0,8.0)
  71 ylim(0.9985, 1.001)
  72 grid(True)
  73 
  74 # Lower inset plot
  75 ax2 = axes([0.42, 0.25, .45, .25])
  76 plot((w/pi)*nyq_rate, absolute(h), linewidth=2)
  77 xlim(12.0, 20.0)
  78 ylim(0.0, 0.0025)
  79 grid(True)
  80 
  81 #------------------------------------------------
  82 # Plot the original and filtered signals.
  83 #------------------------------------------------
  84 
  85 # The phase delay of the filtered signal.
  86 delay = 0.5 * (N-1) / sample_rate
  87 
  88 figure(3)
  89 # Plot the original signal.
  90 plot(t, x)
  91 # Plot the filtered signal, shifted to compensate for the phase delay.
  92 plot(t-delay, filtered_x, 'r-')
  93 # Plot just the "good" part of the filtered signal.  The first N-1
  94 # samples are "corrupted" by the initial conditions.
  95 plot(t[N-1:]-delay, filtered_x[N-1:], 'g', linewidth=4)
  96 
  97 xlabel('t')
  98 grid(True)
  99 
 100 show()

The plots generated by the above code:

The final plots shows the original signal (thin blue line), the filtered signal (shifted by the appropriate phase delay to align with the original signal; thin red line), and the "good" part of the filtered signal (heavy green line). The "good part" is the part of the signal that is not affected by the initial conditions.

Cookbook/FIRFilter (last edited 2011-06-06 04:57:12 by Warren Weckesser)