在LZ给出的SOS matrix和Scale Factors参数中已包含有滤波器的信息。这里给出2种方法形成滤波器参数,这是一个18阶的IIR滤波器,一种是分解成9个2阶IIR滤波器串接,另一种是18阶IIR滤波器:
SOS=[1 -1.753478595924798 1 1 -1.715929423633677 0.989596471431086;...
1 -1.753478595924798 1 1 -1.770907559180061 0.990608353245765;...
1 -1.753478595924798 1 1 -1.702996452119821 0.970500292633477;...
1 -1.753478595924798 1 1 -1.752339882018184 0.973017861982348;...
1 -1.753478595924798 1 1 -1.696640367244520 0.955631016351785;...
1 -1.753478595924798 1 1 -1.733851156227253 0.958455470476445;...
1 -1.753478595924798 1 1 -1.697331544768989 0.946633764727194;...
1 -1.753478595924798 1 1 -1.717331241642967 0.948461822736893;...
1 -1.753478595924798 1 1 -1.704617601062175 0.944269642097509];
Sv=[0.994654917981474
0.994654917981474
0.985575452484598
0.985575452484598
0.978354972045005
0.978354972045005
0.973726850565062
0.973726850565062
0.972134821048755
1.000000000000000];
Tdb=zeros(1,501);
for k=1: 9
B(k,1:3)=SOS(k,1:3)*Sv(k);;
A(k,1:3)=SOS(k,4:6);
[db, mag, pha, w]=freqz_m1(B(k,:),A(k,:));
DB(k,:)=db;
Tdb=Tdb+db;
end
b0=Sv(10);
Tdb=Tdb+20*log10(b0);
figure
plot(w/pi*50,Tdb); grid; axis([0 30 -400 5]);
xlabel('Frequency (kHz)'); ylabel('Magnitude (dB)');
fprintf('B=\n');
for k=1 : 9
fprintf('%4d %5.6f %5.6f %5.6f\n',k,B(k,:));
end
fprintf('\n\n');
fprintf('A=\n');
for k=1 : 9
fprintf('%4d %5.6f %5.6f %5.6f\n',k,A(k,:));
end
fprintf('\n\n');
[b, a]=cas2dir(b0, B, A);
figure
[db, mag, pha, w]=freqz_m1(b,a);
plot(w/pi*50,db); grid; axis([0 30 -400 5]);
xlabel('Frequency (kHz)'); ylabel('Magnitude (dB)');
fprintf('b=\n');
fprintf('%5.6f %5.6f %5.6f %5.6f %5.6f\n',b)
fprintf('\n\n');
fprintf('a=\n');
fprintf('%5.6f %5.6f %5.6f %5.6f %5.6f\n',a)
fprintf('\n\n');
其中用到freqz_m1为
function [db, mag, pha, w]=freqz_m1(b,a);
%Modified version of freqz subroutine
[H,w]=freqz(b,a,1000,'whole');
H=(H(1:501))'; w=(w(1:501))';
mag=abs(H);
db=20*log10(mag+eps);
pha=angle(H);
9个2阶IIR滤波器串接和18阶IIR滤波器的响应图如下,它们的参数为:
B=
1 0.994655 -1.744106 0.994655
2 0.994655 -1.744106 0.994655
3 0.985575 -1.728185 0.985575
4 0.985575 -1.728185 0.985575
5 0.978355 -1.715525 0.978355
6 0.978355 -1.715525 0.978355
7 0.973727 -1.707409 0.973727
8 0.973727 -1.707409 0.973727
9 0.972135 -1.704618 0.972135
A=
1 1.000000 -1.715929 0.989596
2 1.000000 -1.770908 0.990608
3 1.000000 -1.702996 0.970500
4 1.000000 -1.752340 0.973018
5 1.000000 -1.696640 0.955631
6 1.000000 -1.733851 0.958455
7 1.000000 -1.697332 0.946634
8 1.000000 -1.717331 0.948462
9 1.000000 -1.704618 0.944270
b=
0.847849 -13.380160 101.477935 -491.012808 1697.382412
-4449.362152 9161.799827 -15148.109667 20383.171259 -22485.628984
20383.171259 -15148.109667 9161.799827 -4449.362152 1697.382412
-491.012808 101.477935 -13.380160 0.847849
a=
1.000000 -15.491945 115.341293 -547.875022 1859.309795
-4784.778030 9672.660185 -15701.246921 20742.817747 -22466.281798
19995.763845 -14590.650556 8664.764207 -4131.835186 1547.759375
-439.646743 89.223268 -11.552358 0.718847
从图上可以看出9个2阶IIR滤波器串接的响应曲线要比18阶IIR滤波器的响应曲线好,这是因为18阶IIR滤波器对系数的精度和计算的精度都有更高的要求,而在PC机中的双精度已不能满足其要求,而产生的误差造成的。
[本帖最后由 songzy41 于 2009-6-20 06:42 编辑]
qz33a.jpg
(32.72 KB, 下载次数: 5)
2009-6-20 06:40 上传
9个2阶IIR滤波器串接的响应曲线
qz33b.jpg
(33.04 KB, 下载次数: 0)
2009-6-20 06:40 上传
18阶IIR滤波器的响应曲线