107

UPDATE:

I found a Scipy Recipe based in this question! So, for anyone interested, go straight to: Contents » Signal processing » Butterworth Bandpass

I'm having a hard time to achieve what seemed initially a simple task of implementing a Butterworth band-pass filter for 1-D numpy array (time-series).

The parameters I have to include are the sample_rate, cutoff frequencies IN HERTZ and possibly order (other parameters, like attenuation, natural frequency, etc. are more obscure to me, so any "default" value would do).

What I have now is this, which seems to work as a high-pass filter but I'm no way sure if I'm doing it right:

def butter_highpass(interval, sampling_rate, cutoff, order=5):
    nyq = sampling_rate * 0.5

    stopfreq = float(cutoff)
    cornerfreq = 0.4 * stopfreq  # (?)

    ws = cornerfreq/nyq
    wp = stopfreq/nyq

    # for bandpass:
    # wp = [0.2, 0.5], ws = [0.1, 0.6]

    N, wn = scipy.signal.buttord(wp, ws, 3, 16)   # (?)

    # for hardcoded order:
    # N = order

    b, a = scipy.signal.butter(N, wn, btype='high')   # should 'high' be here for bandpass?
    sf = scipy.signal.lfilter(b, a, interval)
    return sf

enter image description here

The docs and examples are confusing and obscure, but I'd like to implement the form presented in the commend marked as "for bandpass". The question marks in the comments show where I just copy-pasted some example without understanding what is happening.

I am no electrical engineering or scientist, just a medical equipment designer needing to perform some rather straightforward bandpass filtering on EMG signals.

Dale K
  • 25,246
  • 15
  • 42
  • 71
heltonbiker
  • 26,657
  • 28
  • 137
  • 252
  • I've tried something at dsp.stackexchange, but they focus too much (more than I can handle) in conceptual issues of engineering and not so much in using the scipy functions. – heltonbiker Aug 23 '12 at 14:11

3 Answers3

143

You could skip the use of buttord, and instead just pick an order for the filter and see if it meets your filtering criterion. To generate the filter coefficients for a bandpass filter, give butter() the filter order, the cutoff frequencies Wn=[lowcut, highcut], the sampling rate fs (expressed in the same units as the cutoff frequencies) and the band type btype="band".

Here's a script that defines a couple convenience functions for working with a Butterworth bandpass filter. When run as a script, it makes two plots. One shows the frequency response at several filter orders for the same sampling rate and cutoff frequencies. The other plot demonstrates the effect of the filter (with order=6) on a sample time series.

from scipy.signal import butter, lfilter

def butter_bandpass(lowcut, highcut, fs, order=5):
    return butter(order, [lowcut, highcut], fs=fs, btype='band')

def butter_bandpass_filter(data, lowcut, highcut, fs, order=5):
    b, a = butter_bandpass(lowcut, highcut, fs, order=order)
    y = lfilter(b, a, data)
    return y


if __name__ == "__main__":
    import numpy as np
    import matplotlib.pyplot as plt
    from scipy.signal import freqz

    # Sample rate and desired cutoff frequencies (in Hz).
    fs = 5000.0
    lowcut = 500.0
    highcut = 1250.0

    # Plot the frequency response for a few different orders.
    plt.figure(1)
    plt.clf()
    for order in [3, 6, 9]:
        b, a = butter_bandpass(lowcut, highcut, fs, order=order)
        w, h = freqz(b, a, fs=fs, worN=2000)
        plt.plot(w, abs(h), label="order = %d" % order)

    plt.plot([0, 0.5 * fs], [np.sqrt(0.5), np.sqrt(0.5)],
             '--', label='sqrt(0.5)')
    plt.xlabel('Frequency (Hz)')
    plt.ylabel('Gain')
    plt.grid(True)
    plt.legend(loc='best')

    # Filter a noisy signal.
    T = 0.05
    nsamples = T * fs
    t = np.arange(0, nsamples) / fs
    a = 0.02
    f0 = 600.0
    x = 0.1 * np.sin(2 * np.pi * 1.2 * np.sqrt(t))
    x += 0.01 * np.cos(2 * np.pi * 312 * t + 0.1)
    x += a * np.cos(2 * np.pi * f0 * t + .11)
    x += 0.03 * np.cos(2 * np.pi * 2000 * t)
    plt.figure(2)
    plt.clf()
    plt.plot(t, x, label='Noisy signal')

    y = butter_bandpass_filter(x, lowcut, highcut, fs, order=6)
    plt.plot(t, y, label='Filtered signal (%g Hz)' % f0)
    plt.xlabel('time (seconds)')
    plt.hlines([-a, a], 0, T, linestyles='--')
    plt.grid(True)
    plt.axis('tight')
    plt.legend(loc='upper left')

    plt.show()

