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

Warren Weckesser warren.weckesser at gmail.com
Tue Apr 15 19:25:27 EDT 2014


On Tue, Apr 15, 2014 at 6:33 PM, John Krasting - NOAA Federal <
john.krasting at noaa.gov> wrote:

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

John,

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)).  If you set padlen to the same value used in
matlab, the results match:

 In [47]: filtfilt(b, a, data, padlen=9)
Out[47]:
array([ 0.39783929,  0.28627384,  0.20296955,  0.16608724,  0.18270068,
        0.24843613,  0.3489044 ,  0.46199836,  0.5619308 ,  0.6263523 ,
        0.64394184,  0.61639284,  0.5535745 ,  0.46801407,  0.37318813])


Warren



> 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
>>
>
>
> _______________________________________________
> 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/8ba2bca0/attachment.html>


More information about the SciPy-User mailing list