数字信号处理 实验 冲激响应不变法IIR数字滤波器设计.doc

数字信号处理实验报告
实验名称：冲激响应不变法IIR数字滤波器设计
学号： 姓名：
评语：
成绩：
一、 实验目的
1、掌握构成一个频率响应与给定的滤波特性相接近的模拟滤波器的设计原理。
2、掌握用冲激响应不变法设计IIR数字滤波器的基本原理和算法。
3、了解数字滤波器和模拟滤波器的频率响应特性，掌握相应的计算方法，分析用冲激响应不变法获得的数字滤波器频率响应特性中出现的混叠现象。
二、实验原理与计算方法
1、冲激响应不变法设计IIR数字滤波器的基本原理和算法
采用冲激响应不变法设计数字滤波器，就是使其单位样值响应与相应的模拟滤波器的冲激响应在抽样点处的量值相等，即
　　　　　　 (1)
其中T为抽样周期。
因此用冲激响应不变法设计IIR数字滤波器的基本步骤，就是首先根据设计要求确定相应的模拟滤波器的传递函数，经Laplace反变换求出冲激响应后，对它进行抽样得到的等于数字滤波器的单位样值响应，再经z变换所得就是数字滤波器的传递函数。
如果模拟滤波器的传递函数的N个极点都是单极点，则可以将写成部分分式展开的形式
　　　　　　　　　　　　　　　 (2)
那么，经Laplace反变换求出的模拟滤波器的冲激响应为

相对应的数字滤波器的单位样值响应为

对上式作z变换，得
(3)
由上面的推导可见，只要模拟滤波器的传递函数的N个极点都是单极点，当已经求出各个极点值和部分分式的系数后，则可以从模拟滤波器的传递函数的表达式（2）直接得到数字滤波器的传递函数的表达式（3）。
2、Butterworth和Chebyshev模拟滤波器的设计方法
（1）Butterworh原型
MATLAB提供了函数[z,p,k]=buttap(N)用来设计N阶归一化的Butterworth模拟低通滤波器；函数[b,a]=u_buttap(N,Omegac)给出未归一化的Butterworth模拟低通滤波器原型；再利用函数[C,B,A]=sdir2cas(b,a)得到级联形式的N阶Butterworth模拟低通滤波器原型。>>N=3;OmegaC=0.5;
>> [b,a]=u_buttap(N,Omegac);
>> [C,B,A]=sdir2cas(b,a)
（2）按给定技术指标设计Butterworth模拟低通滤波器
函数[b,a]=afd_butt(Wp,Ws,Rp,As)用来实现按给定技术指标设计Butterworth模拟低通滤波器；函数[db,mag,pha,w]=freqs_m(b,a,0.5*pi)绘制模拟滤波器频响特性，得出衰减值,幅频特性,相频特性和自变量w；函数[H,w]=freqs(b,a,0.5*pi)绘制模拟滤波器频响特性,给出复函数H,再由abs(H)和angle(H)求出幅频和相频特性；[ha,x,t]=impulse(b,a)求出模拟滤波器冲激响应。
>>Wp=0.2*pi;Ws=0.3*pi;Rp=7;As=16;
>> [b,a]=afd_butt(Wp,Ws,Rp,As)
>> [db,mag,pha,w]=freqs_m(b,a,0.5*pi)
>> [ha,x,t]=impulse(b,a)
（3）数字滤波器设计－Butterworth
函数[b,a]=imp_incr(c,d,T)用来实现冲激响应不变的映射，c、d分别为模拟滤波器系统函数Ha(s)的分子、分母的系数；T为抽样周期；
%数字滤波器技术指标
>>wp=0.2*pi;ws=0.3*pi;Rp=1;As=15;
%对应的模拟滤波器技术指标
>>T=1;Wp=wp/T;Ws=ws/T;
>>[cs,cd]=afd_butt(Wp,Ws,Rp,As);
>> [b,a]=imp_invr(cs,ds,T)
%并联形式的系统函数系数
>>[C,B,A]=dir2par(b,a)
函数[H,w]=freqz(b,a)绘制数字滤波器频响特性,给出复函数H,再由abs(H)和angle(H)求出幅频和相频特性。
（4）Chebyshev原型
MATLAB提供了函数[z,p,k]=cheblap(N,Rp)用来设计N阶通带波动为Rp的归一化Chebyshev模拟低通滤波器；函数[b,a]=u_chb1ap(N,Omegac)给出未归一化的Chebyshev模拟低通滤波器原型。
（5）按给定技术指标设计Chebysheb模拟低通滤波器
函数[b,a]=afd_chb1(Wp,Ws,Rp,As)用来实现按给定技术指标设计Chebyshev模拟低通滤波器；函数[db,mag,pha,w]=freqs_m(b,a,0.5*pi)绘制频响特性；[ha,x,t]=impulse(b,a)求出模拟滤波器冲激响应。
>>Wp=0.2*pi;Ws=0.3*pi;Rp=1;As=16;
>> [b,a]=afd_chb1(Wp,Ws,Rp,As)
>> [db,mag,pha,w]=freqs_m(b,a,0.5*pi)
>> [ha,x,t]=impulse(b,a)
（6）数字滤波器设计－Chebyshev
函数[b,a]=imp_incr(c,d,T)用来实现冲激响应不变的映射，c、d分别为模拟滤波器系统函数Ha(s)的分子、分母的系数；T为抽样周期；
%数字滤波器技术指标
>>wp=0.2*pi;ws=0.3*pi;Rp=1;As=15;
%对应的模拟滤波器技术指标
>>T=1;Wp=wp/T;Ws=ws/T;
>>[cs,cd]=afd_chb1(Wp,Ws,Rp,As);
>> [b,a]=imp_invr(cs,ds,T)
%并联形式的系统函数系数
>>[C,B,A]=dir2par(b,a)
3、数字滤波器和模拟滤波器的频率响应特性
当已经求出模拟滤波器的传递函数和数字滤波器的传递函数后，模拟滤波器的频率响应为
　 　　　　　　　　　　　 　　　 (4)
而数字滤波器的频率响应为
(5)
根据连续时间信号的Laplace变换和抽样后所得的序列的z变换之间的关系,有

