利用脉冲响应不变法设计IIR数字滤波器

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

一、 实验目的

1.掌握利用脉冲响应不变法设计IIR数字滤波器的原理及具体方法。 2.加深理解数字滤波器和模拟滤波器之间的技术指标转化。

3.掌握脉冲响应不变法设计IIR数字滤波器的优缺点及适用范围。

二、 实验设备与环境

计算机、MATLAB软件环境。

三、 实验基础理论

1.基本原理

从时域响应出发，使数字滤波器的单位脉冲响应h(n)模仿模拟滤波器的单位冲击响应ha(t),h(n)等于ha(t)的取样值。 2.变换方法

思 路：拉 式 反 变 换时 域 采 样z变 换????ha(s)??????ha(t)??????ha(nT)?h(n)????H（z)（1）将Ha(s) 进行部分分式展开

Ha(s)??k?1NAks?Pk

（2）对Ha(s)进行拉式反变换

ha(t)??Akepktu(t)

k?1N（3）对ha(t)时域采样得到h(n)

h(n)??Akek?1NpknTu(nT)??Ak(epkT)u(n)

k?1N（4）对h(n)进行Z变换

Akh(z)??pkT?11?ezk?13.设计步骤

N

（1） 确定数字滤波器性能指标fst?1.5kHz

（2） 将数字滤波器频率指标转换成相应的模拟滤波器频率指标

?p??pT

?st??stT

（3） 根据指标?p,?st,Rp,As设计模拟滤波器Ha(s)将Ha(s)展成部分分式形式

Ha(s)??k?1NAks?pk

（4） 把模拟极点pk转换成数字极点epkT，得到数字滤波器

H(z)??AkpkT?1z k?11?eN可见Ha(s)至H(z)间的变换关系为

11z??skT?1s?sk1?ezz?eskT

在MATLAB中有两种方法可以实现上述变换。

方法1：利用residue函数和residuez函数实现脉冲响应不变法，这两个函数的使用方法如下： [r,p,k]=residue(b,a) [b,a]=residue(r,p,k) 实现多项式形式

bMsM?bM?1sM?1???b0H(s)?aNsN?aN?1sN?1??a0

和部分分式形式

rNr1r2H(s)??????k(s)s?p1s?p2s?pN

之间的转换。

[r,p,k]=residuez(b,a) [b,a]=residuez(r,p,k) 实现多项式形式

b0?b1z?1???bMz?MH(z)?a0?a1z?1??aNz?N

和部分分式形式

H(z)?rNr1r2?1?????k?kz??12?1?1?11?p1z1?p2z1?pNz

之间的转换。

方法2：MATLAB提供了impinvar函数采用脉冲响应不变法实现模拟滤波器到数字滤波器的变换，其使用方法如下：

[bz,az]=impinvar(b,a,fs)采用脉冲响应不变法将模拟滤波器系统函数的系数向量b和a转换成为数字滤波器系统函数的系统向量bz和az,fs为采样率。 [bz,az]=impinvar(b,a)采样频率默认为1的情况向下，采用脉冲响应不变法将模拟滤波器变换为数字滤波器。

四、 实验内容

1.设采样频率,采用脉冲响应不变法设计一个三阶巴特沃斯数字低通滤波器，其3dB截止频率。

实验代码与实验结果： N=3; fs=4000; fc=1000; Wc=0.25*pi;

[b,a]=butter(N,Wc,'s'); [bz,az]=impinvar(b,a); w=[0:500]*pi/500; [H,w]=freqz(bz,az); subplot 221; plot(w/pi,abs(H)); xlabel('\\omega(\\pi)');

ylabel('|H(e^j^\\omega)|(dB)'); grid on; subplot 222;

plot(w/pi,20*log10(abs(H))); xlabel('\\omega(\\pi)');

ylabel('|H(e^j^\\omega)|(dB)'); grid on; subplot 223;

plot(w/pi,angle(H)/pi); xlabel('\\omega(\\pi)'); ylabel('|H(e^j^\\omega)|'); grid on; subplot 224;

