Serious Problem with FiltFilt Function

Hello everyone,

I designed a Bandpass Filter using the filterDesigner app. It is an 18th order IIR filter.

Then I export the filter both as df2sos object and "SOS and G" matrices.

Afterwards, I try filtering my data with:

y = filtfilt(Hd,x);

It gives the error: Not enough input arguments.

Well. Why don't I use:

y = filtfilt(Hd.sosMatrix,Hd.scaleValues,x);
% or
y = filtfilt(SOS,G,x);

These do not work for my filter; they result in unstable filters. Even though my originally designed filter is stable, they give unstable results. Here is the MATLAB page for this issue.

So, how to use FiltFilt function with my df2sos object or my "SOS and G" matrices, without getting an unstable result? I'd really appreaciate some help.

Regards,

Berk

Additional Details:

Here is the code to generate my filter object:

Fstop1 = 0.0001;  % First Stopband Frequency
Fpass1 = 0.0004;  % First Passband Frequency
Fpass2 = 0.08;    % Second Passband Frequency
Fstop2 = 0.12;    % Second Stopband Frequency
Astop1 = 50;      % First Stopband Attenuation (dB)
Apass  = 0.5;     % Passband Ripple (dB)
Astop2 = 50;      % Second Stopband Attenuation (dB)

h = fdesign.bandpass('fst1,fp1,fp2,fst2,ast1,ap,ast2', Fstop1, Fpass1, ...
    Fpass2, Fstop2, Astop1, Apass, Astop2);

Hd = designfilt(h, 'cheby2', ...
    'MatchExactly', 'stopband', ...
    'SystemObject', true);

and the screenshot from the filterDesigner app:

Answer 1

[Star Strider]

I have no idea what function was used to design the filter, since ‘FilterDesignerToolbox’ is nowhere to be found in the MATLAB or Signal Processing Toolbox documentation.

Use designfilt to create the filter. See the filtfilt documentation section on d for help with that.

[Berk Kucukoglu]

Hello Star Strider. I meant the filterDesigner app, not the toolbox. Thank you, edited now.

Can you once again check the question, considering the information I have added?

[Star Strider]

I have no idea what you want to do with respect to the ‘fdesign’ call. It is not possible to use them together, so just use one or the other.

Fstop1 = 0.0001;  % First Stopband Frequency
Fpass1 = 0.0004;  % First Passband Frequency
Fpass2 = 0.08;    % Second Passband Frequency
Fstop2 = 0.12;    % Second Stopband Frequency
Astop1 = 50;      % First Stopband Attenuation (dB)
Apass  = 0.5;     % Passband Ripple (dB)
Astop2 = 50;      % Second Stopband Attenuation (dB)

Hd = designfilt('lowpassiir', 'FilterOrder',20, 'StopbandFrequency',Fstop1,... 
    'StopbandAttenuation',60, 'SampleRate',1000);

figure
freqz(Hd)

This functionality is not available on remote platforms.

I cannot run this online (it will not launch freqz for this filter), however it works from a script, producing a Chebyshev II lowpass filter.

[Berk Kucukoglu]

Okay, I did it your way. It still did not work.

I think the problem is not about the way I design my filter. It is about the:

filtfilt(); % Zero-phase filtering

function. Here is how I designed my filter this time:

Fstop1 = 0.0004;  % First Stopband Frequency
Fpass1 = 0.0010;  % First Passband Frequency
Fpass2 = 0.0800;    % Second Passband Frequency
Fstop2 = 0.1200;    % Second Stopband Frequency
Astop1 = 40;      % First Stopband Attenuation (dB)
Apass  = 0.1;     % Passband Ripple (dB)
Astop2 = 40;      % Second Stopband Attenuation (dB)

bpFilt = designfilt('bandpassiir',...
         'PassbandFrequency1',Fpass1,'PassbandFrequency2',Fpass2,...
         'StopbandFrequency1',Fstop1,'StopbandFrequency2',Fstop2,...
         'PassbandRipple',Apass,'StopbandAttenuation1',Astop1,...
         'StopbandAttenuation2',Astop2,'DesignMethod','cheby2',...
         'MatchExactly','passband');

Now here comes the interesting part. I can use two different functions to filter:

y = filter(bpFilt,x); % normal filtering
y = filtfilt(bpFilt,x); % zero-phase filtering

and here are the results:

There seems to be problem with filtfilt function (this happens when the filter order is high).

Hence the topic "Serious Problem with Filtfilt Function".

Note: Those who are curious about the input signal "x" can generate it themselves with the code

Fs = 1e6;
Ts = 1./Fs;
n = 0:1:1*Fs;
f_max = Fs/10;
f = linspace(0,f_max,length(n));
x = cos(2*pi.*f.*n.*Ts); % Basic chirp

[Paul]

