[SciPy-User] scipy.signal vs Matlab: filtfilt and reflection
John Krasting - NOAA Federal
john.krasting at noaa.gov
Tue Apr 15 19:34:11 EDT 2014
Hi Warren -
Thanks for your reply. Yes, when I set padlen=9 in my case I am able to
reproduce Matlab's results. Thanks for helping me see the difference!
-John K.
--------------------------------
Dr. John Krasting
Physical Scientist (NOAA Federal)
NOAA/Geophysical Fluid Dynamics Laboratory
Biogeochemistry, Ecosystems, and Climate Group
Princeton University Forrestal Campus
201 Forrestal Road
Princeton, NJ 08540
P. (609) 452-5359
F. (609) 987-5063
On Tue, Apr 15, 2014 at 7:25 PM, Warren Weckesser <
warren.weckesser at gmail.com> wrote:
>
>
> 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
>>
>>
>
> _______________________________________________
> 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/d88e381c/attachment.html>
More information about the SciPy-User
mailing list