3
$\begingroup$

I observed that filtfilt suffers from an undesired behavior when I provide IIR bandpass filters having steep transition bands. Specifically, the output signal exhibits excessive transient response; nevertheless, this behaviour does not emerge if I use a 'home-made' version of zero phase filtering based on filter. I guess that it is not a numerical issue involving the filter coefficients or structure, nor the filter function.

% design bandpass filter having transition bandwidth of 200 Hz (Fs = 8000)
bp = designfilt('bandpassiir', 'StopbandFrequency1', 50, 'PassbandFrequency1', 250,...
    'PassbandFrequency2', 3600, 'StopbandFrequency2', 3800, 'StopbandAttenuation1', 30,...
    'PassbandRipple', 0.1, 'StopbandAttenuation2', 30, 'SampleRate', 8000, 'DesignMethod', 'cheby2');
% check stability
assert(isstable(bp),'Unstable filter');
% apply filtfilt to a random (white) long input signal; output signal shows an undesirable transient
x = randn(2^20,1);
y = filtfilt(bp,x);
% apply 'home-made' filtfilt to the same input; output signal shows a more acceptable transient
y2 = flipud(filter(bp,flipud(filter(bp,x))));
% compare effects
figure; semilogy(abs(y-y2));

As a rule of thumb, the effect grows as the transition bands get narrower, while it tends to vanish as they get broader.

Where is the problem? Have I missed some recommendations or hints in the function's help?

