实验5 脉冲响应不变法设计IIR数字滤波器

1、 实验5 脉冲响应不变法设计IIR数字滤波器一、实验目的1.掌握利用脉冲响应不变法设计IIR数字滤波器的原理及具体方法。2.加深理解数字滤波器和模拟滤波器之间的技术指标转化。3.掌握脉冲响应不变法设计IIR数字滤波器的优缺点及适用范围。二、实验设备及环境计算机、matlab软件环境。三、实验基础理论1.基本原理从时域响应出发，是数字录波器的单位脉冲响应h(n)模仿模拟滤波器的单位冲击响应ha(t)，h(n)等于ha(t)的取样值。 2变换方法Has拉氏反变换hat时域采样hanT=h(n)z变换Hz将Ha(s)进行部分分式展开Has=k=1NAks-pk (2)对Ha(s)进行拉式变换 hat

2、=k=1NAkepkut(3)对ha(t)时域采样得到h(n)hn=k=1NAkepknTu(nT)=k=1NAkepkTnun(4)对h(n)进行z变换Hz=k=1NAk1-epkz-13.设计步骤 (1)确定数字滤波器性能指标p，st，Rp和As。 (2)将数字滤波器频率指标转换成相应的模拟滤波器频率指标p=pT st=stT (3)根据指标p，st，Rp和As设计模拟滤波器Ha(s)。 (4)将Ha(s)展成部分分式形式Has=k=1NAks-pk (5)将模拟极点Pk转换成数字极点ePkT，得到数字滤波器Has=k=1NAk1-epkTz-1可见Ha(s)至H(z)间的变换关系为1s-

3、sk<=>11-eskTz-1=zz-eskT在MATLAB中有两种方法可以实现上述变换。方法1：利用residue函数和residuez函数实现脉冲响应不变法，这两个函数的使用方法如下r,p,k=residue(b,a)b,a=reisdue(r,p,k)实现多项式形式Hs=bMsM+bM-1sM-1+b0aNsN+aN-1sN-1+a0和部分分式形式Hs=r1s-p1+r2s-p2+rNs-pN+ks之间的转换。r,p,k=residuez(b,a)b,a=reisduez(r,p,k)实现多项式形式Hz=b0+b1z-1+bMz-Ma0+a1z-1+aNz-N和部分分式形式H

4、z=r11-p1z-1+r21-p2z-1+rk1-pkz-1+k1+k2z-1+之间的转换。方法2：MATLAB提供了impinvar函数采用脉冲响应不变法实现模拟滤波器到数字滤波器的变换，其使用方法如下bz,az=impinvar(b,a,fs)采用脉冲响应不变法将模拟滤波器系统函数的系数向量b和a变换成为数字滤波器系统函数的系统向量bz和az，fs为采样频率。bz,az=impinvar(b,a)采样频率默认为1的情况下，采用脉冲响应不变法将模拟滤波器变换为数字滤波器。四、实验内容1.设采用率为fs=4kHz,采用脉冲响应不变法设计一个三阶巴特沃斯数字低通滤波器，其3dB截止频率为fc=

