实验五 IIR滤波器设计
一、实验目的
1. 掌握冲激响应不变法和双线性变换法设计IIR数字滤波器的原理和方法；
2. 观察双线性变换法和冲激响应不变法设计的滤波器的频域特性，了解双线性变换法和冲激响应不变法的特点和区别。
二、实验仪器
1. 计算机
2. MATLAB软件
三、预习要求
理解冲激响应不变法，双线性变换法，掌握常用 matlab 程序。
四、实验原理
1. IIR数字滤波器一般为线性移不变的因果离散系统，N阶IIR滤波器的传递函数为：
，式中系数
至少有一个不为0。
IIR数字滤波器的设计通常利用模拟滤波器作为原型滤波器，直接由模拟滤波器的频率响经过冲激响应不变法或双线性变换法转换成IIR数字滤波器。
2. 常用模拟原型滤波器
（1）巴特沃斯滤波器
巴特沃斯滤波器通带和阻带都单调衰减，所有结束的幅度函数-3dB点为同一点。其幅度响应为：
，式中N为滤波器阶数，
通带截止频率。
（2）切比雪夫Ⅰ型滤波器
切比雪夫Ⅰ型滤波器在通带呈现等波纹特性，阻带单调衰减，其幅度响应为：
，式中N为滤波器阶数；
，表示通带波纹大小，
越大，波纹越大；
为截止频率，不一定为3dB带宽；
为N阶Chebyshev多项式。
（3）切比雪夫Ⅱ型滤波器
切比雪夫Ⅱ型滤波器在阻带呈现等波纹特性，通带单调衰减，其幅度响应为：
，式中
是阻带衰减达到一定数值时的最低频率。
（4）椭圆滤波器
椭圆滤波器在通带和阻带都呈现波纹特性，在带内均匀波动，具有最快的滚降。
，式中
为椭圆函数。
3. 冲激响应不变法
所谓冲激响应不变法就是使数字滤波器的单位冲激响应序列
等于模拟滤波器的单位冲激响应和
的采样值，即：
，其中，T为采样周期。
令
为模拟系统传递函数，且
，则冲激响应不变法得到的数字滤波器传递函数
。
冲激响应不变法的特点是：
（1）时域逼近特性良好；
（2）模拟频率
和数字频率
呈线性关系
。
（3）存在频率混叠效应，故只适用于限带的模拟滤波器。
在 MATLAB 中，可用函数 impinvar 实现从模拟滤波器到数字滤波器的冲激响应不变映射，调用格式为：
[bz,az]=impinvar(b,a,fs)
[bz,az]=impinvar(b,a)
其中，b、a 分别为模拟滤波器的分子和分母多项式系数向量；fs为采样频率（Hz），缺省值 fs=1Hz；bz、az分别为数字滤波器分子和分母多项式系数向量。
举例：若Ha(s)=
，则b=[1,1];a=[1,5,6]。
4. 双线性变换法
由于
平面和
平面的单值双线性映射关系为
，其中 T 为采样周期。因此，若已知模拟滤波器的传递函数，将上式代入即可得到数字滤波器的系统函数
。
在双线性变换中，模拟角频率和数字角频率的变换关系为：
，可见，模拟频率
和数字频率
之间的变换关系为非线性的。
双线性变换法的特点是：
（1）由于映射是单值变换关系，避免了冲激响应不变法的频率混叠现象；
（2）当频率
增加时，变换关系存在着严重的非线性。
在 MATLAB 中，可用函数 bilinear 实现从模拟滤波器到数字滤波器的双线性变换映射，调用格式为：[bz,az]=bilinear(b,a,fs)
5. 数字滤波器设计步骤
（1）将技术指标转换为模拟低通原型设计性能指标。
（2）估计满足性能指标的模拟低通性能阶数和截止频率。
利用 MATLAB 中 buttord、cheb1ord、cheb2ord、ellipord 等函数，调用格式如：
其中，
和
分别是滤波器的通带截止频率
和阻止截止频率
，单位为
；
和
分别是通带最大衰减系数
和阻带最小衰减系数
，单位为dB；‘s’表示为模拟滤波器；函数返回值 n 为模拟滤波器的最小阶数；
为模拟滤波器的截止频率(-3dB 频率)，rad/s。函数适用低通、高通、带通、带阻滤波器。
（3）设计模拟低通原型滤波器。利用 MATLAB 中 buttap、cheb1ap、
cheb2ap、elliap等函数，调用格式为：[z,p,k]=buttap(n)。
采用上述函数所得到原型滤波器的传递函数为零点、极点、增益表达式，需要和函数[b,a]=zp2tf(z,p,k)配合使用，以转化为多项式形式。
（4）由模拟低通原型经频率变换获得模拟低通、高通、带通或带阻滤波器。
利用 MATLAB 中 lp2lp、lp2hp、lp2bp、lp2bs 等函数，调用格式如：
[c1,d1]=lp2lp(c,d,Wn)。
（5）利用冲激响应不变法或双线性不变法，实现模拟滤波器到数字滤波器的映射。
五、实验内容
1. 利用冲激响应不变法，设计一数字低通滤波器，并画出各滤波器的幅频和相频特性曲线。滤波器指标
、
、
、
，采样频率fs=1kHz，要求：
function [db,mag,pha,w]=freqz_m(b,a);%对FIR滤波器，b=h，a=[1]
[H,w]=freqz(b,a,1000,'whole');
H=(H(1:1:501))';
w=(w(1:1:501))';
mag=abs(H);
db=20*log10((mag+eps)/max(mag));