grd=grpdelay(bz,az,w); plot(w/pi,grd);

xlabel('\\omega(\\pi)'); ylabel('|H(e^j^\\omega)|'); grid on;

1|H(ej?)|(dB)|H(ej?)|(dB)0-10-20-30-4000.5?(?)10.5000.5?(?)110.5|H(ej?)|420-0.5-100.5?(?)1|H(e)|j?0-200.5?(?)12.设采样频率，设计数字低通滤波器，满足如下指标 通带截止频率：，通带波动： 阻带截止频率：，阻带衰减：

要求分别采用巴特沃斯、切比雪夫I型、切比雪夫II型和椭圆模拟原型滤波器及脉冲响应不变法进行设计。结合实验结果，分别讨论采用上述方法设计的数字滤波器是否都能满足给定指标要求，分析脉冲响应不变法设计IIR数字滤波器的优缺点及适用范围。 实验代码与实验结果： （1）巴特沃斯 fs=10000; fp=1000; fst=1500; Wp=0.1*pi; Ws=0.15*pi; Rp=1; As=15;

N=ceil((log10((10^(Rp/10)-1)/(10^(As/10)-1)))/(2*log10(Wp/Ws))); Wc=Wp/((10^(Rp/10)-1)^(1/(2*N))); [b,a]=butter(N,Wc,'s'); [bz,az]=impinvar(b,a); w=[0:500]*pi/500; [H,w]=freqz(bz,az); subplot 221; plot(w/pi,abs(H));

xlabel('\\omega(\\pi)');

ylabel('|H(e^j^\\omega)|(dB)'); grid on; subplot 222;

plot(w/pi,20*log10(abs(H)));

xlabel('\\omega(\\pi)');ylabel('|H(e^j^\\omega)|(dB)' ); grid on;subplot 223;

plot(w/pi,angle(H)/pi); xlabel('\\omega(\\pi)'); ylabel('|H(e^j^\\omega)|'); grid on; subplot 224;

grd=grpdelay(bz,az,w);

plot(w/pi,grd);xlabel('\\omega(\\pi)'); ylabel('|H(e^j^\\omega)|'); grid on;

