MATLAB仿真级联型IIR滤波器系数及运算字长的量化性能
本例对IIR滤波器的差分方程
y ( n ) = ∑ i = 0 M x ( n − i ) b ( i ) + ∑ l = 1 N y ( n − 1 ) a ( l ) y(n) = \sum_{i=0}^{M}x(n-i)b(i)+\sum_{l=1}^{N}y(n-1)a(l) y(n)=i=0∑Mx(n−i)b(i)+l=1∑Ny(n−1)a(l)
的滤波器系数及运算字长的量化性能进行仿真,要求绘出系数分别为8比特、12别特,运算字长分别为12别特、24比特的滤波器幅频响应曲线。采用级联型结构对FPGA直接型IIR滤波器的实现_小小低头哥的博客-CSDN博客(上文实现了MATLAB和FPGA对直接型IIR滤波器的实现,这次实现级联型IIR滤波器)中所述的IIR滤波器进行FPGA实现,并仿真测试测试FPGA实现后的滤波效果。其中IIR滤波器量化位数为12比特、系统时钟频率为2KHZ、数据速率为2KHZ、输入数据位数为12比特。
1.将直接型IIR滤波器系数转换成级联型IIR滤波器系数
级联型IIR滤波器,实际是将一个级数较大的IIR滤波器分解成多个技术小于或等于3的直接型IIR滤波器,其中每个IIR滤波器均可以看成独立的结构,只是将前一级的IIR滤波器的输出作为后一级IIR滤波器的输入而已。其结构如图2所示。
在MATLAB仿真中,也是采用将多阶的滤波器分解为多个二阶的滤波器。具体代码如下:
function [b0,B,A]=My_E5_51_dir2cas(b,a)
%将直接型IIR滤波器系数转换为级联型IIR滤波器系数
%b0:增益系数
%B:包含因子系数bk的K行3列矩阵
%A:包含因子系数ak的K行3列矩阵
%a:直接型IIR滤波器分母系数向量
%b:直接型IIR滤波器分子系数向量
%计算增益系数
b0 = b(1)/a(1);
b = b/b(1);
a = a/a(1);
%将分子、分母系数的长度补齐并进行计算
M = length(b); N = length(a);
if M>N
a = [a zeros(1,M-N)];
N = M ;
elseif N > M
b = [b zeros(1,N-M)];
M = N;
end
%级联型IIR滤波器系数向量的初始化
%如果长度为偶数,但是零、极点个数为奇数,也就是阶数为奇数个
%将一个多阶IIR滤波器分解成多个二阶滤波器,
%一个长度为N(阶数为N-1)的滤波器,最多可以出floor(N/2)也就是K个二阶IIR滤波器
K = floor(N/2);
B = zeros(K,3); %每一个二阶滤波器有两个根,那么对应的多项式就有三个系数
A = zeros(K,3); %所以是三列
if K*2 == N %此时阶数为奇数个,最后分解出的二阶滤波器只有一个零、极点,所以手动给
b = [b 0]; %最后一个二阶滤波器加一个值均为0的零极点
a = [a 0];
end
%根据多项式系数,利用roots函数求出所有的根
%从所有根中两两成对取出,组成二阶IIR滤波器的系数的根
%利用cplxpair函数按实部从小到大的顺序成对排序(感觉可以不必这一步)
broots = cplxpair(roots(b));%求出关于分子的所有根
aroots = cplxpair(roots(a));%求出关于分母的所有根
%取出复共轭对的根,并变换成多项式系数
for i=1:2:2*K %每次遍历一次分子、分母二阶滤波器,进行数据处理
Brow = broots(i:1:i+1,:);%分母的两个根
Brow = real(poly(Brow));%逆运算,重新按照两个根生成多项式系数
B(fix(i+1)/2,:) = Brow; %在对应行(对应二阶滤波器) 保存三个多项式系数
Arow = aroots(i:1:i+1,:);%分子的两个根
Arow = real(poly(Arow));%逆运算,重新按照两个根生成多项式系数
A(fix(i+1)/2,:) = Arow; %在对应行(对应二阶滤波器) 保存三个多项式系数
end
代码中计算增益系数一步,可以理解为:本来直接型IIR滤波器公式为:
H
(
Z
)
=
∑
i
=
0
M
b
i
z
−
i
∑
l
=
0
N
a
l
z
l
H(Z)=\frac{\sum_{i=0}^{M}b_iz^{-i}}{\sum_{l=0}^{N}a_lz^{l}}
H(Z)=∑l=0Nalzl∑i=0Mbiz−i
经过增益系数提取后,就变为
H
(
Z
)
=
b
1
a
1
1
+
∑
i
=
1
M
b
i
b
1
z
−
i
1
+
∑
l
=
1
N
a
l
a
1
z
l
H(Z)=\frac{b_1}{a_1}\frac{1+\sum_{i=1}^{M}{\frac{b_i}{b_1}}z^{-i}}{1+\sum_{l=1}^{N}\frac{a_l}{a_1}z^{l}}
H(Z)=a1b11+∑l=1Na1alzl1+∑i=1Mb1biz−i
接下来就是要把
1
+
∑
i
=
1
M
b
i
b
1
z
−
i
1
+
∑
l
=
1
N
a
l
a
1
z
l
\frac{1+\sum_{i=1}^{M}{\frac{b_i}{b_1}}z^{-i}}{1+\sum_{l=1}^{N}\frac{a_l}{a_1}z^{l}}
1+∑l=1Na1alzl1+∑i=1Mb1biz−i
这一部分写成多个二阶滤波器的乘积形式,比如
(
z
−
n
1
)
(
z
−
n
2
)
(
z
−
m
1
)
(
z
−
m
2
)
∗
.
.
.
∗
(
z
−
n
M
−
1
)
(
z
−
n
M
)
(
z
−
m
N
−
1
)
(
z
−
m
N
−
2
)
\frac{(z-n_1)(z-n_2)}{(z-m_1)(z-m_2)}*...*\frac{(z-n_{M-1})(z-n_{M})}{(z-m_{N-1})(z-m_{N-2})}
(z−m1)(z−m2)(z−n1)(z−n2)∗...∗(z−mN−1)(z−mN−2)(z−nM−1)(z−nM)
这一步是通过roots函数求出有所根,再每隔两个根遍历一次,使用poly函数进行逆运算,重新计算对应二阶滤波器直接型的系数。最终得到
1
+
d
1
z
−
1
+
d
2
z
−
2
1
+
c
1
z
−
1
+
c
2
z
−
2
∗
.
.
.
∗
1
+
d
M
−
1
z
−
1
+
d
M
z
−
2
1
+
c
N
−
1
z
−
1
+
c
N
−
2
z
−
2
\frac{1+d_1z^{-1}+d_2z^{-2}}{1+c_1z^{-1}+c_2z^{-2}}*...*\frac{1+d_{M-1}z^{-1}+d_{M}z^{-2}}{1+c_{N-1}z^{-1}+c_{N-2}z^{-2}}
1+c1z−1+c2z−21+d1z−1+d2z−2∗...∗1+cN−1z−1+cN−2z−21+dM−1z−1+dMz−2
式中每个二阶滤波器分子分母的系数分别保存在A、B的每一行中,比如A{1,:}保存第一个二阶滤波器的分子系数{1,d1,d2},B{1,:}保存第一个二阶滤波器的分母系数{1,c1,c2}。
1.1 测试结果
输入下面代码
[b,a] = cheby2(7,60,0.5);
[b0,B,A] = My_E5_51_dir2cas(b,a);
运行后可获得级联型IIR滤波器系数,即
2.对级联型IIR滤波器系数进行量化
由上一步可知,本例的8级IIR滤波器可等效为3个3级IIR滤波器和一个2级IIR滤波器级联组成的滤波器。由于滤波器的增益为0.0145,理论上可将此增益分配给任意一个IIR滤波器,本例分配给第一个级数为3的IIR滤波器,这样有利于降低运算过程中的数据字长。与直接型IIR滤波器系数量化处理过程相似,在进行FPGA实现前必须对级联型IIR滤波器系数进行量化处理。
function [Qb,Qa]=E5_52_Qcoe(b,a,Qcoe);
%b:系统的分子系数
%a:系统的分母系数
%Qnum:量化位数
m=max(max(abs(a),abs(b))); %获取滤波器系数向量中绝对值最大的数
Qm=floor(log2(m/a(1))); %取系数中最大值与a(1)的整数倍
if Qm<log2(m/a(1))
Qm=Qm+1;
end
Qm=2^Qm; %获取量化基准值
Qb=round(b/Qm*(2^(Qcoe-1))) %四舍五入截尾
Qa=round(a/Qm*(2^(Qcoe-1))) %四舍五入截尾
IIR滤波器系数及运算字长的仿真_小小低头哥的博客-CSDN博客本文对此代码有详细介绍
2.1 测试结果
测试代码为
[B(1,:),A(1,:)]=E5_52_Qcoe(B(1,:),A(1,:),12)
[B(2,:),A(2,:)]=E5_52_Qcoe(B(2,:),A(2,:),12)
[B(3,:),A(3,:)]=E5_52_Qcoe(B(3,:),A(3,:),12)
[B(4,:),A(4,:)]=E5_52_Qcoe(b0*B(4,:),A(4,:),12)
测试结果为