Here are the plots that are generated by this script:

Frequency response for several filter orders

enter image description here

Warren Weckesser
  • 110,654
  • 19
  • 194
  • 214
  • 1
    Do you know why the filtered output always starts at value zero? Is it possible to match it with the actual input value `x[0]`? I tried similar stuff with Cheby1 low pass filter, and I got the same problem. – LWZ Apr 03 '13 at 00:10
  • 2
    @LWZ: Use the function `scipy.signal.lfilter_zi`, and the `zi` argument to `lfilter`. For details, see the docstring for `lfilter_zi`. TL;DR? Just change `y = lfilter(b, a, data)` to `zi = lfilter_zi(b, a); y, zo = lfilter(b, a, data, zi=zi*data[0])`. (But this might not make a difference with a bandpass or high pass filter.) – Warren Weckesser Apr 03 '13 at 03:32
  • 1
    I noticed that there is 180 degree phase shift in the output of `scipy.signal.lfiter()` wrt to the original signal and the `signal.filtfilt()` output, why is that? Should I use `filtfilt()` instead if timing is important to me? – Jason Apr 21 '17 at 02:43
  • 1
    That's the phase delay of the filter at that frequency. The phase delay a sinusoid through a Butterworth filter depends nonlinearly on the frequency. For zero phase delay, yes, you can used `filtfilt()`. My answer [here](http://stackoverflow.com/questions/28536191/how-to-filter-smooth-with-scipy-numpy/28541805#28541805) includes an example of using `filtfilt()` to avoid a lag induced by the filter. – Warren Weckesser Apr 21 '17 at 03:52
  • @WarrenWeckesser Sorry I'm still bit lost, the signal processing field is quite new to me. I tried using `filtfilt()` instead, and plotted the power spectrum, and it seems the result is not a bandpass any more. Am I missing something fundamental here? – Jason Apr 21 '17 at 06:46
  • 1
    Hey Jason, I recommend asking questions about signal processing theory over at https://dsp.stackexchange.com/. If you have a question about some code you wrote that isn't working as expected, you can start a new question here on stackoverflow. – Warren Weckesser Apr 21 '17 at 12:21
  • @WarrenWeckesser Thanks for the reply, I think I found what's causing me the problems, the period I'm sampling is too short, and if I increase the time length the response functions show up clearer and the band pass filters work as expected. Again `filtfilt()` has no phase delay which is what I'm after. – Jason Apr 24 '17 at 02:28
  • Can I use this on an audio? I want to filter out hold-music ringtone from an audio, so that final audio will have just audio conversation and nothing else. – kRazzy R Dec 13 '17 at 14:33
  • I think it was better to leave the input variables out of your function definitions, so it would work after inputing any value – FabioSpaghetti Aug 08 '19 at 08:02
  • @WarrenWeckesser I need your help with the code, what does f0 indicate? and why do we add three functions with different frequencies , like f0 to x ? after the line you defined f0 – FabioSpaghetti Aug 09 '19 at 08:37
  • @FabioSpaghetti, the code where `f0` is used is only there to generate some sample noisy data to be filtered. `f0` is 600; it is the frequency of the sample data that I want to retain after filtering. The other values added to `x` are frequencies that I want the filter to eliminate. If you are trying to apply the filtering code in this answer to a signal of your own, you can ignore this whole section in which `x` is generated. Instead of generating `x` like I do here, you will set `x` to be your input signal. – Warren Weckesser Aug 09 '19 at 13:36
  • @WarrenWeckesser Thank you so much ! well I did but It is so hard to maintain the original shape of the signal while filtering only very high frequency peaks – FabioSpaghetti Aug 09 '19 at 13:52
  • @WarrenWeckesser NameError: name 'lfilter_zi' is not defined – FabioSpaghetti Oct 31 '19 at 17:31
  • @WarrenWeckesser if the lfilter_zi doesn't work for high pass filter, what should I do to avoid zeroing the average? – FabioSpaghetti Oct 31 '19 at 18:11
  • any tip for this? https://stackoverflow.com/q/60866193/5025009 – seralouk Mar 26 '20 at 12:47