1)Bd()|j?0.5e|H(000.51?(?)10.5)|j?e(0|H-0.5-100.51?(?)

（2）切比雪夫I型 fs=10000;

0)Bd-50)|(j?e|H(-100-15000.51?(?)2015)|j?e10|H(5000.51?(?)

fp=1000; fst=1500; Rp=1; As=15; Ap=Rp;

epc=sqrt(10^(Ap/10)-1); Wp=2*pi*fp; Wc=Wp;

Wst=2*pi*fst;

N=ceil(acosh(sqrt(10^(0.1*As)-1)/epc)/acosh(Wst/Wc)); [b a]=cheby1(N,Rp,Wc,'s'); [bz az]=impinvar(b,a,fs); w=[0:500]*pi/500; [H w]=freqz(bz,az); subplot 221; plot(w/pi,abs(H)); xlabel('\\omega(\\pi)');

ylabel('|H(e^j^\\omega)|(dB)'); grid on; subplot 222;

plot(w/pi,20*log10(abs(H))); xlabel('\\omega(\\pi)');

ylabel('|H(e^j^\\omega)|(dB)'); grid on; subplot 223;

plot(w/pi,angle(H)/pi); xlabel('\\omega(\\pi)'); ylabel('|H(e^j^\\omega)|'); grid on; subplot 224;

grd=grpdelay(bz,az,w); plot(w/pi,grd);

xlabel('\\omega(\\pi)'); ylabel('|H(e^j^\\omega)|'); grid on;

1.550|H(ej?)|(dB)1|H(ej?)|(dB)00.5?(?)100.5-500-10000.5?(?)110.51510|H(ej?)|0-0.5-100.5?(?)1|H(e)|j?5000.5?(?)1

（3）切比雪夫II型 fs=10000; fp=1000; fst=1500; Rp=1; As=15; Ap=Rp;

epc=sqrt(10^(Ap/10)-1); Wp=2*pi*fp; Wst=2*pi*fst; Wc=Wst;

N=ceil(acosh(sqrt(10^(0.1*As)-1)/epc)/acosh(Wst/Wp)); [b a]=cheby2(N,As,Wc,'s'); [bz az]=impinvar(b,a,fs) w=[0:500]*pi/500; [H w]=freqz(bz,az); subplot 221; plot(w/pi,abs(H)); xlabel('\\omega(\\pi)');

ylabel('|H(e^j^\\omega)|(dB)'); grid on; subplot 222;

plot(w/pi,20*log10(abs(H)));

xlabel('\\omega(\\pi)');

ylabel('|H(e^j^\\omega)|(dB)' ); grid on;subplot 223;

plot(w/pi,angle(H)/pi); xlabel('\\omega(\\pi)'); ylabel('|H(e^j^\\omega)|'); grid on; subplot 224;

grd=grpdelay(bz,az,w); plot(w/pi,grd);

xlabel('\\omega(\\pi)'); ylabel('|H(e^j^\\omega)|'); grid on;

1.5)B1(d|)j?(eH|0.5000.51?(?)10.5|)j?(e0H|-0.5-100.51?(?)

5)B0(d|)(eH|-5-1000.51?(?)5|)j?(e0H|-500.51?(?)j?

（4）椭圆模拟原型滤波器 fs=10000; fp=1000; fst=1500; Rp=1; As=15; Ap=Rp;

epc=sqrt(10^(Ap/10)-1); Wp=2*pi*fp; Wst=2*pi*fst; Wc=Wp;

A=10^(As/20);

k1=epc/(sqrt(A^2-1)); k=Wp/Wst;

N=ceil(ellipke(k)*ellipke(sqrt(1-k1^2))/ellipke(k1)/ellipke(sqrt(1-k^2)));

[bz az]=ellip(N,Rp,As,fp/fs*pi) w=[0:500]*pi/500; [H w]=freqz(bz,az); subplot 221; plot(w/pi,abs(H)); xlabel('\\omega(\\pi)');

ylabel('|H(e^j^\\omega)|(dB)'); grid on; subplot 222;

plot(w/pi,20*log10(abs(H))); xlabel('\\omega(\\pi)');

ylabel('|H(e^j^\\omega)|(dB)'); grid on; subplot 223;

plot(w/pi,angle(H)/pi); xlabel('\\omega(\\pi)'); ylabel('|H(e^j^\\omega)|'); grid on; subplot 224;

grd=grpdelay(bz,az,w); plot(w/pi,grd);

xlabel('\\omega(\\pi)'); ylabel('|H(e^j^\\omega)|'); grid on;

1.5|H(ej?)|(dB)50|H(ej?)|(dB)100.5-50000.5?(?)1-10000.5?(?)110.5|H(ej?)|15100-0.5-100.5?(?)1|H(e)|j?5000.5?(?)1

用脉冲响应不变法设计IIR数字滤波器的优缺点及适用范围： 优点：运算简单，容易实现。可利用模拟滤波器设计技术。

缺点：由于发生混叠效应是阻带增益降低，容易不满足最初要求。且不能用FFT技术。不能够设计高通和带阻滤波器。

适用范围：不能够设计高通和带阻滤波器。用于设计规格化选频滤波器，h(n)无限长,用于模拟滤波器数字化。低阶次更容易实现。

五、实验心得与体会

通过这次实验，我掌握了利用MATLAB实现脉冲响应不变法设计IIR数字滤波器的具体方法。加深了对其设计原理的理解和认识，在实验过程中使我对数字滤波器和模拟滤波器之间的技术指标转化有了深刻理解，并且由直观的实验结果和理论课知识相结合使我掌握了脉冲响应不变法设计IIR数字滤波器的整个过程以及这种方法的优缺点及适用范围。所以，这次实验使我受益匪浅。

本文来源：https://www.bwwdw.com/article/eti6.html