IIR滤波器系数及运算字长的仿真

IIR滤波器系数及运算字长的仿真

对IIR滤波器系统差分方程式所表示的IIR滤波器系数及运算字长进行量化性能仿真，要求绘出系数分别为8比特、12比特，运算字长分别为12比特、24比特的滤波器幅频响应曲线。

  IIR滤波器在运算过程中无法做到全精度运算，因为IIR滤波器是一个存在反馈环节的闭环系统，且中间过程存在除法运算。如果要实现全精度运算，则在运算过程中寄存器所需的字长将十分长，因此在进行FPGA实现之前，必须通过仿真确定滤波器系数及运算字长，并且不同结构所对应的IIR滤波器运算字长需要分别通过仿真来确定。

1.对直接型IIR滤波器系数进行量化

  对于一个例子，如采用cheby2函数设计一个阶数为7（级数为8）的低通IIR滤波器，抽样频率为2000HZ、截止频率为500HZ、阻带衰减为60dB，可在MATLAB命令窗中直接输入下面的命令：

[b,a] = cheby2(7,60,0.5); %得到系数
%12比特量化
m = max(max(abs(a),abs(b)));
Qb = round(b/m*(2^(12-1)-1));
Qa = round(a/m*(2^(12-1)-1));

  如果把系数带入到IIR滤波器系统差分方程式中，可得到
$$Qa(1)y(n)=Qb(1)x(n)+Qb(8)x(n-7)+Qb(2)x(n-1)+Qb(7)x(n-6)+Qb(3)x(n-2)+Qb(6)x(n-5)+Qb(4)x[n-3]+Qb(5)x(n-4)-[Qa(2)y(n-1)+Qa(3)y(n-2)+Qa(4)y(n-3)+Qa(5)y(n-4)+Qa(6)y(n-5)+Qa(7)y(n-6)+Qa(8)y(n-7)]$$

$$900y(n)=13[x(n)+x(n-7)]+38[x(n-1)+x(n-6)]+74[x(n-2)+x(n-5)]+99[x[n-3]+x(n-4)]-[-1623y(n-1)+2047y(n-2)-1427y(n-3)+725y(n-4)-215y(n-5)+42y(n-6)-3y(n-7)]$$
不难发现，左边乘了一个常系数900，每次计算完右边的值后，最后还要除以900，求出最后结果。但是除法运算是十分耗费资源的，如果能除以2的指数倍，就可通过移位的方法近似实现除法运算。所以在量化时一般有意把量化后的IIR滤波器分母系数的第一项设置为2的整数冥次方的形式。其命令为

[b,a] = cheby2(7,60,0.5);
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^(12-1)-1));%四舍五入截尾
Qa = round(a/Qm*(2^(12-1)-1));%四舍五入截尾

最后得到Qa（1）等于512，让我们来分析一下为什么这个操作能得到这个结果。

Qm = floor(log2(m/a(1)));    %获取系数中最大值与a(1)的整数倍
if Qm < log2(m/a(1))
    Qm = Qm + 1;
end

这一段代码就相当于是

Qm = ceil(log2(m/a(1))); 

用MATLAB求出IIR滤波器系数的a(1)始终为1，并假设
$$2^{L1}\le m<2^{L1+1}其中L1为整数$$
则
$$L1\le Qm\le L1+1$$
而Qm=2^Qm就相当于Qm=2^(L1+1)||L1，得到的Qm$\ge$ m且Qm<2*m,可以确定的是其为一个整数，范围在m-2m之间，假设为L1+1(大部分情况也都这样)。则
$$round(a/Qm\ast (2^{(}12-1)-1))=round(a/2^{L1+1}\ast (2^{11}-1))=a\ast 2^{11-(L1+1)}$$
由于a(1)=1，所以得到的Qa(1)=2^(11-(L1+1))。也就是得到了想要的常数。如例子中得到m=2.2735,L1+1=2,理论值Qa(1)=2^(11-2)=512。与仿真值相同。

2.MATLAB仿真程序

fs = 2000;      %抽样频率
fc = 500;       %阻带截止频率
Rs = 60;        %阻带衰减（dB）
N = 7;          %IIR滤波器阶数
Qcoe = [8 12];  %IIR滤波器系数字长
Qout = [12 24]; %IIR滤波器输出字长
delta = [1,zeros(1,511)];   %单位冲激响应作为输入信号
[b,a] = cheby2(N,Rs,fc/(fs/2)); %设计切比雪夫Ⅱ型低通滤波器
%对IIR滤波器系数进行量化，四舍五入截尾
m = max(max(abs(a),abs(b)));
Qm = floor(log2(m/a(1)));
if Qm < log2(m/a(1))
    Qm = Qm + 1;
