[SciPy-User] scipy.signal vs Matlab: filtfilt and reflection

Tue Apr 15 18:33:41 EDT 2014

Hi Warren -

Thanks for your reply.  Sorry for being mistaken, I now see that filtfilt
does the odd padding by default.

I am still seeing differences between Matlab and scipy at the second
decimal place.  This concerns me.  Has anybody else seen this?

I took a sample 15 numbers (pulled from a random dist.), built a
Butterworth filter, and ran filtfilt on the data using Matlab and scipy.
 The results are in the same ballpark, but I am still seeing differences
between Matlab and scipy at the second decimal place.  This concerns me
since I expected to see differences at the lower order bits.  Has anybody
else seen this or know why this is the case?

I've pasted my 15 random numbers and output from both Matlab and scipy
below.

Thanks,
John K.

PS - I tried taking the Butterworth filter from scip

---- scipy v0.13.3, numpy v.1.6.2 ----

data = numpy.array([0.401808033751942, 0.075966691690842,
0.239916153553658, 0.123318934835166, 0.183907788282417, 0.239952525664903,
0.417267069084370, 0.049654430325742, 0.902716109915281, 0.944787189721646,
0.490864092468080, 0.489252638400019, 0.337719409821377, 0.900053846417662,
0.369246781120215])
[
b,a] = scipy.signal.butter(3,0.2)

>>> b
array([ 0.01809893,  0.0542968 ,  0.0542968 ,  0.01809893])

>>> b[0]
0.018098933007514431
>>> b[1]
0.054296799022543286
>>> b[2]
0.054296799022543286
>>> b[3]
0.018098933007514431

>>> a
array([ 1.        , -1.76004188,  1.18289326, -0.27805992])

>>> a[0]
1.0
>>> a[1]
-1.7600418803431686
>>> a[2]
1.1828932620378303
>>> a[3]
-0.27805991763454624

scipy.signal.filtfilt(b,a,data)

array([ 0.40421981,  0.29292145,  0.20773242,  0.1682317 ,  0.18250909,
        0.24674659,  0.34676359,  0.46045534,  0.56190086,  0.62844665,
        0.64819638,  0.62192427,  0.55831894,  0.46881537,  0.36653733])

---- MATLAB v.7.9.0.529 ----

data = [0.401808033751942, 0.075966691690842, 0.239916153553658,
0.123318934835166, 0.183907788282417, 0.239952525664903, 0.417267069084370,
0.049654430325742, 0.902716109915281, 0.944787189721646, 0.490864092468080,
0.489252638400019, 0.337719409821377, 0.900053846417662, 0.369246781120215]

>> [b,a] = butter(3,0.2)

b =

   0.018098933007514   0.054296799022543   0.054296799022543
0.018098933007514

a =

   1.000000000000000  -1.760041880343169   1.182893262037831
 -0.278059917634546

>> filtfilt(b,a,data)

ans =

  Columns 1 through 5

   0.397839293822166   0.286273839374148   0.202969550095165
0.166087236796221   0.182700682986968

  Columns 6 through 10

   0.248436130460523   0.348904396853114   0.461998356964929
0.561930803937072   0.626352300416756

  Columns 11 through 15

   0.643941843065195   0.616392843959916   0.553574502943643
0.468014071436841   0.373188133985282

On Tue, Apr 15, 2014 at 4:38 PM, Warren Weckesser <
warren.weckesser at gmail.com> wrote:

> On 4/15/14, John Krasting - NOAA Federal <john.krasting at noaa.gov> wrote:
> > Hi Scipy Users -
> >
> > Am I correct in reading that filtfilt in scipy.signal (v. 13.0) does not
> > extrapolate data at the beginning and the end of a time series when using
> > the filtfilt function?
> >
> > Based on the Matlab code (
> >
> http://chmielowski.eu/POLITECHNIKA/Dydaktyka/AUTOMATYKA/AutoLab/Matlab/TOOLBOX/SIGNAL/FILTFILT.M
> ),
> > I see that there is a block of code that does the reflection:
> >
> > % Extrapolate beginning and end of data sequence using a "reflection
> > % method".  Slopes of original and extrapolated sequences match at
> > % the end points.
> > % This reduces end effects.
> >     y = [2*x(1)-x((nfact+1):-1:2);x;2*x(len)-x((len-1):-1:len-nfact)];
> >
> > I did not see this in the source for scipy's filtfilt function.   Are
> there
> > any plans to add this as a possible option to the function?
> >
> > Thanks,
> > John K.
> >
>
>
> That matlab code makes an odd extension of the data.  E.g.
>
>
> octave:35> len = 8
> len =  8
> octave:36> x = [1:len]';
> octave:37> x'
> ans =
>
>    1   2   3   4   5   6   7   8
>
> octave:38> nfact = 3;
> octave:39> y =
> [2*x(1)-x((nfact+1):-1:2);x;2*x(len)-x((len-1):-1:len-nfact)];
> octave:40> y'
> ans =
>
>    -2   -1    0    1    2    3    4    5    6    7    8    9   10   11
>
>
> That type of extension is the default behavior of filtfilt
> (
> http://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.filtfilt.html
> ).
>  You can change the type of the extension and the length of the
> padding with the padtype and padlen arguments.
>
> Warren
> _______________________________________________
> SciPy-User mailing list
> SciPy-User at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20140415/c0d5fb0e/attachment.html>