63

The filter design method in accepted answer is correct, but it has a flaw. SciPy bandpass filters designed with b, a are unstable and may result in erroneous filters at higher filter orders.

Instead, use sos (second-order sections) output of filter design.

from scipy.signal import butter, sosfilt, sosfreqz

def butter_bandpass(lowcut, highcut, fs, order=5):
        nyq = 0.5 * fs
        low = lowcut / nyq
        high = highcut / nyq
        sos = butter(order, [low, high], analog=False, btype='band', output='sos')
        return sos

def butter_bandpass_filter(data, lowcut, highcut, fs, order=5):
        sos = butter_bandpass(lowcut, highcut, fs, order=order)
        y = sosfilt(sos, data)
        return y

Also, you can plot frequency response by changing 

b, a = butter_bandpass(lowcut, highcut, fs, order=order)
w, h = freqz(b, a, worN=2000)

to

sos = butter_bandpass(lowcut, highcut, fs, order=order)
w, h = sosfreqz(sos, worN=2000)
Mike
  • 19,114
  • 12
  • 59
  • 91
user13107
  • 3,239
  • 4
  • 34
  • 54
  • 7
    +1 because this is now the better way to go in many cases. As in the comments on the accepted answer, it's also possible to eliminate the phase delay using forward-backward filtering. Just replace `sosfilt` with `sosfiltfilt`. – Mike Mar 27 '18 at 02:10
  • @Mike and user13107 Does the same bug affect high-pass and low-pass Butterworth filters as well? And is the solution the same? – dewarrn1 Jan 25 '19 at 22:41
  • 8
    @dewarrn1 It's not actually correct to call it a "bug"; the algorithm is correctly implemented, but it's inherently unstable so it's just a bad choice of algorithm. But yes, it affects *any* filter at higher order — not just high- or low-pass and not just Butterworth filters, but also others [like Chebyshev](https://www.mathworks.com/help/signal/ref/cheby1.html#bucqqye-1) and so on. Anyway, in general, it's best to just always choose `sos` output, because that will always avoid the instability. And unless you need real-time processing, you should just always use `sosfiltfilt`. – Mike Jan 26 '19 at 13:53
  • 1
    Sorry I hadn't noticed this answer long ago! @user13107, yes, the transfer function (or 'ba') representation of a linear filter has some serious numerical issues when the order of filter is large. Even relatively low order filters can have problems when the desired bandwidth is small compared to the sampling frequency. My original answer was written before the SOS representation was added to SciPy, and before the `fs` argument was added to many of the functions in `scipy.signal`. The answer is long overdue for an update. – Warren Weckesser Aug 09 '19 at 13:43
  • any help for this? https://stackoverflow.com/q/60866193/5025009 – seralouk Mar 26 '20 at 12:47
  • @Mike if sosfiltfilt can remove phase delay, why not use it instead of sosfilt at all time? Why is real-time processing an exception? – John Jul 14 '22 at 03:03
  • @John Both `filtfilt` and `sosfiltfilt` filters the data forwards *and* backwards in time. This means that you have to have the complete signal in order to use `*filtfilt`. Put another way, if you want to know what the filtered value of your signal is at the moment you get the unfiltered value, `*filtfilt` needs data from the future. – Mike Jul 14 '22 at 13:25
  • @Mike but how do you know the length of the data from the future? Say we have two multi-sine signals (same freq. components) with different duration and we do offline filtering using *filtfilt, even if we use the same filter, we will get different filtered data due to different duration? – John Jul 14 '22 at 14:07
  • @John I don't think I understand this comment. Generally, yes you can expect to get different results for different signals, even if they represent identical underlying functions — regardless of what type of filter you use. The point is that for real-time filtering, you *can't* know the future, so you can't use `*filtfilt`. For post-processing, you just need the future *of a particular point*, which you're assumed to already have. – Mike Jul 15 '22 at 16:12
  • @Mike For post-processing using *filtfilt, to get the value of filtered signal at point N, I only need the value of input signal at point N + 1? – John Jul 18 '22 at 00:54
  • No, you basically need the complete signal. You *could* use `*filtfilt` as soon as you get point N, but when you include point N+1, the filtered point N would change — and again when you include point N+2, and so on. The changes should get smaller as you add more points, but there will still be some changes for a very long time, depending on the filter used. – Mike Jul 18 '22 at 09:51
