数字信号处理仿真实验六

6 IIR数字滤波器的设计

6.1 实验目的

1、掌握脉冲响应不变法和双线性变换法设计IIR数字滤波器的原理和方法。

2、观察双线性变换法和脉冲响应不变法设计的滤波器的频率特性，了解双线性变换法和脉冲响应不变法的特点和区别。

6.2 实验原理

6.2.1 脉冲响应不变法

所谓脉冲响应不变法就是使数字滤波器的单位脉冲响应序列h(n)等于模拟滤波器的单位冲激响应和(t)的采样值，即：，其中，T为采样周期。

在MATLAB中，可用函数impinvar实现从模拟滤波器到数字滤波器的脉冲响应不变映射，调用格式为：

[b,a]=impinvar(c,d,fs)

[b,a]=impinvar(c,d)

其中，c、d分别为模拟滤波器的分子和分母多项式系数向量；fs为采样频率（Hz），缺省值fs=1Hz；b、a分别为数字滤波器分子和分母多项式系数向量。

6.2.2 双线性变换法：

由于s平面和z平面的单值双线性映射关系为s＝，其中T为采样周期。

因此，若已知模拟滤波器的传递函数，将上式代入即可得到数字滤波器的系统函数H(z)。在双线性变换中，模拟角频率和数字角频率的变换关系为：

可见，Ω和w之间的变换关系为非线性的。

在MATLAB中，可用函数bilinear实现从模拟滤波器到数字滤波器的双线性变换映射，调用格式为：[b,a]=bilinear(c,d,fs)

6.2.3 设计步骤：

(1) 定技术指标转换为模拟低通原型设计性能指标。

(2) 估计满足性能指标的模拟低通性能阶数和截止频率。

利用MATLAB中buttord、cheb1ord、cheb2ord、ellipord等函数，调用格式如：

其中，Wp为带通边界频率，rad/s；Ws为阻带边界频率，rad/s；Rp为带通波动，dB；Rs为阻带衰减，dB；‘s’表示为模拟滤波器；函数返回值n为模拟滤波器的最小阶数；Wn为模拟滤波器的截止频率(-3dB频率)，rad/s。函数适用低通、高通、带通、带阻滤波器。

(3) 设计模拟低通原型。利用MATLAB中buttap、cheb1ap、cheb2ap、elliap等函数，调用格式如[z,p,k]=buttap(n)。采用上述函数所得到原型滤波器的传递函数为零点、极点、增益表达式，需要和函数[c,d]=zp2tf(z,p,k)配合使用，以转化为多项式形式。

(4) 由模拟低通原型经频率变换获得模拟低通、高通、带通或带阻滤波器。利用MATLAB中lp2lp、lp2hp、lp2bp、lp2bs等函数，调用格式如[c1,d1]=lp2lp(c,d,Wn)。

(5) 利用脉冲响应不变法或双线性不变法，实现模拟滤波器到数字滤波器的映射。本实验用到的特殊函数：[db,mag,pha,grd,w]=freqz_m(b,a)，计算幅频和相频响应。

6.3 实验内容以及步骤

用双线性变换法设计一个Chebyshev1型数字带通滤波器，设计指标为T=1ms,Rp=1dB,Wp1=0.35π，Wp2=0.65π,Rs=60dB,Ws1=0.2πWs2=0.8π。按实验步骤附上所设计滤波器的H（z）及相应的幅频特性曲线定性分析得到的图形，判断设计是否满足要求。

Ts=0.001;

Fs=1/Ts;

Rp=1;Rs=60;

wp1=0.35*pi;

wp2=0.60*pi;

ws1=0.2*pi;

ws2=0.8*pi;

Wp1=(2/Ts)*tan(wp1/2);

Wp2=(2/Ts)*tan(wp2/2);

Wp=[Wp1,Wp2];

Ws1=(2/Ts)*tan(ws1/2);

Ws2=(2/Ts)*tan(ws2/2);

Ws=[Ws1,Ws2];

BW=Wp2-Wp1;

Omegaw0=sqrt(Wp1*Wp2);

[N,OmegaC]=cheb1ord(Wp,Ws,Rp,Rs,'s');

[z0,p0,k0]=cheb1ap(N,Rp);

AnalogB=k0*real(poly(z0));

AnalogA=real(poly(p0));

[BandB,BandA]=lp2bp(AnalogB,AnalogA,Omegaw0,BW);

[DigitalB,DigitalA]=bilinear(BandB,BandA,Fs);

[sos,G]=tf2sos(DigitalB,DigitalA);

[Hz,Wz]=freqz(DigitalB,DigitalA,1024,'whole');

dbHz=20*log10((abs(Hz)+eps)/max(abs(Hz)));

grd=grpdelay(DigitalB,DigitalA,Wz);

subplot(2,3,1);

plot(Wz/pi,abs(Hz));

title('幅频响应');

xlabel(''),ylabel('幅度:|Hz|');

axis([0,1,0,1.1]);

set(gca,'XTickMode','manual','XTick',[0,ws1/pi,wp1/pi,wp2/pi,ws2/pi,1]);

grid;

subplot(2,3,4);

plot(Wz/pi,dbHz);

title('模值(dB)');

xlabel('频率（单位：\pi）');

ylabel('分贝（dB）');

axis([0,1,-40,5]);

set(gca,'YTickmode','manual','YTick',[-50,-30,-2,0]);

set(gca,'YTickLabelMode','manual','YTickLabels',['50';'30';'?2';'?0']);

grid;

subplot(2,3,2);

plot(Wz/pi,angle(Hz)/pi);

title('相频响应');

xlabel('');

ylabel('单位：\pi');

axis([0,1,-1,1]);

grid;

subplot(2,3,5);

title('零极点图');

ylabel('单位：\dB');

xlabel('单位：\pi');

zplane(DigitalB,DigitalA);

axis([-1.1,1.1,-1.1,1.1]);

subplot(2,3,3);

plot(Wz/pi,grd);

title('群延迟');

ylabel('样本');

axis([0,1,0,8]);

set(gca,'YTickmode','manual','YTick',[0:0.5:10]);

grid;

set(gcf,'color','w');

6.4 实验结果分析