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:

[Paul]

@Star Strider

Why did you delete the original Answer to this question? I think it contained useful information for anyone who was following this question or might come across it in the future.

Answer 1

[Paul]

The filter is not unstable. Rather, the issue is an interaction between how filtfilt works and the implementation of the biquadratic sections of the filter as designed.

Design the filter

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);

This call no longer works (assuming it did when this question was first posted)

try
Hd = designfilt(h, 'cheby2', ...
    'MatchExactly', 'stopband', ...
    'SystemObject', true);
catch ME
    ME
end

ME = 
  MException with properties:

    identifier: 'signal:designfilt:InputsMustBePVP'
       message: 'All inputs except the filter response must be in the form of name-value pairs.'
         cause: {}
         stack: [5x1 struct]
    Correction: []

But this does, so we'll go with this

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

Define an input for testing

rng(100)
N = 2e5;
x = randn(1,N);

Hd is a dsp.SOSFilter object

whos Hd

  Name      Size            Bytes  Class            Attributes

  Hd        1x1                 8  dsp.SOSFilter              

which does not have filtfilt() method (I wonder why not?)

methods(Hd)

Methods for class dsp.SOSFilter:

SOSFilter           firtype             impz                islinphase          measure             phasedelay          scale               stepz               
cascade             freqrespest         impzlength          ismaxphase          noisepsd            phasez              scalecheck          tf                  
clone               freqrespopts        info                isminphase          noisepsdopts        realizemdl          scaleopts           zerophase           
coeffs              freqz               isLocked            isreal              order               release             specifyall          zpk                 
cost                freqzmr             isallpass           issos               outputDelay         reorder             ss                  zplane              
cumsec              grpdelay            isfir               isstable            parallel            reset               step                

Static methods:

helpFilterAnalysis  helpFixedPoint      

Call "methods('handle')" for methods of dsp.SOSFilter inherited from handle.

which is why this call doesn't work (I think it resolves to filtfilt, which doesn't work with these input arguments)

try
    y = filtfilt(Hd,x);
catch ME
    ME
end

ME = 
  MException with properties:

    identifier: 'MATLAB:narginchk:notEnoughInputs'
       message: 'Not enough input arguments.'
         cause: {}
         stack: [5x1 struct]
    Correction: []

The properties of Hd are

Hd

Hd = 
  dsp.SOSFilter with properties:

            Structure: 'Direct form II'
    CoefficientSource: 'Property'
            Numerator: [8x3 double]
          Denominator: [8x3 double]
       HasScaleValues: true
          ScaleValues: [0.7309 0.7309 0.6503 0.6503 0.4907 0.4907 0.2475 0.2475 1]

  Use get to show all properties

which is why this call won't work

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

Instead, we can call filtfilt like so

y = filtfilt([Hd.Numerator,Hd.Denominator],Hd.ScaleValues,x);
figure
plot(y)

The filtfilt operation is not unstable; it just has a massive transient at the start that takes a long time to damp out.

This issue is also discussed here, which includes a link to this Question and Answer from the @MathWorks Support Team that shows two suggested workarounds. Here is the approach using the second of the two suggestions. The other is based on using reorder.

The approach is to convert the second order sections to zero/pole/gain form and then convert back to sos. Without checking the details, the idea is that this process results in better "balancing" amongst all of the sections that eliminates the transient.

Do the conversion and then filtfilt

[z,p,k] = sos2zp([Hd.Numerator,Hd.Denominator],prod(Hd.ScaleValues));
[Hdsos,g] = zp2sos(z,p,k);
y1 = filtfilt(Hdsos,g,x);

The transient is eliminated

figure
plot(y1)

and the output has the expected response relative to the input, at least until the gain gets extremely low.

[h,w] = freqz(Hdsos,8192); h = h*g;
figure
[h1,w1] = tfestimate(x,y1);
plot(w1,db(abs(h1)),w,db(abs(h).^2)),grid