[Octave-bug-tracker] [bug #47038] signal package: Ringing, possibly unwanted, with filtfilt

[Dan Sebald]

Follow-up Comment #2, bug #47038 (project octave):

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.


    _______________________________________________________

Reply to this item at:

  <http://savannah.gnu.org/bugs/?47038>

_______________________________________________
  Message sent via/by Savannah
  http://savannah.gnu.org/