Serious Problem with FiltFilt Function

56 次查看（过去 30 天）
Berk Kucukoglu
Berk Kucukoglu 2021-5-2
编辑： Paul 2025-1-12
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:
  1 个评论
Paul
Paul 2025-1-11
编辑：Paul 2025-1-12
@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.

请先登录，再进行评论。

请先登录，再回答此问题。

回答（1 个）

Paul
Paul 2025-1-12
编辑：Paul 2025-1-12
Ran in:
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

类别

在 Help Center 和 File Exchange 中查找有关 Digital Filtering 的更多信息

标签

产品


版本

R2020b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by