4

For a bandpass filter, ws is a tuple containing the lower and upper corner frequencies. These represent the digital frequency where the filter response is 3 dB less than the passband.

wp is a tuple containing the stop band digital frequencies. They represent the location where the maximum attenuation begins.

gpass is the maximum attenutation in the passband in dB while gstop is the attentuation in the stopbands.

Say, for example, you wanted to design a filter for a sampling rate of 8000 samples/sec having corner frequencies of 300 and 3100 Hz. The Nyquist frequency is the sample rate divided by two, or in this example, 4000 Hz. The equivalent digital frequency is 1.0. The two corner frequencies are then 300/4000 and 3100/4000.

Now lets say you wanted the stopbands to be down 30 dB +/- 100 Hz from the corner frequencies. Thus, your stopbands would start at 200 and 3200 Hz resulting in the digital frequencies of 200/4000 and 3200/4000.

To create your filter, you'd call buttord as

fs = 8000.0
fso2 = fs/2
N,wn = scipy.signal.buttord(ws=[300/fso2,3100/fso2], wp=[200/fs02,3200/fs02],
   gpass=0.0, gstop=30.0)

The length of the resulting filter will be dependent upon the depth of the stop bands and the steepness of the response curve which is determined by the difference between the corner frequency and stopband frequency.

sizzzzlerz
  • 4,277
  • 3
  • 27
  • 35
  • I tried to implement it, but something is still missing. One thing is that `gpass=0.0` raises a division by zero error, so I changed it to 0.1 and the error stopped. Besides that, the docs for `butter` say: `Passband and stopband edge frequencies, normalized from 0 to 1 (1 corresponds to pi radians / sample).` I'm in doubt if your answer did the calculations right, so I'm still working on that and will give some feedback soon. – heltonbiker Aug 23 '12 at 16:53
  • (also, although my `ws` and `wp` have two elements each, the filter only performs low or high pass (via `btype` argument), but not band-pass) – heltonbiker Aug 23 '12 at 17:53
  • 1
    According to the documentation at http://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.buttord.html, buttord designs low, high, and band pass filters. As far as gpass, I guess buttord doesn't allow 0 dB attenuation in the passband. Set it to some non-zero value then. – sizzzzlerz Aug 23 '12 at 18:10