将(4)和(5)式代入上式，得模拟滤波器和数字滤波器的频率响应之间的关系为
　　　 　 (6)
上式表明，数字滤波器的频谱是模拟滤波器频谱的周期延拓，其周期为抽样角频率。如果模拟滤波器的频谱带宽不是有限的，数字滤波器的频谱在高频区（二分之一抽样角频率附近）就会造成频谱的混叠，如果频率高时幅度频谱衰减较小，这种混叠现象将相当明显。
由(6)式可以看出，当抽样周期T很小时，数字滤波器的幅度频谱有极高的增益，为避免这样的问题，可将(3)式改写为
　　 　　　　　　　　　　　　　(7)

实验内容及结果
（1）设计一个Butterworth滤波器，满足如下级数指标：
通带边界频率wp=0.4π，通带衰减函数ap=0.5dB，
阻带边界频率ws=0.6π，阻带衰减函数as=50dB，
采用冲激响应不变法，选取T=2，记录所得的模拟滤波器的阶数N，求出有理函数形式的系统函数，画出模拟滤波器和数字滤波器的频率响应的幅频和相频特性曲线以及单位抽样响应h(n)和冲激响应ha(t)的图像。
a、 模拟Butterworth滤波器
实验代码：
Wp=0.4.*pi;
Ws=0.6.*pi;
Rp=0.5;
As=50; %对应的模拟滤波器技术指标
[b,a]=afd_butt(Wp,Ws,Rp,As); %实现按给定技术指标设计
[db,mag,pha,w]=freqs_m(b,a,0.5*pi); %绘制模拟滤波器频响特性
[ha,x,t]=impulse(b,a); %求出模拟滤波器冲激响应
[H,w]=freqs(b,a); %绘制模拟滤波器频响特性
subplot(3,1,1);
plot(abs(H));
title('模拟滤波器的幅频特性')
subplot(3,1,2);
plot(angle(H));
title('模拟滤波器的相频特性')
subplot(3,1,3);
plot(ha);
title('模拟滤波器冲激响应')
实验截图：

b、数字Butterworth滤波器
实验代码：
wp=0.4*pi;
ws=0.6*pi;
Rp=0.5;
As=50; %对应的模拟滤波器技术指标
T=2;
Wp=wp/T;
Ws=ws/T;
[cs,cd]=afd_butt(Wp,Ws,Rp,As);
[b,a]=imp_invr(cs,cd,T); %并联形式的系统函数系数
[C,B,A]=dir2par(b,a);
[H,w]=freqz(b,a) ; %函数[H,w]=freqz(b,a)绘制数字滤波器频响特性,
%给出复函数H,再由abs(H)和angle(H)求出幅频和相频特性。[hn,t]=impz(b,a);
subplot(3,1,1);
plot(abs(H));
title('数字Butterworth幅频特性')
subplot(3,1,2);
plot(angle(H));
title('数字Butterworth相频特性')
subplot(3,1,3);
plot(hn);
title(' hn抽样响应')
实验截图：

（2）将（1）中所得的模拟滤波器改为Chebyshev逼近，记录滤波器的阶数N，描绘出模拟滤波器和数字滤波器的幅频和相频特性曲线以及冲激响应ha(t)和单位抽样响应h(n)的图像。
a、模拟Chebyshev滤波器
实验代码：
Wp=0.4.*pi;
Ws=0.6.*pi;
Rp=0.5;
As=50; %对应的模拟滤波器技术指标
[b,a]=afd_chb1(Wp,Ws,Rp,As); %实现按给定技术指标设计
[db,mag,pha,w]=freqs_m(b,a,0.5*pi); %绘制模拟滤波器频响特性
[ha,x,t]=impulse(b,a); %求出模拟滤波器冲激响应
[H,w]=freqs(b,a); %绘制模拟滤波器频响特性
subplot(3,1,1);
plot(abs(H));
title('模拟滤波器的幅频特性')
subplot(3,1,2);
plot(angle(H));
title('模拟滤波器的相频特性')
subplot(3,1,3);
plot(ha);
title('模拟滤波器冲激响应')
实验截图：

b、数字Chebyshev滤波器
实验代码：
wp=0.4*pi;
ws=0.6*pi;
Rp=0.5;
As=50; %对应的模拟滤波器技术指标
T=2;
Wp=wp/T;
Ws=ws/T;
[cs,cd]=afd_chb1(Wp,Ws,Rp,As);
[b,a]=imp_invr(cs,cd,T); %并联形式的系统函数系数
[C,B,A]=dir2par(b,a);
[H,w]=freqz(b,a); %函数[H,w]=freqz(b,a)绘制数字滤波器频响特性,
%给出复函数H,再由abs(H)和angle(H)求出幅频和相频特性。
[hn,t]=impz(b,a);
subplot(3,1,1);
plot(abs(H));
title('数字Butterworth幅频特性');
subplot(3,1,2);
plot(angle(H));
title('数字Butterworth相频特性');
subplot(3,1,3);
plot(hn);
title('hn抽样响应')
实验截图：