filtfilt giving unexpected results

Question

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()

This is the eeg data

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?

Answer 1

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:

but it's supposed to look like this:

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)