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