This is an archival dump of old wiki content --- see scipy.org for current material.
Please see http://scipy-cookbook.readthedocs.org/

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

The plots generated by the above code:

fir_demo_taps.png fir_demo_freq_resp.png fir_demo_signals.png

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.

SciPy: Cookbook/FIRFilter (last edited 2015-10-24 17:48:23 by anonymous)