实验五 IIR数字滤波器设计及实现

一、实验目的

? 熟悉用脉冲响应不变法和双线性变换法设计IIR数字滤波器的原理与方法。

? 学会调用MATLAB信号处理工具箱中滤波器设计函数设计IIR数字滤波器，学会根据滤波需求确定滤波器指标参数。

二、实验原理和案例分析

1．模拟低通滤波器设计

模拟低通滤波器的逼近有巴特沃思型、切比雪夫型、和考尔型，分别用如下的函数实现。

[Z,P,K]=buttap(n)；返回一个n阶、巴特沃思型归一化的模拟低通滤波器的零极点增益模型

[Z,P,K]=cheb1ap(n,Rp)；n阶、通带内的最大衰减Rp、切比雪夫Ⅰ型 [Z,P,K]=cheb2ap(n,Rs)；n阶、阻带内的最小衰减Rs、切比雪夫Ⅱ型 [Z,P,K]=elliap(n,Rp,Rs)；n阶、通带内的最大衰减Rp、阻带内的最小衰减Rs、考尔型

2．模拟低通滤波器阶数n的选择函数

滤波器阶数的选择在整个滤波器的设计中占有十分重要的地位和作用。根据需要选择合适的滤波器阶数，MATLAB工具箱中提供了对应于各类模拟低通滤波器的阶数选择函数，如巴特沃思型的buttord、切比雪夫型的cheb1ord、cheb2ord和考尔型的ellipord。这些函数的调用格式大同小异，仅以buttord为例加以说明。

[n,Wn]=buttord(Wp,Ws,Rp,Rs,’s’)

输入参数：Wp通带截止频率，Ws阻带截止频率，Rp通带最大衰减，Rs阻带最小衰减。

输出参数：n为符合要求的滤波器最小阶数，Wn为巴特沃思型模拟低通滤波器3dB截止频率。 ’s’：表示模拟域。

3．零极点增益模型到传递函数模型的转换 [num,den]=zp2tf(Z,P,K)

输入参数：Z，P，K分别表示系统的零极点增益模型的零点、极点和增益； 输出参数：num，den分别为同一系统传递函数模型的分子和分母多项式系数。 4．模拟域的频率变换

将归一化的模拟低通滤波器转换成所需要类型（低通、高通、带通和带阻）的模拟滤波器，可分别用如下命令实现：

[b,a]=lp2lp(Bap,Aap,Wn)；把传递函数形式的归一化模拟低通滤波器原型转换成3dB截止频率为Wn的同型低通滤波器。

[b,a]=lp2hp(Bap,Aap,Wn)；转换成高通

[b,a]=lp2bp(Bap,Aap,W0,Bw)；转换成带通，W0：中心频率，Bw：带宽

[b,a]=lp2bs(Bap,Aap,W0,Bw)；转换成带阻 5．模拟滤波器数字化

[bz,az]=bilinear(b,a,Fs)：采用双线性变换法的映射关系。其中，Fs是采样频率。

[bz,az]=impinvar(b,a,Fs)：采用冲击响应不变法的映射关系。

【实例5-1】设计模拟巴特沃斯低通滤波器，Wp=300Hz，Rp=1dB，Ws=800Hz，Rs=20dB。

%模拟低通滤波器技术指标 Rp=1; Rs=20; fp=300; fs=800; Wap=2*pi*fp; Was=2*pi*fs;

%设计模拟巴特沃斯低通滤波器 [N,Wac]=buttord(Wap,Was,Rp,Rs,'s')

[b,a]=butter(N,Wac,'s') %设计模拟巴特沃斯低通滤波器，Wap为通带边界频率，

rad/s；Was为阻带边界频率，rad/s;Rp为通带最大衰减，dB；Rs为阻带最小衰减，dB；'s'表示为模拟滤波器；函数返回值N为模拟滤波器的最小阶数；Wac为模拟滤波器的截止频率(-3dB频率)，rad/s；b、a分别为模拟滤波器的系统函数分子和分母多项式系数向量；

[H,W]=freqs(b,a); %求模拟滤波器的频率响应 %绘制频响幅度谱

plot(W/2/pi,20*log10(abs(H))); %横轴为频率，单位：HZ;纵轴频响幅度，单位：dB

axis([0,1500,-50,0]);

xlabel('频率Hz');ylabel(' H幅值dB');

【实例5-2】设计模拟巴特沃斯高通滤波器，Wp=800Hz，Rp=1dB，Ws=300Hz，Rs=20dB。

%模拟高通滤波器技术指标 Rp=1; Rs=20; fp=800; fs=300; Wap=2*pi*fp; Was=2*pi*fs;

[N,Wc]=buttord(Wap,Was,Rp,Rs,'s')

[b,a]=butter(N,Wc,'high','s') %设计模拟巴特沃斯高通滤波器，Ωp为通带边界频率，rad/s；Ωs为阻带边界频率，rad/s；αp为通带最大衰减，dB；αs为阻带最小衰减，dB；'s'表示为模拟滤波器；函数返回值N为模拟滤波器的最小阶数；Ωc为模拟滤波器的截止频率(-3dB频率)，rad/s；b、a分别为模拟高通滤波器的系统函数分子和分母多项式系数向量； [H,W]=freqs(b,a);

plot(W/2/pi,20*log10(abs(H))); axis([0,1000,-50,0]);

xlabel('频率Hz');ylabel(' H幅值dB');