MATLAB仿真级联型IIR滤波器系数及运算字长的量化性能

MATLAB仿真级联型IIR滤波器系数及运算字长的量化性能

本例对IIR滤波器的差分方程
$$y(n)=\sum _{i=0}^{M}x(n-i)b(i)+\sum _{l=1}^{N}y(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滤波器系数

图1 直接型结构

  级联型IIR滤波器，实际是将一个级数较大的IIR滤波器分解成多个技术小于或等于3的直接型IIR滤波器，其中每个IIR滤波器均可以看成独立的结构，只是将前一级的IIR滤波器的输出作为后一级IIR滤波器的输入而已。其结构如图2所示。

图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)=\frac{{\sum }_{i=0}^M{b}_i{z}^{-i}}{{\sum }_{l=0}^N{a}_l{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}$$
接下来就是要把
$$\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}$$
这一部分写成多个二阶滤波器的乘积形式,比如
$$\frac{(z-n_1）(z-n_2)}{(z-m_1)(z-m_2)}\ast ...\ast \frac{(z-n_{M-1}）(z-n_M)}{(z-m_{N-1})(z-m_{N-2})}$$
这一步是通过roots函数求出有所根，再每隔两个根遍历一次，使用poly函数进行逆运算，重新计算对应二阶滤波器直接型的系数。最终得到
$$\frac{1+d_1{z}^{-1}+d_2{z}^{-2}}{1+c_1{z}^{-1}+c_2{z}^{-2}}\ast ...\ast \frac{1+d_{M-1}{z}^{-1}+d_M{z}^{-2}}{1+c_{N-1}{z}^{-1}+c_{N-2}{z}^{-2}}$$
式中每个二阶滤波器分子分母的系数分别保存在A、B的每一行中，比如A{1,:}保存第一个二阶滤波器的分子系数{1,d₁,d₂}，B{1,:}保存第一个二阶滤波器的分母系数{1,c₁,c₂}。

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)

测试结果为