实验四 IIR数字滤波器的设计

电气与信息工程学院
数字信号处理
实验报告
学生姓名
班级电子信息工程
学号
指导教师
2019.12
实验四 IIR 数字滤波器的设计
一、实验目的：
1. 掌握双线性变换法及脉冲响应不变法设计IIR 数字滤波器的具体设计方法及其原理，熟悉用双线性变换法及脉冲响应不变法设计低
通、高通和带通IIR 数字滤波器的MATLAB 编程。

2. 观察双线性变换及脉冲响应不变法设计的滤波器的频域特性，了解双线性变换法及脉冲响应不变法的特点。

3. 熟悉巴特沃思滤波器、切比雪夫滤波器和椭圆滤波器的频率特
性。

二、实验原理:
1． 脉冲响应不变法
用数字滤波器的单位脉冲响应序列 模仿模拟滤波器的冲激响应 ,让 正好等于 的采样值，即 ，其中 为采样间隔，如果以 及 分别表示 的拉式变换及 的Z 变换，则
)2(1)(m T j s H T z H m a e z sT ∑∞-∞==+=π
2．双线性变换法
S 平面与z 平面之间满足以下映射关系：
);(,2121,11211ωωσj re z j s s T s T z z z T s =+=-+
=+-⋅=-- s 平面的虚轴单值地映射于z 平面的单位圆上，s 平面的左半平面完全映射到z 平面的单位圆内。

双线性变换不存在混叠问题。

双线性变换是一种非线性变换，这种非线性引起的幅频特性畸变可通过预畸而得到校正。

三、实验内容及步骤：
实验中有关变量的定义:
fc 通带边界频率；
fr阻带边界频率；
δ通带波动；
At 最小阻带衰减；
fs采样频率；
T采样周期
上机实验内容：
（1）fc=0.3KHz,δ=0.8dB,fr=0.2KHz, At=20dB,T=1ms;
设计一个切比雪夫高通滤波器，观察其通带损耗和阻带衰减是否满足要求。

MATLAB源程序：
wp=2*1000*tan(2*pi*300/(2*1000));
ws=2*1000*tan(2*pi*200/(2*1000));
[N,wn]=cheb1ord(wp,ws,0.8,20,'s'); %给定通带（wp）和阻带(ws)边界角频率，通带波动波动0.8，阻带最小衰减20dB，求出最低阶数和通带滤波器的通带边界频率Wn
[B,A]=cheby1(N,0.5,wn,'high','s');%给定通带（wp）和阻带(ws)边
界角频率，通带波动
[num,den]=bilinear(B,A,1000);
[h,w]=freqz(num,den);
f=w/(2*pi)*1000;
plot(f,20*log10(abs(h)));
axis([0,500,-80,10]);
grid;xlabel('频率');ylabel('幅度/dB')
程序结果
num = 0.0304 -0.1218 0.1827 -0.1218 0.0304 den = 1 1.3834 1.4721 0.8012 0.2286
系统函数：
1234
1234
0.0304 -0.1218z 0.1827z-0.1218z0.0304z H(z)=
1.0000+1.3834z+1.4721z+ 0.8012z+0.2286z
----
----
++
幅频响应图：
分析：由图可知，切比雪夫滤波器幅频响应是通带波纹，阻带单调衰减的。

δ=0.8，fr=0.2kHz,At=30dB,满足设计要求
（2）fc=0.2kHz,δ=1dB,fr=0.3kHz,At=25dB,T=1ms;分别用脉冲响应不变法及双线性变换法设计一Butterworth数字低通滤波器，观察所设计数字滤波器的幅频特性曲线，记录带宽和衰减量，检查是否满足要求。

比较这两种方法的优缺点。

MATLAB源程序：
T = 0.001;fs = 1000;fc = 200;fr = 300;
wp1 = 2*pi*fc;wr1 = 2*pi*fr;
[N1,wn1] = buttord(wp1,wr1,1,25,'s')
[B1,A1] = butter(N1,wn1,'s');
[num1,den1] = impinvar(B1,A1,fs);%脉冲响应不变法
[h1,w] = freqz(num1,den1);
wp2 = 2*fs*tan(2*pi*fc/(2*fs))
wr2 = 2*fs*tan(2*pi*fr/(2*fs))
[N2,wn2] = buttord(wp2,wr2,1,25,'s')
[B2,A2] = butter(N2,wn2,'s');
[num2,den2] = bilinear(B2,A2,fs);%双线性变换法
[h2,w] = freqz(num2,den2);
f = w/(2*pi)*fs;
plot(f,20*log10(abs(h1)),'-.',f,20*log10(abs(h2)),'-'); axis([0,500,-100,10]);grid;xlabel('频率/Hz ');ylabel('幅度/dB')
title('巴特沃思数字低通滤波器');
legend('脉冲相应不变法','双线性变换法',1);
结果分析：
脉冲响应不变法的低通滤波器系统函数:
num1 -2.3647 0.0002 0.0153 0.0995 0.1444 0.0611
0.0075 0.0002 3.6569 0
den1 1 -1.9199 2.5324 -2.2053 1.3868 -0.6309
0.2045 -0.0450 0.0060 -0.0004
z
987654321876543210.0004z 0.0060z 0.0450z 0.2045z 0.6309z 1.3869z 2.2053z 2.5324z 1.9199z 1 3.6569z 0.0002z 0.0075z 0.0611z 0.1444z 0.0995z 0.0153z 0.0002z 2.3647------------------+-+-+-+-++++++++-=)(z H 双线性变换法设计的低通滤波器系统函数:
num2 0.0179 0.1072 0.2681 0.3575 0.2681
0.1072 0.0179
den2 1 -0.6019 0.9130 -0.2989 0.1501 -0.0208 0.0025
6
543216
543210025.00208.01501.02989.09130.06019.010176.01072.02681.03575.02681.00.10720179.0)(------------+-+-+-++++++=z z z z z z z z z z z z z H
分析：脉冲响应不变法的频率变化是线性的，数字滤波器频谱响应出现了混叠，影响了过渡带的衰减特性，并且无传输零点；双线性变化法的频率响应是非线性的，因而消除了频谱混叠，在f=500Hz 出有一个传输零点。

脉冲响应不变法的一个重要特点是频率坐标的变换是线性的，
ω＝Ω
Τ，ω与Ω是线性关系：在某些场合，要求数字滤波器在时域上能模仿模拟滤波器的功能时，如要实现时域冲激响应的模仿，一般使用脉冲响应不变法。

脉冲响应不变法的最大缺点：有频谱周期延拓效应，因此只能用于带限的频响特性，如衰减特性很好的低通或带通,而高频衰减越大，频响的混淆效应越小，至于高通和带阻滤波器,由于它们在高频部分不衰减，因此将完全混淆在低频响应中，此时可增加一保护滤波器，滤掉高于的频带，再用脉冲响应不变法转换为数字滤波器，这会增加设计的复杂性和滤波器阶数，只有在一定要满足频率线性关系或保持网络瞬态响应时才采用。

双线性变换法的主要优点是S平面与Z平面一一单值对应，s平面的虚轴(整个jΩ)对应于Z平面单位圆的一周，S平面的Ω=0处对应于Z平面的ω=0处，Ω= ∞处对应于Z平面的ω= π处,即数字滤波器的频率响应终止于折叠频率处，所以双线性变换不存在混迭效应。

双线性变换缺点: Ω与ω成非线性关系，导致：
a. 数字滤波器的幅频响应相对于模拟滤波器的幅频响应有畸变，(使数字滤波器与模拟滤波器在响应与频率的对应关系上发生畸变)。

b. 线性相位模拟滤波器经双线性变换后，得到的数字滤波器为非线性相位。

c.要求模拟滤波器的幅频响应必须是分段恒定的，故双线性变换只能用于设计低通、高通、带通、带阻等选频滤波器。

(3)利用双线性变换法分别设计满足下列指标的巴特沃思型、切比雪夫型和椭圆型数字低通滤波器，并作图验证设计结果：fc=1.2kHz ，δ≤0.5dB ,fr=2kHz ，At≥40dB, fs=8kHz，比较这种滤波器的阶数。

MATLAB源程序：
clear all;
wc=2*pi*1200;wr=2*pi*2000;rp=0.5;rs=40;fs=8000;
w1=2*fs*tan(wc/(2*fs));
w2=2*fs*tan(wr/(2*fs));
[Nb,wn]=buttord(w1,w2,rp,rs,'s') %巴特沃思
[B,A]=butter(Nb,wn,'s');
[num1,den1]=bilinear(B,A,fs);
[h1,w]=freqz(num1,den1);
[Nc,wn]=cheb1ord(w1,w2,rp,rs,'s') %切比雪夫
[B,A]=cheby1(Nc,rp,wn,'s');
[num2,den2]=bilinear(B,A,fs);
[h2,w]=freqz(num2,den2);
[Ne,wn]=ellipord(w1,w2,rp,rs,'s') %椭圆型
[B,A]=ellip(Ne,rp,rs,wn,'low','s');
[num3,den3]=bilinear(B,A,fs);
[h3,w]=freqz(num3,den3);
f=w/(2*pi)*fs;
plot(f,20*log10(abs(h1)),'-',f,20*log10(abs(h2)),'--
',f,20*log10(abs(h3)),':');
axis([0,4000,-100,10]);grid;
xlabel('Frequency in Hz'); ylabel('Gain in dB');
title('三种数字低通滤波器');
legend('巴特沃思数字低通滤波器','切比雪夫数字低通滤波器','椭圆数字低通滤波器',3);
巴特沃思数字低通滤波器的系统函数系数：
num1= 0.0032 0.0129 0.0302 0.0453 0.0453 0.0302
0.0129 0.0032 0.0003
den1= -2.7996 4.4581 -4.5412 3.2404 -1.6330 0.5780 -0.1370 0.0197 -0.0013
切比雪夫数字低通滤波器的系统函数系数：
num2= 0.0026 0.0132 0.0264 0.0264 0.0132 0.0026 den2= 1 -2.9775 4.2932 -3.5124 1.6145 -0.3334 椭圆数字低通滤波器的系统函数系数：
num3= 0.03887 0.0363 0.0665 0.0363 0.0389
den3= 1 -2.1444 2.3658 -1.3250 0.3332
程序结果图：
分析：设计结果表明，巴特沃思数字低通滤波器、切比雪夫数字低通滤波器、椭圆数字低通滤波器的阶数分别是9、5、4阶。

可见，对于给定的阶数，椭圆数字低通滤波器的阶数最少(换言之，对于给定的阶数，过渡带最窄)，就这一点来说，他是最优滤波器。

由图表明，巴特沃思数字低通滤波器过渡带最宽，幅频响应单调下降；椭圆数字低通滤波器过渡带最窄，并具有等波纹的通带和阻带响应；切比雪夫数字低通滤波器的过渡带介于两者之间。

(4）分别用脉冲响应不变法及双线性变换法设计一巴特沃思型数字带通滤波器，已知，其等效的模拟滤波器指标δ<3dB，2kHz<f≤3kHz ；At≥5dB ，f≥6kHz ；At≥20dB ，f≤1.5kHz 。

MATLAB源程序：
wp1 = 2*pi*2000;wp2 = 2*pi*3000;
ws1= 2*pi*1500;ws2= 2*pi*6000;
[N1,wn1] = buttord([wp1 wp2],[ws1 ws2],3, 20 ,'s');%求巴特沃思滤波器的阶数
[B1,A1] = butter(N1,wn1,'s');%给定阶数和边界频率设计滤波器[num1,den1] = impinvar(B1,A1,30000);%脉冲相应不变法
[h1,w] = freqz(num1,den1);
w1=2*30000*tan(2*pi*2000/(2*30000));
w2=2*30000*tan(2*pi*3000/(2*30000));
wr1=2*30000*tan(2*pi*1500/(2*30000));
wr2=2*30000*tan(2*pi*6000/(2*30000));
[N,wn]=buttord([w1 w2],[wr1 wr2],3,20,'s');%求巴特沃思滤波器的阶数
[B,A]=butter(N,wn,'s');
[num,den]=bilinear(B,A,30000);%双线性变化法
[h2,w]=freqz(num,den);
f=w/(2*pi)*30000;
plot(f,20*log10(abs(h1)),'-.',f,20*log10(abs(h2)),'-'); axis([0,15000,-60,10]);
xlabel('Frequency in Hz'); ylabel('Gain in dB');
grid;
title('巴特沃思数字带通滤波器');
legend('脉冲相应不变法','双线性变换法',1);
脉冲相应不变法设计的巴特沃思数字带通滤波器系统函数的分子、分母多项式系数：
num1= -1.5158 0.0057 -0.0122 0.0025 0.0089 -0.0049 den1= 1 -4.8056 10.2376 -12.2625 8.7012 -3.4719
0.6145
双线性变换法设计的巴特沃思数字带通滤波器系统函数的分子、分母多项式的系数：
num= 0.0014 0 -0.0042 7.1054 0.0042 5.7732 -0.0014 den= 1 -4.8071 10.2473 -12.2838 8.7245 -3.4849 0.6176
（5）利用双线性变换法设计满足下列指标的切比雪夫型数字带阻滤波器，并作图验证设计结果：当kHz f kHz 21≤<时，dB At 18≥；当Hz f 500≤以及kHz f 3≥时，dB 3≤δ；采样频率kHz f s 10=。

MATLAB 源程序：
w1=2*10000*tan(2*pi*1000/(2*10000));
w2=2*10000*tan(2*pi*2000/(2*10000));
wr1=2*10000*tan(2*pi*500/(2*10000));
wr2=2*10000*tan(2*pi*3000/(2*10000));
[N,wn]=cheb1ord([wr1 wr2],[w1 w2],3,18,'s');%计算阶数
[B,A]=cheby1(N,3,wn,'stop','s');%给定阶数和参数设计滤波器
[num,den]=bilinear(B,A,10000);%双线性变化法
[h,w]=freqz(num,den);%频率响应
f=w/(2*pi)*10000;
plot(f,20*log10(abs(h)));
axis([0,5000,-120,10]);
grid;xlabel('频率/Hz');ylabel('幅度/dB')
title('切比雪夫数字带阻滤波器');
程序结果图：
四、实验思考题
1．双线性变换法中Ω和ω之间的关系是非线性的，在实验中你注意到这种非线性关系了吗？从哪几种数字滤波器的幅频特性曲线中可以观察到这种非线性关系？
答：在双线性变化法中，模拟频率与数字频率不再是线性关系，所以一个线性相位模拟滤波器经双线性变换后，得到的数字滤波器不再保持原有的线性相位了，在每一幅使用了双线性变换的图中，可以看到在采样频率一半处，幅度为零，这显然不是线性变换能够产生的，这是由于双线性变换将模拟域中的无穷远点映射到了改点处。

2. 能否利用公式完成脉冲响应不变法的数字滤波器设计？为什么？
答：IIR数字滤波器的设计实际上是求解滤波器的系数,
a b，它是数
k k
学上的一种逼近问题，即在规定意义上（通常采用最小均方误差准则）去逼近系统的特性。

如果在S平面上去逼近，就得到模拟滤波器；如果在z平面上去逼近，就得到数字滤波器。

但是它的缺点是，存在频率混叠效应，故只适用于阻带的模拟滤波器。