How to do a Biquad (SOS) filtering

Question

Jay 2024-1-24

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2073501-how-to-do-a-biquad-sos-filtering

编辑： Jay 2024-1-24

I have the following filters which I want to implement with cascade of multiple second-order sections (SOS).

%% Filter_1
h = fdesign.bandpass('N,F3dB1,F3dB2', 10, 100e3, 120e3, 400e3);
Hd = design(h, 'butter');
Hd = convert(Hd, 'df2sos');
%% Filter_2
[be,ae] = ellip(4,1,40,0.2);
f = 0:0.001:0.2;
g = grpdelay(be,ae,f,2);
a = max(g)-g;
[num,den]=iirgrpdelay(8, f, [0 0.2], a);
[sos,g] = tf2sos(num,den);

So there are 9 SOS filters (e.g., 5 from the Filter_1 and 4 from the Filter_2) and I want to cascade all these. Is there a way to filter these nine filters all together? Or am I supposed to filter Filter 1 first and then apply Filter 2 later? Any example for each case?

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

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

Answer 1

Paul 2024-1-24

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2073501-how-to-do-a-biquad-sos-filtering#answer_1395896

编辑：Paul 2024-1-24

在 MATLAB Online 中打开

Hi Jay,

I think it works like this if you want to cascade the filters together into a single, sos filter.

%% Filter_1
h = fdesign.bandpass('N,F3dB1,F3dB2', 10, 100e3, 120e3, 400e3);
Hd = design(h, 'butter');
Hd = convert(Hd, 'df2sos');
sos1 = Hd.sosMatrix;
g1   = prod(Hd.ScaleValues)
g1 = 5.9796e-05
%% Filter_2
[be,ae] = ellip(4,1,40,0.2);
f = 0:0.001:0.2;
g = grpdelay(be,ae,f,2);
a = max(g)-g;
[num,den]=iirgrpdelay(8, f, [0 0.2], a);
[sos2,g2] = tf2sos(num,den);
%% Filter_3 is the combined filter
sos3 = [sos1;sos2];
g3   = g1*g2
g3 = 8.7612e-06

Frequency response Hd compared to sos1

[hd,wd] = freqz(Hd);

h1 = freqz(sos1,wd); h1 = h1*g1;

figure

subplot(211);plot(wd,db(abs([hd h1])));grid

subplot(212);plot(wd,180/pi*angle([hd h1]));grid

Frequency response of sos2 and sos3

h2 = freqz(sos2,wd); h2 = h2*g2;
h3 = freqz(sos3,wd); h3 = h3*g3;

Compare response of product to the product of the responses

figure

subplot(211),plot(wd,db(abs([hd.*h2 h3]))),grid

subplot(212);plot(wd,180/pi*angle([hd.*h2 h3])),grid

It appears that sos3,g3 captures the combined response of the the cascade of the original Hd and sos2.

However, a standard sos filter has the poles and zeros paired in a particular way and the sections in a particular order. It's likely that sos3 does not meet these rules. If this is important, one option would be to convert sos3 to another form and then convert back to sos

[z3,p3,k3] = sos2zp(sos3,g3);
[sos4,g4] = zp2sos(z3,p3,k3);
[g3 g4]
ans = 1×2
1.0e-05 *

    0.8761    0.8761
[sos3 sos4]
ans = 9×12
    1.0000         0   -1.0000    1.0000    0.5773    0.9107    1.0000         0   -1.0000    1.0000   -1.1868    0.3521
    1.0000         0   -1.0000    1.0000    0.0144    0.9066    1.0000   -2.0000    1.0000    1.0000   -1.6930    0.7251
    1.0000         0   -1.0000    1.0000    0.4432    0.7776    1.0000    2.0000    1.0000    1.0000    0.2735    0.7265
    1.0000         0   -1.0000    1.0000    0.1146    0.7716    1.0000   -2.0000    1.0000    1.0000   -1.6178    0.7310
    1.0000         0   -1.0000    1.0000    0.2735    0.7265    1.0000    2.0000    1.0000    1.0000    0.1146    0.7716
    1.0000   -3.3703    2.8398    1.0000   -1.1868    0.3521    1.0000   -3.3703    2.8398    1.0000    0.4432    0.7776
    1.0000   -2.3348    1.3791    1.0000   -1.6930    0.7251    1.0000   -2.3348    1.3791    1.0000   -1.5449    0.7850
    1.0000   -2.2132    1.3681    1.0000   -1.6178    0.7310    1.0000   -2.2132    1.3681    1.0000    0.0144    0.9066
    1.0000   -1.9680    1.2739    1.0000   -1.5449    0.7850    1.0000   -1.9680    1.2739    1.0000    0.5773    0.9107

The pole/zero pairing came out different as did the section ordering. But sos4 is, mathematically, the same filter as sos3.

h4 = freqz(sos4,wd); h4 = h4*g4;

figure

subplot(211),plot(wd,db(abs([h3 h4]))),grid

subplot(212);plot(wd,180/pi*angle([h3 h4])),grid

Forward filtering a signal through Hd, and then filtering the output of Hd through sos2 might be tricky. You'd have to figure out how to get the final states out of the Hd object, and the figure out how to transform those states as initial conditions in the filtering with sos2. Keep in mind that sosfilt, as might be used for the sos2 filtering, doesn't even accept initial conditions as an input argument.