bugGNU Octave - Bugs: bug #47038, [octave forge] (signal) Ringing,...

 
 

bug #47038: [octave forge] (signal) Ringing, possibly unwanted, with filtfilt

Submitter:  Alexander Klein <matalex>
Submitted:  Mon 01 Feb 2016 12:23:16 PM UTC
   
 
Category:  Octave Package Severity:  3 - Normal
Priority:  1 - Later Item Group:  Inaccurate Result
Status:  Need Info Assigned to:  None
Originator Name:  Open/Closed:  * Open
Release:  * dev Operating System:  * Any
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Wed 10 Feb 2016 09:49:30 AM UTC, comment #2: 

I'd say this definitely needs to be looked at.  On my system I had a filtfilt.m version in my path from probably 5+ years ago and it is extremely simple, i.e., essentially calling filter in both directions after tagging some zeros at the end of the input signal to capture the filter transient and use that as the start of the reverse filter.  Here is the comment from that old version:

## TODO: In Matlab filtfilt `reduces filter startup transients by carefully
## TODO:    choosing initial conditions, and by prepending onto the input
## TODO:    sequence a short, reflected piece of the input sequence'.
## TODO:    Once filtic is written, use that here.

There is a filtic() function, but I don't see it called anywhere in the latest filtfilt.m.  Not sure it applies.

I experimented a bit, and I found that the main contributing factor is this line:


    v = [2*x(1,c)-x((lrefl+1):-1:2,c); x(:,c);
         2*x(end,c)-x((end-1):-1:end-lrefl,c)]; # a column vector


If I change that line to the following:


    v = [zeros(size(2*x(1,c)-x((lrefl+1):-1:2,c))); x(:,c);
         zeros(size(2*x(end,c)-x((end-1):-1:end-lrefl,c)))]; # a column vector


the transient at the front end goes away.  I guess that is no surprise, but why the 2 of 2*x(1,c)?  Consider if filter length is 2

x(1,1) = 3
x(2,1) = 1.5
x(3,1) = -0.5
...

Then

v(1) = 3 + 3 - (-0.5) = 6.5
v(2) = 3 + 3 - (1.5)  = 4.5
v(3) = 3
v(4) = 1.5
v(5) = -0.5

The point is that there is always some additional bias in this padding meant to prime the filter.

I'm not following the logic of using si*v(1) for the initial state either.  I took out the padding at the front and still see a small transient that isn't at the level of the random signal.  Guess I'd have to read the paper referenced within, but I don't see how this initialization scheme makes sense.

Further note, your example runs very slow because filtfilt.m is using a for-loop whereas filter might not be.  It is clear that the F = flipud ( filter ( b, a, flipud ( filter ( b, a, M ) ) ) ); line is faster than filtfilt.

Dan Sebald <sebald>
Mon 01 Feb 2016 07:24:57 PM UTC, comment #1: 

I observe what you describe with your given example. Can someone test this example in Matlab to see what the comparison looks like there?

Mike Miller <mtmiller>
Group Member
Mon 01 Feb 2016 12:23:16 PM UTC, original submission:  

Filtfilt may introduce strange artefacts into the spectra and time series of filtered signals.

This seems to result from the code computing the initial filter state, and manifests itself as a ringing artefact at both ends of the filtered signal, which does not appear when simply applying "filter" in both directions.

See the attached file for sample code plotting a short demo.

While being far from catastrophic for most applications, the effect is certainly strange enough to cause a certain amount of headache to unprepeared users, as it does not seem to be documented anywhere, and it might break assumptions regarding expectation values towards both ends of time series.

Is this really how we want in this function to work?

Alexander Klein <matalex>

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:
   
   
Comment:
   

Attached Files
file #36229:  source_me added by matalex (777B - application/octet-stream)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by jwe (Updated the item)
  • -email is unavailable- added by sebald (Posted a comment)
  • -email is unavailable- added by matalex (Submitted the item)
  •  

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

    Only group members can vote.

     

    Follow 13 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2019-02-27 mtmiller Assigned tomtmiller None
        Releaseother dev
    2019-02-27 mtmiller Carbon-CopyRemoved 80942 -
    2018-04-03 mtmiller Priority5 - Normal 1 - Later
        StatusPostponed Need Info
    2018-03-21 mtmiller StatusNeed Info Postponed
    2017-08-13 jwe Summarysignal package: Ringing, possibly unwanted, with filtfilt [octave forge] (signal) Ringing, possibly unwanted, with filtfilt
    2016-03-29 mtmiller Assigned toNone mtmiller
    2016-02-01 mtmiller StatusNone Need Info
        Release3.8.2 other
        Operating SystemBSD Any
        SummaryRinging, possibly unwanted, with filtfilt signal package: Ringing, possibly unwanted, with filtfilt
    2016-02-01 matalex Attached File- Added source_me, #36229

    Back to the top

    Powered by Savane 3.13-4448.
    Corresponding source code