0
$\begingroup$

I have designed four IIR bandpass filters in MATLAB. The filters have their poles and zeros inside the ROC in Z-domain.I have checked that the filter is stable using the isstable function in matlab.Also I have plotted the poles and zeros to verify myself that they are within the unit circle.But when I implement it in code using the coefficients, the output is growing after a certain number of input samples. Is the filter stable or not? If not, how to make it stable. I don't know any other method to do this and I need some help.

Thanks.

This is the code to create and check the filetr for stability.


fs=50

%creating the filter 
[A1,B1,C1,D1] = ellip(6,0.1,50,[0.1  2.95 ]/25);
[A2,B2,C2,D2] = ellip(6,0.2,50,[ 3   7 ]/25);
[A3,B3,C3,D3] = ellip(8,0.05,50,[ 7  13 ]/25);
[A4,B4,C4,D4] = ellip(10,0.1,50,[ 13 24 ]/25);

%getting second order coefficients of size Nx6, where N is order of the filter
s1 = ss2sos(A1,B1,C1,D1);
s2 = ss2sos(A2,B2,C2,D2);
s3 = ss2sos(A3,B3,C3,D3);
s4 = ss2sos(A4,B4,C4,D4);

%create cascaded biquad filter from the coefficients
biquad_0_3   = dfilt.df1sos(s1)
biquad_3_7   = dfilt.df1sos(s2)
biquad_7_13  = dfilt.df1sos(s3)
biquad_13_25 = dfilt.df1sos(s4)

%CHECK STABILITY
stable1 = isstable(s1);
stable2 = isstable(s2);
stable3 = isstable(s3);
stable4 = isstable(s4);

zplane(ss2zp(A1,B1,C1,D1))
zplane(ss2zp(A2,B2,C2,D2))
zplane(ss2zp(A3,B3,C3,D3))
zplane(ss2zp(A4,B4,C4,D4))

%plots magnitude response and group delay
  fvtool(biquad_0_3,biquad_3_7,biquad_7_13,biquad_13_25); 
  grpdelay(biquad_0_3)
  grpdelay(biquad_3_7)
  grpdelay(biquad_7_13)
  grpdelay(biquad_13_25) 
$\endgroup$
1
  • 1
    $\begingroup$ Please supply sample plots of the input and output signals as well as the code of the acutal usage of the filters. $\endgroup$
    – Max
    Mar 18, 2022 at 10:52

2 Answers 2

5
$\begingroup$

Mostly likely numerical problems.

These are steep filters with high order and very low cut off frequency. Technically they are stable since the poles are inside the unit circle but they are EXTREMELY close to the unit circle, so it's vulnerable to any type of numerical noise.

So it depends a lot on how the filter operation is implemented. Anything that uses "transfer function" representation ([b ,a]) will not work here. As long as you keep this as cascaded second order section form this is can still work, but the implementation has do be done very carefully.

Here is a simple example: if we take the poles from the first filter, turn them into a denominator polynomial and back into poles again, we get poles that are outside of the unit circle.

% bandpass
[z,p,k] = ellip(6,0.1,50,[0.1  2.95 ]/25);
% to transfer function and back
p1 = roots(poly(p));
% count unstable poles
fprintf('   Original: %d poles outside the unit circle\n',sum(abs(p)>=1));
fprintf('TF and back: %d poles outside the unit circle\n',sum(abs(p1)>=1));
$\endgroup$
1
  • $\begingroup$ The last section of MATLAB's butterworth filter page has a nice plot showing this instability using the transfer function form. $\endgroup$
    – Ash
    Mar 18, 2022 at 20:27
4
$\begingroup$

The [B, A] = ellip(6,0.1,50,[0.1, 2.95]/25); filter has three poles outside the unit circle.

The [B, A] = ellip(6,0.2,50,[3 ,7]/25); filter is stable.

The [B, A] = ellip(8,0.05,50,[7, 13]/25); filter is stable.

The [B, A] = ellip(10,0.1,50,[13, 24]/25); filter is stable.

$\endgroup$
2
  • $\begingroup$ That's the problem with Matlab's implementation. You can do [z,p,k] = ellip(6,0.1,50,[0.1, 2.95]/25); instead and the poles stay where they belong. $\endgroup$
    – Hilmar
    Mar 18, 2022 at 14:18
  • $\begingroup$ Using [B, A] = ellip(6,0.1,50,[0.1, 2.95]/25); and [Z,P,Q] = tf2zpl(B,A); produces a P vector showing three poles outside the unit circle. $\endgroup$ Mar 19, 2022 at 11:03

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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