fir filter python tutorial

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