This cookbook recipe demonstrates the use of scipy.signal.butter to create a bandpass Butterworth filter. scipy.signal.freqz is used to compute the frequency response, and scipy.signal.lfilter is used to apply the filter to a signal. (This code was originally given in an answer to a question at stackoverflow.com.)
1 from scipy.signal import butter, lfilter
2
3
4 def butter_bandpass(lowcut, highcut, fs, order=5):
5 nyq = 0.5 * fs
6 low = lowcut / nyq
7 high = highcut / nyq
8 b, a = butter(order, [low, high], btype='band')
9 return b, a
10
11
12 def butter_bandpass_filter(data, lowcut, highcut, fs, order=5):
13 b, a = butter_bandpass(lowcut, highcut, fs, order=order)
14 y = lfilter(b, a, data)
15 return y
16
17
18 if __name__ == "__main__":
19 import numpy as np
20 import matplotlib.pyplot as plt
21 from scipy.signal import freqz
22
23 # Sample rate and desired cutoff frequencies (in Hz).
24 fs = 5000.0
25 lowcut = 500.0
26 highcut = 1250.0
27
28 # Plot the frequency response for a few different orders.
29 plt.figure(1)
30 plt.clf()
31 for order in [3, 6, 9]:
32 b, a = butter_bandpass(lowcut, highcut, fs, order=order)
33 w, h = freqz(b, a, worN=2000)
34 plt.plot((fs * 0.5 / np.pi) * w, abs(h), label="order = %d" % order)
35
36 plt.plot([0, 0.5 * fs], [np.sqrt(0.5), np.sqrt(0.5)],
37 '--', label='sqrt(0.5)')
38 plt.xlabel('Frequency (Hz)')
39 plt.ylabel('Gain')
40 plt.grid(True)
41 plt.legend(loc='best')
42
43 # Filter a noisy signal.
44 T = 0.05
45 nsamples = T * fs
46 t = np.linspace(0, T, nsamples, endpoint=False)
47 a = 0.02
48 f0 = 600.0
49 x = 0.1 * np.sin(2 * np.pi * 1.2 * np.sqrt(t))
50 x += 0.01 * np.cos(2 * np.pi * 312 * t + 0.1)
51 x += a * np.cos(2 * np.pi * f0 * t + .11)
52 x += 0.03 * np.cos(2 * np.pi * 2000 * t)
53 plt.figure(2)
54 plt.clf()
55 plt.plot(t, x, label='Noisy signal')
56
57 y = butter_bandpass_filter(x, lowcut, highcut, fs, order=6)
58 plt.plot(t, y, label='Filtered signal (%g Hz)' % f0)
59 plt.xlabel('time (seconds)')
60 plt.hlines([-a, a], 0, T, linestyles='--')
61 plt.grid(True)
62 plt.axis('tight')
63 plt.legend(loc='upper left')
64
65 plt.show()
The plots generated by the above code:
