Matlab——IIR滤波器系数的量化

lit_sang

已于 2022-07-20 20:08:34 修改

阅读量2.3k

点赞数

文章标签： matlab 算法

于 2022-07-20 20:00:06 首次发布

本文链接：https://blog.csdn.net/lit_sang/article/details/125899110

版权

1 基本概念

1.1 直接型

直接型结构是根据系统函数 $H\left ( z \right )$ 直接得到，并没有经过任何的重新排列。N阶的IIR系统的差分方程为

$y\left ( n \right )=\sum_{i=0}^{M}b_{i}x\left ( n-i \right )+\sum_{i=1}^{N}a_{i}y\left ( n-i \right )$

其对应的系统函数为

$H\left ( z \right )=\frac{\sum_{i=0}^{M}b_{i}z^{-i}}{1-\sum_{i=0}^{N}a_{i}z^{-i}}$

1.2 级联型

假设 $N\geq M$ ，利用因式分解，可以将系统函数分解成多个子系统级联的形式，即将 $H\left ( z \right )$ 表示为

$H\left ( z \right )=H_{1}\left ( z \right )H_{2}\left ( z \right )...H_{k}\left ( z \right )=\prod_{k=1}^{K}H_{k}\left ( z \right )$

其中， $H_{k}\left ( z \right )$ 是一阶或二阶子系统。二阶子系统 $H_{k}\left ( z \right )$ 一般形式为

$H_{k}\left ( z \right )=\frac{b_{k0}+b_{k1}z^{-1}+b_{k2}z^{-2}}{1-a_{k1}z^{-1}-a_{k2}z^{-2}}$ k=1,2...K

2 Matlab量化系数

2.1 量化系数前后，零极点图的对比

n=16;
fs=2000;
[b,a]=ellip(4,0.2,40,[300 500]/(fs/2));
figure(1),zplane(b,a);xlabel('原系统'); %原始未量化滤波器
%直接型系数的量化
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;
bd = round(b/Qm*(2^(n-1)-1))
ad = round(a/Qm*(2^(n-1)-1))
figure(2),zplane(bd,ad);xlabel('直接型');%直接型系数的量化

[bc,ac]=df2cf(b,a); bc=real(bc);ac=real(ac); %分解为级联型确保系数为实数
bc16=round(bc/Qm*(2^(n-1)-1));ac16=round(ac/Qm*(2^(n-1)-1));%对直接型系数的量化
bc1=[0 0 0];bc2=bc1;bc3=bc1;bc4=bc1;ac1=bc1;ac2=bc1;ac3=bc1;ac4=bc1;
for i=1:3 
    bc1(i)=bc16(1,i);
    bc2(i)=bc16(2,i);
    bc3(i)=bc16(3,i);
    bc4(i)=bc16(4,i); 
    ac1(i)=ac16(1,i);
    ac2(i)=ac16(2,i);
    ac3(i)=ac16(3,i);
    ac4(i)=ac16(4,i);
end;
bc22=conv(bc1,bc2);bc23=conv(bc22,bc3);bc=conv(bc23,bc4); %重新组合系数
ac22=conv(ac1,ac2);ac23=conv(ac22,ac3);ac=conv(ac23,ac4);
figure(3),zplane(bc,ac);xlabel('级联型') ;

function [bc,ac]=df2cf(b,a)
N=length(b)-1;
z=roots(b)
p=roots(a)
for k=1:N/2
bc(k,:)=poly(z(2*k-1:2*k)); 
ac(k,:)=poly(p(2*k-1:2*k));
end
bc(1,:)=b(1)*bc(1,:);

n — Filter order

Filter order, specified as an integer scalar. For bandpass and bandstop designs, n represents one-half the filter order.

r = roots(p) 以列向量的形式返回 p 表示的多项式的根。输入 p 是一个包含 n+1 多项式系数的向量，以 xn 系数开头。0 系数表示方程中不存在的中间幂。例如：p = [3 2 -2] 表示多项式 3x2+2x−2。

p = poly(r)（其中 r 是向量）返回多项式的系数，其中多项式的根是 r 的元素。

w = conv(u,v) 返回向量 u 和 v 的卷积。如果 u 和 v 是多项式系数的向量，对其卷积与将这两个多项式相乘等效。

2.2 量化系数前后，幅频特性的对比

n=26;
fs=2000;
wp=[693 701]*2/fs;
ws=[685 709]*2/fs;
Rp=0.2;
Rs=60;
[N,Wn]=ellipord(wp,ws,Rp,Rs);
[b,a]=ellip(N,Rp,Rs,Wn,'bandpass');%系统的实现
[H,w]=freqz(b,a,4096);
figure(1),plot((w/pi)*(fs/2),20*log10(abs(H)));%未量化系统的频率特性
xlabel('原系统');axis([0,fs,-80,10]);
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;
bd = round(b/Qm*(2^(n-1)-1))
ad = round(a/Qm*(2^(n-1)-1))
[H1,w]=freqz(bd,ad,4096);%直接型的量化
figure(2),plot(w,20*log10(abs(H1)),'r',w,20*log10(abs(H)),'b');
xlabel('直接型');axis([0,3.14,-80,10]);

[bc,ac]=df2cf(b,a);
bcm= round(bc/Qm*(2^(n-1)-1))
acm= round(ac/Qm*(2^(n-1)-1)) %级联型及其量化
x=[1 zeros(1,4095)];
t=0:4095;%设置单位冲激信号,求单位冲激响应
h1=filter(bcm(1,:),acm(1,:),x);h2=filter(bcm(2,:),acm(2,:),h1);
h3=filter(bcm(3,:),acm(3,:),h2);h4=filter(bcm(4,:),acm(4,:),h3);
h5=filter(bcm(5,:),acm(5,:),h4);
[H2,w]=freqz(h5,1,4096);
figure(3),plot(w,20*log10(abs(H2)),'r',w,20*log10(abs(H)),'b');
xlabel('级联型');axis([0,3.14,-80,10]);