3
$\begingroup$

I am trying to filter EEG signals using butterworth filter and filtfilt. I have gone through a lot of documentation and these 2 commands seem sufficient for filtering. However, the results are bizarre.

from scipy.signal import butter, filtfilt
import sys, pickle
from numpy import *
import matplotlib.pyplot as plt

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

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

with open(sys.argv[1]) as eeg:
    eeg_data = pickle.load(eeg)
    eeg_data = eeg_data[:,3]

fs = 128
lowcut = 1
highcut = 40.0

y = butter_bandpass_filter(eeg_data, lowcut, highcut, fs, order = 9)
plt.plot(y,'r')
plt.show()

enter image description here

This is the eeg data enter image description here

Also, isn't the output in time domain? It seems incorrect irrespective of the domain. The order is too high i.e. 1e28 for the first 1000 points. I checked for the other points and the order is 1e17, even for points after 4000. Is this what I should expect?

$\endgroup$
  • 1
    $\begingroup$ Have a look at the filter coefficients $a$ and $b$. Maybe you could post them so we can see if there's anything wrong with the filter. $\endgroup$ – Matt L. Jul 7 '14 at 7:16
  • $\begingroup$ BTW, I don't know ScyPi but are you sure that "order=order" is OK to use (in butter_bandpass)? Here things might go very wrong. $\endgroup$ – Matt L. Jul 7 '14 at 7:20
  • $\begingroup$ Typically, you would like to have an even order for bandpass filters (when they are based on straightforward transformations as in this case). However, I am not sure how scipy deals with this, so that you get an 18:th-order filter when inserting 9 (there should be some warning if odd numbers doesn't work, so probably you get an 18:th-order filter). Now, assuming that you have an 18:th-order IIR filter with a high concentration of poles close to DC, it is not completely unlikely that you run into some sort of start-up behavior. Plot from 2000 samples. $\endgroup$ – Oscar Jul 7 '14 at 9:13
  • $\begingroup$ @MattL. order=order is fine in python, the first is the argument name, and the second is a variable. @user2497484 My guess would be that 9th order is too high and you're running into numerical problems. Can you try lower orders first? $\endgroup$ – endolith Jul 7 '14 at 14:04
  • $\begingroup$ Why does your eeg_data have a DC/mean value around 2350? Older scipy versions will produce inaccurate coefficients for higher orders, but I tried your script with eeg_data = randn(10000)*20+2350 and it worked ok. Does this work for you also? Maybe something is wrong with the data itself? $\endgroup$ – endolith Jul 7 '14 at 14:36
8
$\begingroup$

So the issue is that your filter order is too high. There are 2 problems with this:

  1. SciPy has a bug that generates inaccurate filters at high orders.
  2. On any platform, higher-order filters cannot be done in a single stage.

SciPy bug:

SciPy was previously generating prototype filters as a list of poles and zeros, then converting them to transfer functions, then transforming them to their new frequency and band-type, then converting them back into poles and zeros. This is bad because transfer function representation has numerical problems and loses precision. So your filter currently has these poles and zeros:

Z plane pole-zero plotcloseup of 1+0j

but it's supposed to look like this:

Z plane pole-zero plotcloseup of 1+0j

This probably isn't your problem, since these poles and zeros are still approximately the bandpass filter you asked for. At even higher orders the poles can be so wrong they jump outside the unit circle and make the filter explode. This has been fixed in SciPy 0.14.0, as long as you use output='zpk'.

Second-order sections:

When you actually do the filtering, you have to use a transfer function, not a list of poles and zeros, so your high-order filter would suffer the same numerical corruption if you tried to do the entire filter in a single stage. What we do instead is (starting from the list of poles and zeros, not from the inaccurate transfer function) break up the filter into second-order sections, which each have minimal numerical problems, and cascade the sections together to get the final equivalent output.

This is true for any platform, hardware or software, but Matlab and Octave have functions like sosfilt to do the work for you. SciPy does not, but it will in the future.

For now, you can break up the sections manually and use lfilter repeatedly with each section to get the desired output or copy the code from that pull request.

Update

SciPy has SOS now. So the best general-purpose filtering method is to generate your filter in SOS format, something like:

sos = butter(order, [low, high], btype='band', output='sos')

and then filter using sosfiltfilt:

y = sosfiltfilt(sos, data)
$\endgroup$
  • 1
    $\begingroup$ Thank you for your insight on why the higher order wasn't working. In my case, the combination of 'filtfilt' with order = 6 worked the best. 'lfilter' however, was giving trouble even with lower orders. $\endgroup$ – user2497484 Jul 11 '14 at 9:22
  • 1
    $\begingroup$ As an update the sos function is now available in scipy. I ended to use sosfiltfilt because sosfilt created a lot of artifacts at the beginning of the signal. $\endgroup$ – hadrienj Aug 4 '17 at 9:13
  • $\begingroup$ @hhh you understand the difference between sosfilt and sosfiltfilt though, right? $\endgroup$ – endolith Aug 4 '17 at 13:37

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.