This is an archival dump of old wiki content --- see scipy.org for current material.

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:

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)