If you zoom in on the filtfilt() plot you'll see it looks similar to the filter() plot. The difference appears to be in how filtfilt() computes the initial conditions of the filter states, which is why there are huge transients at the start and end of the sequence. Why filtfilt() computes very, very large initial conditions is still a mystery, at least to me.

[Star Strider]

@Paul — Thank you!

Padding the ends of the signal with a ones vector (usually ones(1,10) or more for a row vector) multiplied by the first and last values (respectively) of the signal reduces the transients filtfilt generates. (Those elements can later be deleted from the filtered signal.) At least, this has been my experience.

[Berk Kucukoglu]

Thank you for the answer @Paul

I think there is something wrong with the Filtfilt() function. Looking at the impulse response of the designed filter, such an enormous transient does not seem reasonable. And if it is reasonable, I really would enjoy a brief explanation on this.

@Star Strider thank you for the tip. I might try this. However, I am still not comfortable using a function that I think is not working properly. And if it is working properly, then I am not comfortable using your tip, haha!

Thank you people for the answers. I hope some MATLAB official will look into the issue; by fixing this or by teaching us some maths, depending on the answer.

[Paul]

@Berk Kucukoglu

I don't believe that the transient is related to the impulse response per se. Rather it is due to this feature of filtfilt:

"filtfilt minimizes start-up and ending transients by matching initial conditions"

filtfilt() uses non-zero initial condtions in the filter, and it is these initial condtions that cause that transient.

I'm sufficiently curious about this that I'm trying to understand exactly what filtfilt() is doing in how it sets the intial conditions in the filter, which are an input argument to filter() inside filtfilt().

[Paul]

I'm becoming a bit more skeptical about the implementation of filtfilt(). Consider the following, where the filter is a harmless, 2nd order filter.

Fs = 2000;
Ts = 1./Fs;
n = 0:1:1*Fs;
f_max = Fs/10;
f = linspace(0,f_max,length(n));
x = cos(2*pi.*f.*n.*Ts); % Basic chirp
 
d = designfilt('lowpassiir', ...        % Response type
       'PassbandFrequency',50, ...     % Frequency constraints
       'StopbandFrequency',100, ...
       'PassbandRipple',4, ...          % Magnitude constraints
       'StopbandAttenuation',15, ...    % 15 for second order filter
       'DesignMethod','ellip', ...      % Design method
       'MatchExactly','passband', ...   % Design method options
       'SampleRate',Fs);               % Sample rate
   
[b,a] = tf(d);
[b;a]

ans = 2×3

0.0577 -0.1050 0.0577 1.0000 -1.9067 0.9231

y1 = filtfilt(d,x);
y2 = filtfilt(b,a,x);

plot([y1 ; y2].')

As we see, in this case filtfilt() returns the same result operating on either the digitalFilter object or the numerator and denominator coefficients directly.

Now change a filter design parameter to increase the order to a fourth-order filter

d = designfilt('lowpassiir', ...        % Response type
       'PassbandFrequency',50, ...     % Frequency constraints
       'StopbandFrequency',100, ...
       'PassbandRipple',4, ...          % Magnitude constraints
       'StopbandAttenuation',55, ...    % 15, 55  for second/fourth order filter
       'DesignMethod','ellip', ...      % Design method
       'MatchExactly','passband', ...   % Design method options
       'SampleRate',Fs);               % Sample rate
   
[b,a] = tf(d);
[b;a]

ans = 2×5

0.0012 -0.0039 0.0055 -0.0039 0.0012 1.0000 -3.8985 5.7250 -3.7529 0.9265

y1 = filtfilt(d,x);
y2 = filtfilt(b,a,x);

plot([y1 ; y2].') 

Now there is a clear difference in the transient at the end, and zooming in at the beginning reveals a difference as well.

plot([y1 ; y2].')
set(gca,'XLim',[0 10])

Same filter, different results, and the results seem to be too different to ascribe to round-off effects for what appears to be very simple filter.

[Berk Kucukoglu]

Hello @Paul,

I am very happy to see that this topic got your attention. I am also curious and trying to understand what is happening with filtfilt() function.

I have read your post and I think you are correct. The filter begins to act in a weird way as the filter order goes high.

Also, I would like to share what I have noticed. It might be a small detail but I find it interesting.

Say I have two different chirp signals:

x1 = cos(2*pi.*f.*n.*Ts); % Basic chirp 1
x2 = sin(2*pi.*f.*n.*Ts); % Basic chirp 2

Normally I would not even call these chirps "different". They are identical except a tiny phase shift. Now I filter these signals using the filtfilt() function.

y_filtfilt_1 = filtfilt(bpFilt,x1); % zero-phase filtering, Chirp 1
y_filtfilt_2 = filtfilt(bpFilt,x2); % zero-phase filtering, Chirp 2

And here comes the results...

Look at the difference in scales. I overlapped both of them and zoomed onto the initial transient:

I cannot understand what's wrong. I will keep searching. Probably won't use filtfilt() for a while now. Doing the search mainly for curiosity as well :)