5、1kHz。实验代码及结果：>> b, a = butter(3, 0.5*pi, 's')b = 0 0 0 3.8758a = 1.0000 3.1416 4.9348 3.8758>> bz, az = impinvar(b, a)bz = 0 0.5813 0.2114 0az = 1.0000 -0.3984 0.2475 -0.0432>> w = 0:500*pi/500;>> H, w = freqz(bz, az);>> subplot(211);>> plot(w/pi, abs(H);&

6、gt;> xlabel('omega(pi)');>> ylabel('|H(ejomega)|');>> subplot(212);>> plot(w/pi, angle(H)/pi);>> xlabel('omega(pi)');>> ylabel('Phase of H(ejomega)(pi)');2.设采用频率为fs=10kHz，设计数字低通滤波器，满足如下指标：通带截止频率：fp=1kHz，通带波动：Rp=1dB阻带截止频率：fst=1.5kHz,阻带衰

7、减：As=15dB要求分别采用巴特沃斯、切比雪夫1型、切比雪夫2型和椭圆模拟原型滤波器及脉冲响应不变法进行设计。结合实验结果，分别讨论采用上述方法设计的数字滤波器是否都能满足给定指标要求，分析脉冲响应不变法设计IIR数字滤波器的优缺点及适用范围。实验代码及结果：Butterworth滤波器：>> N, Wn = buttord(0.2*pi, 0.3*pi, 1, 15, 's')N = 6Wn = 0.7087>> b, a = butter(N, Wn, 's')b = 0 0 0 0 0 0 0.1266a = 1.0000 2.7

8、380 3.7484 3.2533 1.8824 0.6905 0.1266>> bz, az = impinvar(b, a)bz = -0.0000 0.0007 0.0105 0.0167 0.0042 0.0001 0az = 1.0000 -3.3443 5.0183 -4.2190 2.0725 -0.5600 0.0647>> w = 0:500*pi/500;H, w = freqz(bz, az);subplot(211);plot(w/pi, abs(H);xlabel('omega(pi)');ylabel('|H(ejom

9、ega)|');subplot(212);plot(w/pi, angle(H)/pi);xlabel('omega(pi)');ylabel('Phase of H(ejomega)(pi)');(2) Chebyshev I型>> N, Wn = cheb1ord(0.2*pi, 0.3*pi, 1, 15, 's')N = 4Wn = 0.6283>> b, a = cheby1(N, 1, Wn, 's')b = 0 0 0 0 0.0383a = 1.0000 0.5987 0.5740

10、0.1842 0.0430>> bz, az = impinvar(b, a)bz = -0.0000 0.0054 0.0181 0.0040 0az = 1.0000 -3.0591 3.8323 -2.2919 0.5495>> w = 0:500*pi/500;>> H, w = freqz(bz, az);>> subplot(211);>> plot(w/pi, abs(H);>> grid on;>> xlabel('omega(pi)');>> ylabel('

11、;|H(ejomega)|');>> subplot(212);>> plot(w/pi, angle(H)/pi);>> grid on;>> xlabel('omega(pi)');>> ylabel('Phase of H(ejomega)(pi)');（3）、Chebyshev II 型：>> N, Wn = cheb2ord(0.2*pi, 0.3*pi, 1, 15, 's')N = 4Wn = 0.8236>> b, a = cheby2(N

12、, 15, Wn, 's')b = 0.1778 0.0000 0.9651 0.0000 0.6547a = 1.0000 2.0708 2.3157 1.3298 0.6547>> bz, az = impinvar(b, a)bz = -0.1904 0.5912 -0.5380 0.2775 0.0224az = 1.0000 -1.9388 1.7388 -0.6968 0.1261>> w = 0:500*pi/500;>> H, w = freqz(bz, az);>> subplot(211);>> p

13、lot(w/pi, abs(H);>> axis tight;>> grid on;>> xlabel('omega(pi)');>> ylabel('|H(ejomega)|');>> subplot(212);>> plot(w/pi, angle(H)/pi);>> axis tight;>> grid on;>> xlabel('omega(pi)');>> ylabel('Phase of H(ejomega)

14、(pi)');（4）椭圆滤波器>> N, Wn = ellipord(0.2*pi, 0.3*pi, 1, 15, 's')N = 3Wn = 0.6283>> b, a = ellip(N, 1, 15, Wn, 's')b = 0 0.2960 -0.0000 0.1873a = 1.0000 0.6161 0.4837 0.1873>> bz, az = impinvar(b, a)bz = 0.2960 -0.4461 0.3063 0az = 1.0000 -2.1192 1.7935 -0.5401>

15、> w = 0:500*pi/500;>> H, w = freqz(bz, az);>> subplot(211);>> plot(w/pi, abs(H);>> axis tight;>> grid on;>> xlabel('omega(pi)');>> ylabel('|H(ejomega)|');>> subplot(212);>> plot(w/pi, angle(H)/pi);>> axis tight;>> grid on;>> xlabel('omega(pi)');>> yla