$\endgroup$
  • $\begingroup$ I can't use your code with my version of Octave. Could you post the filter coefficients? $\endgroup$ – Matt L. Jan 2 '18 at 14:03
  • $\begingroup$ @MattL. I've tried this with Octave as well. I generated an arbitrary 2nd order Chebyshev filter and the phenomenon occurs too, so I believe the problem the OP has is not related to that specific filter. I used [b,a] = cheby2(10,30,[5/800, 380/800]);. $\endgroup$ – Tendero Jan 2 '18 at 14:19
  • $\begingroup$ @Tendero: OK, good to know. $\endgroup$ – Matt L. Jan 2 '18 at 14:31
  • $\begingroup$ @Tendero: "2nd order Chebyshev filter"? That filter has order 20. And it is numerically unstable, because it has poles extremely close to the unit circle. So no wonder that things can go wrong, especially with extremely long input signals (and they're even longer when that reflection method is used). $\endgroup$ – Matt L. Jan 2 '18 at 17:20
  • $\begingroup$ @MattL. Oops, sorry about that, I meant "2nd type", my bad. And regarding the unstability, you may be right - I just put arbitrary values for every parameter just to see what could happen and forgot to check that out. I'll try with another filter and edit my answer. $\endgroup$ – Tendero Jan 2 '18 at 17:36
5
$\begingroup$

It must have to do with the initial conditions used by the function filtfilt.m. The idea is to match initial conditions in a way such that startup and end transients are minimized. This, however, doesn't always seem to work, and it appears that for your filter specifications it actually does more harm than good. As far as I know there is no way to tell filtfilt.m not to mess with initial conditions. What remains in this case is to simply implement a straightforward version, just as you did.

Also take a look at this related answer.

EDIT:

I tried to reproduce your results with Octave, but I don't get any of the effects that you describe. I use a $10^{th}$ order Chebyshev 2 filter (which is the lowest order satisfying your specs), and I used a modified version of filtfilt allowing me to both switch on and off the use of initial conditions and the extrapolation of the input signal. There are two explanations for this: either your version of Matlab's filtfilt has a bug, or your filter is substantially different from mine and it introduces numerical problems (due to poles very close to the unit circle), which happen to manifest themselves when the input signal is extrapolated, and/or when special initial conditions are used (as is the case with filtfilt). If you provide your filter coefficients we can probably solve this mystery.

EDIT 2:

Tendero kindly reproduced the filter coefficients using the Matlab commands provided in the question:

b = [0.682585947577945 -0.012906667355225 -3.354515884719406 0.025370108447479 6.651816094305592 -0.000000000000000 -6.651816094305592 -0.025370108447479 3.354515884719406 0.012906667355225 -0.682585947577945];
a = [1.000000000000000 -0.194376467633997 -4.144930396202897 0.592472615834755 7.013126532356848 -0.691483043645839 -6.033921428993976 0.364443043103532 2.633924737223088 -0.072944768631229 -0.465923573748158];

I used these coefficients in Octave with the function filtfilt.m with a random input signal x = randn(2^20,1) just like in the question. As mentioned above, I used a modified version of filtfilt.m to be able to switch on and off the different measures for reducing transients:

  1. no extrapolation of the input signal, zero initial conditions
  2. no extrapolation of the input signal, optimal initial conditions (as implemented in the standard version of filtfilt.m)
  3. extrapolation of the input signal, zero initial conditions
  4. extrapolation of the input signal, optimal initial conditions

The last combination is the processing implemented in the standard version of filtfilt.m.

The plot below shows the first and last $50$ samples of the $4$ different output signals. For all other samples the signals are virtually identical. enter image description here

As expected, one can see slight differences during startup and at the end, but these differences are very small, and the problem of the OP cannot be reproduced using the original filter coefficients and the Octave function filtfilt.m. There must be either a bug in the Matlab version of filtfilt.m used by the OP, or another problem that is neither related to the filter coefficients nor to the function filtfilt.m.

$\endgroup$
  • 1
    $\begingroup$ You could try to make by splitting the filter in second order sections (does not always work) or, alternatively, search for filtfilthd on the MathWorks FileExchange which offers additional options / alternatives to the reflection approach that is used by Matlab’s filtfilt. $\endgroup$ – user883521 Jan 2 '18 at 17:52
  • $\begingroup$ Matt, I have the coefficients from the OP's question. I'll post them here so maybe you can try something with them. num = 0.6826 -0.0129 -3.3545 0.0254 6.6518 -0.0000 -6.6518 -0.0254 3.3545 0.0129 -0.6826 den = 1.0000 -0.1944 -4.1449 0.5925 7.0131 -0.6915 -6.0339 0.3644 2.6339 -0.0729 -0.4659 $\endgroup$ – Tendero Jan 7 '18 at 15:50
  • $\begingroup$ @MattL. As you had pointed out, I had used a filter that was numerically unstable, so maybe the OP's filter is as well. I deleted it because I thought it was not correct to say 'Hey, the same happens to me but my filter is okay' when it was not. Maybe you can try with those coefficients in Octave and check what happens (I don't have Octave here, just MATLAB, and filtfilt script is much more complicated to understand and read than Octave's version). $\endgroup$ – Tendero Jan 7 '18 at 16:02
  • $\begingroup$ Yes, I copy-pasted the first three lines of the OP's code. Then I did [num,den] = tf(bp); to get the coefficients. $\endgroup$ – Tendero Jan 7 '18 at 16:39
  • $\begingroup$ I believe that the fact that the coefficients are truncated leads to different results. I tried using full-precision coefficients and the frequency response is pretty different to that obtained using just 4 decimal positions (as I posted above) for low frequencies. I don't know how I could give you the full-precision values, though. Any ideas? $\endgroup$ – Tendero Jan 7 '18 at 16:44
2
$\begingroup$

Looking at the script filtfilt.m, there are two things that differ from the straightforward implementation you did in one single line.

First of all, there is what Matt stated in his answer. The function filtfilt uses initial conditions when it calls the function filter (you can see Mathworks' documentation).

The other difference is that when filtfilt calls filter, it does not use x as input but another array created within the function. I don't understand why or what that signal is, but I believe that's the reason why we are getting that weird peak. In Octave's version, what filtfilt does is the following (putting filter delays initialization to zero, i.e. the modified filtfilt):

lx = size(x,1);
lb = length(b);
la = length(a);
n = max(lb, la);
lrefl = 3 * (n - 1);
if la < n, a(n) = 0; endif
if lb < n, b(n) = 0; endif

v = [ 2*x(1) - x((lrefl+1):-1:2);
      x(:);
      2*x(end) - x((end-1):-1:end-lrefl)];

v = filter(b,a,v);                   # forward filter
v = flipud(filter(b,a,flipud(v)));   # reverse filter
y = v((lrefl+1):(lx+lrefl));

Maybe someone can help us understand what that v represents and why it is useful, but that's the only difference I can spot between filtfilt and your 'homemade' implementation, and so what I expect to be the reason for the undesired transient.

$\endgroup$
  • 3
    $\begingroup$ I found the following here: "In addition to the forward-reverse filtering, it attempts to minimize startup transients by adjusting initial conditions to match the DC component of the signal and by prepending several filter lengths of a flipped, reflected copy of the input signal." $\endgroup$ – Matt L. Jan 2 '18 at 15:30
  • $\begingroup$ @MattL. Nice finding! I'll check the references in that link when I get home and see if there is anything there in order to understand the topic in more depth. $\endgroup$ – Tendero Jan 2 '18 at 15:40
  • $\begingroup$ As stated by Matt L. above, the initial v represents the 'padded' signal (padding applied on both sides). In this case it is 'odd' padding, i.e. the first and last few samples are mirrored, both horizontally and vertically, about the end points (hence the 2*x(1) - x(...)). Alternatively one could use 'even' (first & last few samples mirrored about the end point horizontally), 'constant' (constant extrapolation of the end-point values), or 'zero' padding, and the padding length can also be varied. $\endgroup$ – djvg Oct 11 '18 at 7:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.