end
Qm = 2^Qm;
Qb8 = round(b/Qm*(2^(Qcoe(1)-1)-1));
Qa8 = round(a/Qm*(2^(Qcoe(1)-1)-1));
Qb12 = round(b/Qm*(2^(Qcoe(2)-1)-1));
Qa12 = round(a/Qm*(2^(Qcoe(2)-1)-1));
%求理想幅度响应
y=filter(b,a,delta);
%求量化后的幅度响应，Qy=My_E5_32_QuantIIRDirectArith(Qb,Qa,din,Qcoe,Qout)
%为自编的根据系数及输出数据位数的量化函数
%IIR滤波器输出的函数
c8o12=My_E5_32_QuantIIRDirectArith(Qb8,Qa8,delta,Qcoe(1),Qout(1));
c8o24=My_E5_32_QuantIIRDirectArith(Qb8,Qa8,delta,Qcoe(1),Qout(2));
c12o12=My_E5_32_QuantIIRDirectArith(Qb12,Qa12,delta,Qcoe(2),Qout(1));
c12o24=My_E5_32_QuantIIRDirectArith(Qb12,Qa12,delta,Qcoe(2),Qout(2));
%求IIR滤波器输出幅频响应
Fy=20*log10(abs(fft(y)));   Fy=Fy-max(Fy);
Fc8o12=20*log10(abs(fft(c8o12)));   Fc8o12=Fc8o12-max(Fc8o12);
Fc8o24=20*log10(abs(fft(c8o24)));   Fc8o24=Fc8o24-max(Fc8o24);
Fc12o12=20*log10(abs(fft(c12o12))); Fc12o12=Fc12o12-max(Fc12o12);
Fc12o24=20*log10(abs(fft(c12o24))); Fc12o24=Fc12o24-max(Fc12o24);
%设置幅频响应的横坐标，单位为HZ
x_f=[0:(fs/length(Fy)):fs-fs/length(Fy)];
figure(1)
plot(x_f,Fy,'-',x_f,Fc8o12,'.',x_f,Fc8o24,'-');
%axis([0 fs/2 -100 5]);
xlabel("频率（HZ）");ylabel("幅度（dB）");
legend("理想输出",'8比特系数及12比特输出量化结果',"8比特系数24比特输出量化结果");
grid;
figure(2)
plot(x_f,Fy,'-',x_f,Fc12o12,'.',x_f,Fc12o24,'-');
%axis([0 fs/2 -100 5]);
xlabel("频率（HZ）");ylabel("幅度（dB）");
legend("理想输出",'12比特系数及12比特输出量化结果',"比特系数24比特输出量化结果");
grid;

My_E5_32_QuantIIRDirectArith(Qb,Qa,din,Qcoe,Qout)%为自编的根据系数及输出数据位数的量化函数

function Qy=My_E5_32_QuantIIRDirectArith(Qb,Qa,din,Qcoe,Qout);
%直接型IIR滤波器系数及运算字长量化仿真
%Qb=直接型IIR滤波器分子系数向量
%Qa=直接型IIR滤波器分母系数向量
%din=输入原始数据
%Qcoe=IIR滤波器系数的量化位数
%Qout=IIR滤波器的输出数据量化位数
%Qy=IIR滤波器量化输出

