Differences between Python and MATLAB filtfilt function

Question

I'm trying to port some MATLAB code to Python and am running into some strange behaviour. I am implementing a 5th order Butterworth bandpass filter. The sampling rate is 30 Hz.

Running MATLAB R2012b, Spyder 2.2.0 with Python 2.7, SciPy 0.12.0 on Windows 7 x64.

In MATLAB:

[b,a] = butter(5, [0.75*2/30, 5.0*2/30], 'bandpass');
y = filtfilt(b, a, input_signal)

This is the raw signal:

and the filtered signal:

and its power spectrum:

which makes sense since the normalized bandpass frequencies are 0.05 - 0.33.

I found before that SciPy's butter function does not give the same coefficients as MATLAB so I exported the filter coefficients from MATLAB to Python using hdf5write (see here: https://stackoverflow.com/questions/7117797/export-matlab-variable-to-text-for-python-usage)

In Python:

y = signal.filtfilt(b, a, input_signal, padtype = None)

and the output is:

and its power spectrum:

I used padtype = None because by default it is padtype = 'odd'. However, I've tried all the different padding options and they all look more or less the same.

I'm not entirely sure what's going wrong...any help would be greatly appreciated.

EDIT: Added graph for padtype = "odd" and b and a filter coefficients used

Signal with padtype = "odd":

Seems like there is some underlying signal so I don't think it's the impulse response although that large transient at the beginning is strange.

Python filter coefficients (5th order Butterworth filter):

passband = [0.75*2/30, 5.0*2/30]
b, a = scipy.signal.butter(5, passband, 'bandpass')

b, a are arrays of type float64.

b = array([  5.49209388e-03,   0.00000000e+00,  -2.74604694e-02,
    -1.97776791e-17,   5.49209388e-02,   2.47220989e-17,
    -5.49209388e-02,  -1.97776791e-17,   2.74604694e-02,
     0.00000000e+00,  -5.49209388e-03])
a = array([  1.        ,  -6.52098852,  19.50384534, -35.47189804,
    43.65758795, -38.07760914,  23.83047021, -10.55670367,
     3.16710078,  -0.58122912,   0.04945954])

MATLAB filter coefficients

b =   0.0055, 0, -0.0275, 0, 0.0549, 0, -0.0549, 0, 0.0275, 0, -0.0055
a =   1.0000, -6.5210, 19.5038, -35.4719, 43.6576, -38.0776, 23.8305, 
      -10.5567, 3.1671, -0.5812, 0.0495

Turns out the coefficients are about the same, I believe the issues I had before were due to the frequency being much higher (coefficients must be generated differently in MATLAB and Python).

The input signal is a list of float32, although I've tried converting to array with numpy.array and the result is the same.

EDIT: More information about padtypes and Python vs. MATLAB

Python SciPy's filtfilt function includes a parameter called padtype which indicates the type of padding extended on both sides of the signal. This padding serves to reduce transients. Odd and even are descriptors of the type of symmetry these extensions have with the endpoints (http://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.filtfilt.html).

Some illustrations of the difference in padtypes:

padtype = "even"

padtype = "odd"

padtype = None

Based on the results, it appears that MATLAB uses odd padding (which is also the default for Python):

Answer 1

It works for me:

from scipy.io import loadmat
from scipy.signal import butter, filtfilt
from matplotlib.pyplot import plot

signaldata = loadmat('signaldata.mat')

input_signal = signaldata['input_signal'][0]

passband = [0.75*2/30, 5.0*2/30]
b, a = butter(5, passband, 'bandpass')

y = filtfilt(b, a, input_signal)
plot(y)

I don't know what your issue is but maybe you can figure it out from here.

Answer 2

According to the newsgroup https://mail.python.org/pipermail/scipy-user/2014-April/035648.html,

The difference is in the default padding length. In matlab's filtfilt, it is 3*(max(len(a), len(b)) - 1), and in scipy's filtfilt, it is 3*max(len(a), len(b)).

Therefore, I modified the padlen when invoking filtfilt in Python. And the results between MATLAB and Python will be the same.

MATLAB code

input_signal = [1,2,3,4,5,1,2,3,4,5,1,2,3,4,1,2,3,4,5,1,2,3,4,5,1,2,3,4,5,1,2,3,4,5]
[b,a] = butter(5, [0.75*2/30, 5.0*2/30], 'bandpass');
y = filtfilt(b, a, input_signal)

Python with SciPy code

from scipy.signal import butter, filtfilt
input_signal = [1,2,3,4,5,1,2,3,4,5,1,2,3,4,1,2,3,4,5,1,2,3,4,5,1,2,3,4,5,1,2,3,4,5]
[b,a] = butter(5, [0.75*2/30, 5.0*2/30], 'bandpass');
y = filtfilt(b, a, input_signal, padtype = 'odd', padlen=3*(max(len(b),len(a))-1))

Answer 3

I got the same results in Paco's example, but changing the input signal to a vector of ones does not yield same results in MATLAB and Python (Scipy 1.0.0):

MATLAB:

input_signal = ones(10000,1);
[b,a] = butter(5, [0.75*2/30, 5.0*2/30], 'bandpass');
y = filtfilt(b, a, input_signal);
y(1:5)'

ans =

    1.0e-13 *

    -0.306211857054885  -0.468000282200120  -0.612592864567431  -0.736885525141862  -0.838459113089290

Python:

from numpy import ones
from scipy.signal import butter, filtfilt
input_signal = ones((10000,))
[b,a] = butter(5, [0.75*2/30, 5.0*2/30], 'bandpass');
y = filtfilt(b, a, input_signal, padtype = 'odd', padlen=3*(max(len(b),len(a))-1))
y[:5]

Out[49]: 
array([-2.52233927e-14, -1.79580673e-14, -1.07857124e-14, -3.97094772e-15, 2.30161154e-15])

Answer 4

After much trial and error, I found the following two statements to be equivalent:

% Matlab
rf_data_filt = filtfilt(b, a, rf_data);

# Python
rf_data_filt = filtfilt(b, a, rf_data, axis=0, padtype='odd', padlen=3*(max(len(b),len(a))-1))

It seems that three kwargs are necessary: axis=0, padtype='odd', padlen=3*(max(len(b),len(a))-1)