数字信号处理实验5

实验5 IIR 数字滤波器的设计
一、实验目的
（1）加深对脉冲响应不变法、双线性变换法的了解；
（2）掌握MATLAB 进行滤波器设计的子函数。

二、知识点提示
本章节的主要知识点是使用脉冲响应不变法、双线性变换法设计IIR DF ；重点是掌握利用脉冲响应不变法、双线性变换法设计数字低通、高通、带通、带阻滤波器的方法；难点是MATLAB 相关子函数的使用。

三、实验涉及的MATLAB 子函数
1.buttord
功能：确定巴特沃斯滤波器的阶数和3dB 截止频率。

调用格式：①[n,wn]=buttord(wp,ws,Rp,As)
其中，1)(0≤≤s p ωω或，其值为1时表示0.5Fs 。

Rp 为通带最大衰减指标，As 为阻带最小衰减指标。

②[n,wn]=buttord(wp,ws,Rp,As,’s ’)
wp 、ws 可以是实际的频率值或角频率值，wn 将取相同的量纲。

Rp 为通带最大衰减指标，As 为阻带最小衰减指标。

当wp>ws 时为高通滤波器，当wp 、ws 为二元向量时，为带通或带阻滤波器，此时wn 也为二元向量。

2. buttap
功能：巴特沃斯模拟滤波器原型
调用格式：[z,p,k]=buttap(n)
3. lp2lp 、lp2hp 、lp2bp 、lp2bs 模拟域频率变换的子函数
功能：低通原型到模拟低通滤波器的变换、低通原型到模拟高通滤波器的变换、低通原型到模拟带通滤波器的变换、低通原型到模拟带阻滤波器的变换。

调用格式：[bt,at]=lp2lp(b,a,w0)
将传递函数表示的截止频率为1 rad/s 的模拟低通滤波器原型变换成截止频率为w0的低通滤波器。

4. butter
功能：巴特沃斯滤波器设计的完全函数。

调用格式：[b,a]=butter(n,wn,’ftype’)
n为滤波器阶数，wn为滤波器截止频率（0~1）。

（在MATLAB滤波器设计工具函数中，数字频率采用标准化频率，取值范围为0~1之间，标准化频率1对应数字频率π，对应的模拟频率为采样频率的一半。

）’ftype’为滤波器类型：’high’为高通，截止频率wn；’stop’为带阻，截止频率wn=[w1,w2]；缺省时为低通和带通滤波器，低通滤波器时wn为截止频率，带通滤波器时wn=[w1,w2]。

5. impinvar
功能：用脉冲响应不变法实现模拟到数字滤波器的变换。

调用格式：[bd,ad]=impinvar(b,a,Fs)
将模拟滤波器系数b,a变换成数字滤波器系数bd,ad，两者冲激响应不变。

Fs缺省值为1Hz。

6. bilinear
功能：双线性变换，将s域映射到z域的标准方法。

调用格式：①[numd,dend]=bilinear(num,den,fs)
将模拟域传递函数变换为数字域传递函数，fs为取样频率。

②[numd,dend]=bilinear(num,den,fs,fp)
将模拟域传递函数变换为数字域传递函数，fs为取样频率，fp为通带截止频率。

四、实验原理
在MATLAB中，经典设计IIR数字滤波器采用下面步骤：
（1）根据给定的性能指标和方法，首先对设计性能指标中的频率指标进行转换，转换后的频率指标作为模拟滤波器原型设计指标；
wp=0.2*pi;ws=0.3*pi;rp=1;as=15;fs=10000;omgp=wp*fs;omgs=ws*fs;
（2）估计模拟滤波器最小阶数和边界频率，可利用MATLAB工具函数buttord，cheblord 等。

[n,omgc]=buttord(omgp,omgs,rp,as,'s');
（3）设计模拟低通滤波器原型，可利用MATLAB工具函数buttap，cheblap等；
[z,p,k]=buttap(n);
（4）由模拟低通原型经频率变换得到模拟滤波器（低通、高通、带通、带阻），可利用MATLAB工具函数lp2lp，lp2hp，lp2bp，lp2bs等；
[bap,aap]=zp2tf(z,p,k);[ba,aa]=lp2lp(bap,aap,omgc);
【或者将上面的两步合并，利用巴特沃斯设计的完全函数butter实现[ba,aa]=butter(n,omgc,'s');】（5）将模拟滤波器离散化得到IIR数字滤波器，可利用MATLAB工具函数bilinear，
impinvar 等。

[bd,ad]=impinvar(ba,aa,fs);
五、实验内容
1. 利用脉冲响应不变法，用巴特沃斯滤波器原型设计一个低通滤波器，满足：采样频率为10KHz ，通带截止频率、通带最大衰减、阻带截止频率和阻带最小衰减分别为
dB 15,π3.0,dB 1,π2.0====s s p p A R ωω
代码：
wp=0.2*pi;
ws=0.3*pi;
rp=1;
as=15;
fs=10000;
omgp=wp*fs;
omgs=ws*fs;
[n,omgc]=buttord(omgp,omgs,rp,as,'s');
[z,p,k]=buttap(n);
[bap,aap]=zp2tf(z,p,k);
[ba,aa]=lp2lp(bap,aap,omgc);
[bd,ad]=impinvar(ba,aa,fs);
[H,wn] = freqz(bd,ad);
plot(wn/pi,20*log10(abs(H)))
grid on
结果：
-15-14
-13
-12
-11
-10
-9
-8
-7
-6
-5
-4
-3
-2
-1
2. 利用双线性变换法设计巴特沃斯高通数字滤波器，满足：通带截止频率为400Hz ，阻带截止频率为200Hz ，通带衰减小于3dB ，阻带衰减大于15dB ，采样频率为1000Hz 。

（选作）
代码：
wp=0.8*pi;
ws=0.4*pi;
rp=3;
as=15;
fs=1000;
omgp=wp*fs;
omgs=ws*fs;
[n,omgc]=buttord(omgp,omgs,rp,as,'s');
[z,p,k]=buttap(n);
[bap,aap]=zp2tf(z,p,k);
[ba,aa]=lp2hp(bap,aap,omgc);
[numd,dend]=bilinear(ba,aa,fs);
[H,wn] = freqz(numd,dend);
plot(wn/pi,20*log10(abs(H)))
grid on
xlabel('频率(Hz)')
ylabel('幅度')
title('数字滤波器幅频曲线')
频率(Hz)幅度数字滤波器幅频曲线
3. 设计巴特沃斯低通数字滤波器，满足：通带截止频率、通带最大衰减、阻带截止频率和阻带最小衰减分别为dB 15,π3.0,dB 1,π2.0====s s p p A R ωω。

代码：
wp=0.2;
ws=0.3;
rp=1;
as=15;
[n,wn]=buttord(wp,ws,rp,as);
[b,a]=butter(n,wn);
[H,wn] = freqz(b,a);
plot(wn/pi,20*log10(abs(H)));
grid on
xlabel('频率(Hz)')
ylabel('幅度')
title('数字滤波器幅频曲线')
结果：
0.10.2
0.30.4
-15-14
-13
-12
-11
-10
-9
-8
-7-6
-5
-4
-3
-2
-1
频率(Hz)幅度数字滤波器幅频曲线
4. 设计巴特沃斯高通数字滤波器，满足：通带截止频率为400Hz ，阻带截止频率为200Hz ，通带衰减小于3dB ，阻带衰减大于15dB ，采样频率为1000Hz 。

代码：
wp=0.8;
ws=0.4;
rp=3;
as=15;
[n,wn]=buttord(wp,ws,rp,as);
[b,a]=butter(n,wn,'high');
[H,wn] = freqz(b,a);
plot(wn/pi,20*log10(abs(H)))
grid on
xlabel('频率(Hz)')
ylabel('幅度')
title('数字滤波器幅频曲线')
结果：
0.30.40.50.60.7
0.80.91
-15-14
-13
-12
-11
-10
-9
-8
-7-6
-5
-4
-3
-2
-1
频率(Hz)幅度数字滤波器幅频曲线
三、思考题
3.1 为什么脉冲响应不变法不能用于设计数字高通滤波器和带阻滤波器？数字滤波器的频率响应与模拟滤波器的频率响应有何区别？
由于频谱的周期延拓会产生混叠失真。

数字滤波器的频率响应是模拟滤波器的频率响应的周期延拓。

3.2 使用双线性变换法时模拟频率和数字频率有何关系？会带来什么影响？如何解决？ 模拟频率和数字频率是非线性关系。

双线性变换法虽然避免了“频率混叠效应”，但出现了模拟频率和数字频率为一种非线性关系的情形。

根据公式修正。