%初始化输出
Qy = zeros(1,length(din));
N = length(Qb); %存储Qb长度，即系数个数
%根据直接型IIR滤波器求取量化后的输出
for n=1:length(din)     %逐个遍历输入数据，计算相应的输出
    y = 0;  %一个输入数据对应的输出
    for k=1:N  %计算公式中每个x(k)、y(k)对应的结果

        if n>=k         %Rb储存对应的x(n-k+1);分别从din(n)被赋值到din(n-N+1)
            Rb = din(n-k+1); %正好把公式中的x(n)-x(n-N+1)都利用到了
        else %此时k>n 按照IIR公式，x(n)应该乘以Qb(1),x(n-L)应该乘以Qb(1+L),
             %假设length(Qb)=8，此时正好k=6,当计算到第5个输入数据时，即n=5,此时n<k
             %但是当输入只有5个数据时，顶多使用到Qb(1)-Qb(5),之前循环中的if条件中已经利用完了，
             %此时已经不需要再计算了，所以把Rb=0,相当于取消这次计算
             %只有当输入的输入数据量n大于等于系数长度的时候length(Qb),此时x(n)前面length(Qb)-1
             %个输入数据都用得到,才不会跳入到这个else中
            Rb = 0; 
        end

        if n>=k+1       %Ra储存对应的Qy(n-k+1);分别从Qy(n-1)被赋值到Qy(n-N)
            Ra = Qy(n-k);%正好把公式右边的y(n-1)-y(n-N)都利用到了
        else %道理同上
            Ra = 0;
        end
        %上面的两个if语言几乎同步，一般情况下，n>N时
        %在一次循环中，Rb=x(n-L)，则Ra=y(n-L-1)
        %不同的是第二个if语句没有考虑y(n)。

        %知道x(n-L)、y(n-L-1)后，就可以先带入到公式中进行计算了。
        %因为k=1-N时 Rb=x(n)-x(n-N+1)、Ra=y(n-1)-y(n-N),
        %所以Qb=Qb(1)-Qb(length(Qb))、Qa=[Qa(2)-Qb(N),0]
        if k==N     %此时Ra=y(n-N),但是公式中只有顶多计算到y(n-N-1)
            y = y + Rb*Qb(k) - Ra*0;%所以让Ra乘以0，相当于不计算这一项
        else 
            y = y + Rb*Qb(k) - Ra*Qa(k+1);
        end
    end  %一个输入数据的IIR公式计算完毕，此时得到的结果y=Qa(1)*y(n)
    Qy(n) = y/Qa(1);
    Qy(n)=floor(Qy(n)*(2^(Qout-1)-1)) / (2^(Qout-1)-1);   %直接截尾 量化误差
end

3.结果展示

图1

图2

  对此结果树上的分析为：上述仿真程序没有考虑到输入数据的量化位数。从仿真结果来看，系数及运算字长均会对滤波器的性能产生影响。进行8比特系数量化后，系统的滤波性能已受到较大的影响；进行12比特系数量化后，12比特系数及24比特输出量化结果相差不大；12比特系数及24比特系数输出量化后，IIR滤波器的阻带衰减比无限精度运算（理想条件）恶化了10dB。
  需要注意的是，从图中可以看出，输出量化位数对IIR滤波器性能影响不大。仿真程序My_E5_32_QuantIIRDirectArith(Qb,Qa,din,Qcoe,Qout)中对输出数据量化的过程，是在进行一次完整的滤波输出运算后进行量化处理的。在进行滤波运算过程中采用的是无限精度运算方式，中间并没有量化误差，这种情况下，只要输出数据量化位数足够长，就可以完全表示输出数据的范围，则可以有效减少对滤波器性能的影响。IIR滤波器输出的范围由输入数据及滤波器系数决定，因此可以通过仿真对比输入、输出数据的有效范围来确定IIR滤波器输出数据的位数。对于本例涉及到的IIR滤波器，很容易仿真获取输出数据的范围小于等于输入数据的范围，因此输出数据的位数只要与输入数据位数相同既满足要求。

4.总结

  程序中有一处直接截尾没怎么搞懂意思。

Qy(n)=floor(Qy(n)*(2^(Qout-1)-1)) / (2^(Qout-1)-1);   %直接截尾 量化误差

  问了Chatgpt，回答如下：在数字滤波器中，输出信号是由输入信号和滤波器的系数经过运算得到的。由于数字信号是以离散时间和离散幅度的形式表示的，因此在数字滤波器中，输入信号和滤波器的系数都需要进行量化，即用离散的数值来表示。而滤波器输出的量化指的是将滤波器输出信号进行离散化表示的过程。输出信号的幅度和相位可能会受到量化误差的影响，这些误差可能会影响滤波器的性能和精度。因此，在进行滤波器输出时，需要将输出信号进行量化，以确保输出信号的精度和性能。具体地说，在滤波器输出时，需要对输出信号进行离散化表示，并使用离散的数值来表示输出信号的幅度和相位。为了实现这一目的，通常使用定点数表示输出信号，并将输出信号的幅度和相位四舍五入到最接近的量化级别。这样，输出信号就可以在指定的量化位数内保持一定的精度，同时避免因量化误差导致的输出偏差。

  结合上述内容，并带入过几个值，算是有了初步认识，不过还是认识浅薄。先放在这，慢慢思考