数字信号处理实验教材正文

1、数字信号处理实验指导书实验一 离散时间信号与系统的实现与表示由于MATLAB数值计算的特点，用它来分析离散时间信号与系统是很方便的。一、离散时间信号在MATLAB中，可以用一个列向量来表示一个有限长度的序列，例如：在MATLAB中可以表示为：>>x = 2 1 0 2 3 -1 2 3;但是这种表示方法没有包含抽样的位置信息，要完全表示，要用和两个向量表示，例如：（箭头表示第0个抽样点的位置）在MATLAB中可以表示为：>> n=-4 -3 -2 -1 0 1 2 3;>> x = 2 1 0 2 3 -1 2 3;由于内存有限，MATLAB无法表示无限序列

2、。下面介绍几种典型离散信号的表示方法。1、单位脉冲序列（1）在MATLAB中zeros(1,N)函数来实现有限区间的，例如：>>x= zeros(1,N);>>x(1)=1;（2）我们还可以通过编写impseq函数来实现，该函数代码如下：function x,n= impseq(n0,n1,n2)%产生单位脉冲序列x(n)=(n-n0)%x,n= impseq(n0,n1,n2)% 1,n=n0% x =% 1,nn0%n=时间序列n=n1:n2;x=(n-n0)=0该函数产生如下信号，取值区间为，得到序列例1：输入命令：>>x,n= impseq(0,-5

3、,5)得到如下结果：x = 0 0 0 0 0 1 0 0 0 0 0n = -5 -4 -3 -2 -1 0 1 2 3 4 5 2、单位阶跃序列（1）在MATLAB中ones(1,N)函数来实现有限区间的，例如：>>x= ones(1,N);（2）我们还可以通过编写impseq函数来实现，该函数代码如下：function x,n= impseq(n0,n1,n2)%产生单位阶跃序列x(n)=u(n-n0);n1n0n2%x,n= impseq(n0,n1,n2)% 1,nn0% x =% 1,n<n0%n=时间序列n=n1:n2;x=(n-n0)>=0该函数产生如下

4、信号，取值区间为，得到序列例2：输入命令：>>x,n= impseq(0,-5,5)得到如下结果：x = 0 0 0 0 0 1 1 1 1 1 1n = -5 -4 -3 -2 -1 0 1 2 3 4 53、正弦序列sin函数就可以产生正弦波，例如：>>n= 0:N-1;>>x=A*sin(n)4、实指数序列MATLAB实现：>>n= 0:N-1;>>x=a.n5、复指数序列MATLAB实现：>>n= 0:N-1;>>x=exp(lu+j*w)*n);6、随机序列rand(1,N)：产生0,1上均匀分布的随

5、机序列；randn(1,N)：产生均值为0，方差为1的高斯随机序列，即白噪声序列。二、离散时间系统离散时间系统的阶差分方程为：在Matlab中，提供函数filter，用来求解给定输入和差分方程系数时，差分方程的数值解。调用格式为：Y=filter(B,A,X)其中：调用参数，分别为和各项差分项系数；为输入信号序列，注意必须保证系数不为零。输出序列与输入序列具有相同的长度。返回参数为差分方程的数值解例3：已知某线性时不变系统的差分方程如下：当输入序列为时的输出结果。 解：MATLAB程序如下： N=41;a=0.8 -0.44 0.36 0.22; %输入序列X各项差分项系数b=1 0.7 -0

6、.45 -0.6; %输出序列Y各项差分项系数x=1 zeros(1,N-1); %单位脉冲序列n=0:1:N-1; %时间序列y=filter(a,b,x); %输入单位脉冲序列得到单位脉冲响应stem(n,y);xlabel('n');ylabel('幅度');运行结果如下：三、实验内容1、编写程序用来产生下列基本序列，并绘制它们的波形图。（1）单位脉冲序列，起点为0，终点为10，在处有一单位脉冲；（2）单位阶跃序列，起点为0，终点为10，在前为0，后为1；（3）复指数序列，2、已知系统的差分方程为：求系统的单位脉冲响应并显示其波形图（选取）实验二 离散卷积

7、和一、卷积和的定义信号的卷积是针对时域信号处理的一种分析方法。信号的卷积一般用于求取信号通过某系统后的响应。在信号与系统中，我们通常求取某系统的单位冲激响应，所求得的可作为系统的时域表征。任意系统的系统响应可用卷积的方法求得：圆周卷积是数字信号处理的重要内容，由于它也是卷积的一种，必须了解其定义，圆周卷积的条件为卷积长度必须大于每个序列长度和。二、卷积工具箱函数序列的卷积和运算， MTLAB提供的计算两个离散序列卷积的函数conv，其调用方式为：y=conv(f1,f2)其中：调用参数为卷积运算所需的两个序列，返回值是卷积结果。MATLAB conv( )函数不需要给定的非零样值的时间序号，也

8、不返回卷积和序列的时间序号。即的返回参数中只有卷积的结果，没有的取值范围。此外，conv( )假定都是从开始，这就限制了它的应用范围（由离散序列卷积的性质可知，当序列和的起始点都为时，的取值范围为至length(x)+length(h)-2）。因此，要对从任意值开始的序列进行卷积和运算，同时正确标识出函数conv( )的计算结果各量，还须构造序列和的对应序号向量。例1：已知序列：，求其卷积和。解：MATLAB程序如下：n1=-2:2;f1=ones(1,length(k1);n2=0:3;f2=2.n2;f=conv(f1,f2)运行结果为：f = 1 3 7 15 15 14 12 8由于用

9、conv函数求出卷积后没有给出新序列所对应的时间变量。为此，我们可以构造一个实用点函数dconv( )，该函数实现的功能为：它在完成conv函数功能的同时，还产生一个对应新序列的时间变量，实现了序号向量的返回，并且可以绘制出序列和时域波形图。dconv函数MATLAB程序如下所示：function f,n=dconv(f1,f2,n1,n2)%the function of compute f=f1*f2%f: 卷积序列f(n)对应的非零样值向量%n: 序列f(n)的对因序号向量%f1: 序列f1(n)非零样值向量%f2: 序列f2(n)非零样值向量%n1: 序列f1(n)的对应序号向量%n2

10、: 序列f2(n)的对应序号向量f=conv(f1,f2) %计算序列f1与f2的卷积和fn0=n1(1)+n2(1); %计算序列f非零样值的起点位置n3=length(f1)+length(f2)-2; %计算卷积和f的非零样值的宽度n=n0:n0+n3 %确定卷积和f非零样值的序号向量 subplot(2,2,1);stem(n1,f1);title('f1(n)');xlabel('n') ;ylabel('f1(n)') %在子图1绘序列f1(k)时域波形图subplot(2,2,2); stem(n2,f2); title('

11、f2(n)'); xlabel('n'); ylabel('f2(n)'); %在子图2绘序列f2(k)时波形图subplot(2,2,3); stem(n,f); title('f1(n)与f2(n的卷积和f(n)'); xlabel('n'); ylabel('f(n)'); %在子图3绘序列f(k)的波形图例2：对于下面两个离散序列，我们可以调用dconv( )函数来求离散序列的卷积。解：MATLAB命令如下：f1=1 2 -1 -3 -2 4 ;n1=1:6;f2=ones(1,5);n2=-2:2

12、; %f1，f2两个信号的向量表示f,n=dconv(f1,f2,n1,n2) %子函数调用运行程序如下所示：f = 1 3 2 -1 -3 0 -2 -1 2 4n = -1 0 1 2 3 4 5 6 7 8f = 1 3 2 -1 -3 0 -2 -1 2 4n = -1 0 1 2 3 4 5 6 7 8需要注意的是，调用conv( )函数计算序列卷积时，该函数将向量和以外的序列样值均视为零，因此，当序列或为无限长序列时调用conv( )函数就可能出现误差。如果碰到无限长序列时候，我们必须将其截断才能求带入到conv( )函数中。此时，函数将把截断区域外的区间视为零，故conv( )计

13、算出的卷积只有部分是真实的。三、实验内容1、对下面三个序列，用conv函数验证卷积特性（交换律、结合律、分配律） 交换律 结合律 分配律 其中： 2、用dconv( )函数求例1的离散卷积和。实验三 离散系统的z域分析一、z域模型对于一个线性的离散系统，MALAB提供了6种表示方法：传递函数法，零极点增益法，状态空间法，部分分式展开法，二阶部分积和Lattice结构。每种方法都适用于不同的场合，考虑到理论知识的内容，在此尽介绍传递函数法，零极点增益法和带余数的部分分式展开法。1、传递函数法已知，则系统传递函数为：在MATLAB中用两个向量表示该系统。调用格式为：num=b1,b2,bm,bm+

14、1den=1,a1,an-1,an其中： num，den分别为分子多项式与分母多项式的系数，都是按的降幂进行排列的。2、零极点增益法MATLAB信号处理工具箱提供的tf2zp和zp2tf函数可进行系统函数的不同表示形式的转换。正幂有理多项式表示的系统函数为：用零点、极点和常数表示系统函数的零极增益模型为：（1）由传函模型得到零极增益模型MATLAB工具箱提供的tf2zp函数可以实现由系统的传函模型到零极增益模型的转换。调用格式为：z,p,k=tf2zp(b,a)其中：调用参数为，分别表示的分子和分母多项式的系数（从高到底排列）。返回参数为，分别表示零点、极点和增益通过tf2zp求出系统的零极点

15、后可根据极点的位置判断系统的稳定性。（2）由零极增益模型得到传函模型MATLAB工具箱提供的zp2tf 函数可以实现由系统的零极增益模型到传函模型的转换。调用格式为：b,a=zp2tf(z,p,k)其中：调用参数为，分别表示零点、极点和增益返回参数为，分别表示的分子和分母多项式的系数（从高到底排列）。例1：试将系统函数，表示为零极增益形式。解：用正幂有理多项式表示的系统函数为：MATLAB程序如下：%确定一阶因子形式的系统函数b=1 0 0.04 0;a=1 -0.8 0.16 -0.128;z,p,k=tf2zp(b,a);disp(零点);disp(z);disp(极点)；disp(p);

16、disp(常数)；disp(k);程序运行结果为：零点 0 0+0.2000i 0-0.2000i极点 0.8000 0.0000+0.4000i 0.0000-0.4000i常数 1由此可以写出系统函数的零极增益形式：（3）零极点图MATLAB工具箱提供了zplane函数用于绘制系统的零极点图。调用格式为：zplane(b,a)其中：调用参数为，分别表示的分子和分母多项式的系数（从高到底排列）。由zplane函数绘制得到的零极点图中，用符号表示零点；表示极点；符号旁的数字表示零极点的阶数；虚线画的圆表示单位圆。我们可以通过平面的零极点分布来判断系统的稳定性。例2：画出例1系统的零极点图，判别

17、其稳定性解：MATLAB程序如下：b=1 0 0.04 0;a=1 -0.8 0.16 -0.128;zplane(b,a)程序运行结果为：上由图可知，该系统极点全部位于单位圆内，故系统稳定。3、部分分式展开法（1）由一般式得到部分分式展开式若给出的一般式：可以由函数residuez得到的部分分式展开式。调用形式为：r,p,k=resifuez(b,a)其中：调用参数为，分别表示的分子和分母多项式的系数（从高到底排列）。返回参数为，其中为极点向量，为极点对应的系数向量，为余数多项式向量。（2）由部分分式展开式得到一般式若给出的部分分式展开式：可以由函数residuez得到的一般式。调用格式为：

18、b,a=residuez(r,p,k)其中：调用参数为，其中为极点向量，为极点对应的系数向量，为余数多项式向量。返回参数为，分别表示的分子和分母多项式的系数（从高到底排列）。所以上述部分分式展开式的residuez调用参数分别为：r=r1 r2 r3 r4p=p1 p2 p3 p3k=k1 k2同一极点在向量中出现了两次，表示是一个二阶的重极点。例3：试用MATLAB计算的部分分式展开式。解： MATLAB程序如下： b=1.5,0.98,-2.608,1.2,-0.144; a=1,-1.4,0.6,-0.072; r,p,k=residuez(b,a) %部分分式展开 disp(留书)；d

19、isp(r) disp(极点)；disp(p) disp(常数)；disp(k);程序运行结果为：留数0.7000+0.0000i 0.5000-0.0000i 0.3000极点0.6000+0.0000i 0.6000-0.0000i 0.2000常数0 2由此可以写出部分分式展开的结果： 二、z变换和逆z变换变换是时域离散信号和系统分析及设计的重要的数学工具。序列的变换为：对于给定的，使收敛的所有值的集合称为变换的收敛域，即满足：的反变换定义为：MATLAB的符号数学工具箱提供了计算正变换的函数ztrans和计算逆变换的函数iztrans。调用形式为：F=ztrans(f)该调用方式用于求

20、符号函数的变换，返回函数的自变量为；F=ztrans(f,w)该调用方式用于求符号函数的变换，返回函数的自变量为；F=ztrans(f,n,w)该调用方式用于对自变量为的符号函数求变换，返回函数的自变量为。f=iztrans(F)该调用方式用于对自变量为的符号函数求逆变换，返回函数的自变量为；f=iztrans(F,n)该调用方式用于对自变量为的符号函数求逆变换，返回函数的自变量为；f=iztrans(F,w,n)该调用方式用于对自变量为的符号函数求逆变换，返回函数的自变量为。例4：已知序列，求其变换。解：MATLAB程序如下：syms nf=sym('2(-n)'); %定义

21、序列F=ztrans(f) %求z变换程序运行结果为：F =2*z/(2*z-1)即：例5：已知一离散系统的系统函数，求其冲激响应。解：MATLAB程序如下：syms n zH=sym('z/(z2+3*z+2)');h=iztrans(H,n) %求逆z变换程序运行结果为：h =(-1)k-(-2)k，即：对象函数求逆变换，还可以利用residuez( )对象函数作部分分式展开，然后按部分分式展开法求得原函数。在此不做介绍。三、离散系统的频率响应分析若离散系统是稳定的，其系统函数的收敛域应包含单位圆，离散系统的频率响应即为单位圆上（）的系统函数，即：其中：为系统的幅频特性，为

22、系统的相频特性。在MATLAB中，利用freqz( )函数可方便的求出系统的频率响应。调用格式为：freqz(b,a)该调用方式将绘制系统在范围内的幅频特性和相频特性图；其中：分别为系统函数的分子、分母多项式的系数向量。freqz(b,a,whole)该调用方式将绘制系统在范围内的幅频特性和相频特性图。freqz(b,a,N)该调用方式将绘制系统在范围内个频率等分点的幅频特性和相频特性图，的缺省值为512。freqz(b,a,N,whole)该调用方式将绘制系统在范围内个频率等分点的幅频特性和相频特性图。此外，还有如下相类似的四种调用形式。H,w= freqz(b,a)H,w=freqz(b,

23、a,whole) H,w=freqz(b,a,N) H,w=freqz(b,a,N,whole) 其中，返回向量包含了离散系统频率响应在（或）范围内各频率点处的值，返回向量则包含了在（或）范围内个（的缺省值为512）频率等分点。利用这些调用方式MATLAB并不直接绘制系统的频率特性图，但可由向量用abs、angle、plot等函数来绘制幅频特性和相频特性图。例6：已知一离散系统的系统函数，试绘制频率特性图。解：MATLAB程序如下：b=0 1 0;a=1 0.3 0.2;freqz(b,a,'whole')程序运行结果为：四、实验内容1、已知一因果的LTI系统的函数为，求系统的

24、零极增益模型，部分分式模型，并用函数zplane画出系统的零极点分布，判断系统的稳定性。2、已知因果系统的系统函数为：求系统的单位脉冲响应并显示其波形图；若输入为，求系统的零状态响应并显示其波形图。实验四 离散傅立叶变换（DFT）DFT是数字信号处理中最重要的数学工具之一。其实质是对有限长序列频谱的离散化，即通过DFT使时域有限长序列与频域有限长序列相对应，从而可在频率域用计算机进行信号处理。更重要的是DFT有多种快速算法（FFT），可使信号处理速度提高好几倍，使数字信号的实时处理得以实现。因此，DFT即具有重要的理论意义又有广泛的实际应用价值。熟悉DFT的定义、物理意义和重要性质，有助于正确

25、使用DFT解决数字信号处理的实际问题。本次试验主要学习用MATLAB工具箱函数，以序列的时域和频域波形直观的验证DFT的物理意义及频域采样理论。一、DFT1、DFT定义设序列的长度为，则的点离散傅立叶变换对定义为：正变换： （1）逆变换： （2）其中，成为变换区间的长度。把（1）式改写为矩阵乘法运算：其中为序列行向量，是阶矩阵，通常称为旋转因子矩阵。上式用MATLAB的矩阵运算表示为：2、MATLAB程序点DFT的MATLAB程序如下：用矩阵乘法计算N点DFTclear; close allxn=input(请输入序列x=)N=length(xn);n=0:N-1; k=n; nk=n*k;

26、%生成0:N-1*0:N-1方阵WN=exp(-j*2*pi/N);Wnk=WN.nk; %产生旋转因子Xk=xn*Wnk; %计算N点DFT只要输入序列，运行该程序，即可实现的点DFT优点：概念清楚，编程简单缺点：占用内存大，运行速度低，不太实用二、FFT快速傅里叶变换简称FFT，它并不是与DFT不同的另一种变换，而是为了减少DFT计算次数的一种快速有效的算法。FFT算法基本可以分为两类：按时间抽取的FFT（DIT-FFT）算法和按频率抽取的FFT（DIF-FFT）。若，其中为整数，这时候的FFT称为基2FFT。关于FFT算法的具体原理在这里不作介绍，请大家自己查阅相关资料。MATLAB基础

27、部分提供了fft， ifft， fft2和ifft2等快速计算傅立叶变换的函数，它使DFT的运算速度量提高了若干的数量级。信号处理工具箱提供了等有关的变换函数。为了简化程序，后面的例题均直接调用这些函数。和用于计算一维快速正逆傅立叶变换，其调用格式如下：X=fft(x,N)用于计算序列向量的点DFT。默认时fft函数自动按的长度计算DFT。当为2的整数次幂时，fft按基2算法计算，否则用混合基算法。 调用格式相仿。例1：基本序列的离散傅立叶变换计算已知一下序列：复正弦序列：余弦序列：正弦序列：分别对和计算以上序列的点DFT，并绘出幅频特性曲线。解：直接产生序列调用函数求解本题的程序如下：cle

28、ar;close allN=16;N1=8;n=0:N-1;k=0:N1-1; %产生序列x1(n),计算DFTx1(n)x1n=exp(j*pi*n/8); %产生序列x1(n)X1k=fft(x1n,N); %计算N点DFTx1(n)Xk1=fft(x1n,N1); %计算N1点DFTx1(n)%产生序列x2(n),计算DFTx2(n)x2n=cos(pi*n/8); X2k=fft(x2n,N); %计算N点DFTx2(n)Xk2=fft(x2n,N1); %计算N1点DFTx2(n)%产生序列x3(n),计算DFTx3(n)x3n=sin(pi*n/8); X3k=fft(x3n,N)

29、; %计算N点DFTx3(n)Xk3=fft(x3n,N1); %计算N1点DFTx3(n)用，语句画出这三种信号、两种采样序列长度、共六种情况的曲线。图形划分及标注语句从略。运行结果如下：三、验证N点DFT物理意义长度为的序列的傅立叶变换和点离散傅立叶变换的定义式如下：式中，比较上面两式可以得出的傅立叶变换与的离散傅立叶变换的关系式为：简言之，序列点DFT的物理意义就是对在上进行点等间隔采样。为了验证这种关系，程序先对密集采样，绘制出幅频曲线图。然后再分别做8点和16点DFT来验证如上采样关系。为了能直观的看出和之间的采样关系，再的幅度和相位图中同时画出了的幅频和相频曲线。例2：（1），绘出

30、幅频曲线和相频曲线（2）计算并图示的8点DFT（3）计算并图示的16点DFT解：MATLAB程序如下：%计算并图示DFTx(n)clear;close allN1=8;N2=16; %两种FFT变换长度n=0:N1-1;w=2*pi*(0:2047)/2048Xw=(1-exp(-j*4*w)./(1-exp(-j*w) %对x(n)频谱函数采样2048点subplot(3,2,1);plot(w/pi,abs(Xw)title('x(n)的幅频曲线');xlabel('w/pi');ylabel('幅度');subplot(3,2,2);plo

31、t(w/pi,angle(Xw)title('x(n)的相频曲线');axis(0,2,-pi,pi);line(0,2,0,0)xlabel('w/pi');ylabel('相位(rad)');xn=(n>=0)&(n<4) %产生序列x(n)=xnX1k=fft(xn,N1); %计算N1点的DFTx(n)X2k=fft(xn,N2); %计算N2点的DFTx(n)figure(2)k1=0:N1-1subplot(3,2,1);stem(k1,abs(X1k),'.') %画N1点离散频谱幅度title(

32、'N1点DFTx(n)X1k')xlabel('k(w=2pik/N1)');ylabel('X1k幅度');hold onplot(N1/2*w/pi,abs(Ww) %叠加上连续频谱幅度subplot(3,2,2);stem(k1,angle(X1k),'.') %画N1点离散频谱相位title('X1k的相频曲线');axis(0,N1,-pi,pi);line(0,N1,0,0)xlabel('k(w=2pik/N1)');ylabel('相位(rad)');hold on

33、plot(N1/2*w/pi,angle(Xw) %叠加上连续频谱相位subplot(3,2,3);stem(k1,angle(X2k),'.') %画N2点离散频谱相位运行结果为：四、实验内容1、设是长度为的矩形序列，用matlab程序分析FFT取不同长度时的频谱变化2、设完成如下要求：（1）计算其DTFT，并画出区间的波形（2）计算4点DFT，并把结果显示在（1）所画的图形中（3）对补零，计算64点DFT，并显示结果。3、思考：用DFT理论解释例1中为何两种N值下的DFT结果差别如此之大？实验五 IIR数字滤波器的设计IIR数字滤波器设计的主要方法是先设计低通模拟滤波器，进

34、行频率变换，将其转换成相应的（高通的、带通的）模拟滤波器，再转换成为高通、带通或带阻数字滤波器。对设计的全过程的各个步骤，MATLAB都提供了相应的工具箱函数，使IIR数字滤波器设计变得非常简单。一、工具箱函数介绍 IIR数字滤波器的设计步骤可由下图的流程图来表示。滤波器系数模拟低通滤波器原型设计buttap，cheblap，cheb2ap，besselap，ellipap频率变换（变为高通、带通、带阻等）lp2lp，lp2hp，lp2bp，lp2bs求最小阶数buttord，cheblord，cheb2ord，ellipord合为一步的设计函数butter，cheb1，cheb2，ellip

35、，besself模拟数字变换bilinearimpinvar设计指标滤波器系数图1 IIR数字滤波器设计步骤流程图上图清楚的表示了5类20个信号处理工具箱函数的作用。下面以巴特沃斯滤波器设计函数为典型，介绍此流图中函数的功能和用法。其它类型的滤波器设计函数用法可以类推。1、求最小阶数N的函数buttord调用格式为：N,wc= buttord(wp,ws,Ap,As,s)其中：调用参数为滤波器指标；返回参数为巴特沃斯模拟滤波器的阶数及频率参数【注】及均在区间归一化，单位为弧度/秒；去掉最后的变元后，得到的就是数字滤波器，调用格式为： N,wc= buttord(wp,ws,Ap,As)其中：调

36、用参数为滤波器指标；返回参数为巴特沃斯数字滤波器的阶数及频率参数2、模拟低通滤波器原型设计函数buttap调用格式为：z,p,k= buttap(N)其中：调用参数为butterworth滤波器阶数；返回参数分别为滤波器的零点、极点和增益。得到后，很容易求出滤波器系数。3、模拟频率变换函数lp2lp调用格式为： Bt,At= lp2lp(B,A,wo)其中：调用参数为单位截止频率的模拟低通滤波器分母和分子多项式的系数；（单位：弧度/秒）为低通滤波器所期望的截止频率； 返回参数为得到的频率为的低通滤波器的分母和分子多项式的系数。4、模拟数字变换函数模拟数字变换函数包括双线性变换函数bilinea

37、r和脉冲响应不变法函数impinvar两种。（1）双线性变换函数bilinear MATLAB工具箱提供了函数bilinear，使用双线性变换法实现模拟滤波器的数字化。调用格式为：Bz,Az= bilinear (B,A,Fs)其中：调用参数为模拟滤波器传递函数分母和分子多项式的系数；为抽样频率； 返回参数为用双线性变换法设计得到的数字滤波器分母和分子多项式的系数。（2）脉冲响应不变法函数impinvarMATLAB工具箱提供了函数impinvar，使用脉冲响应不变法实现模拟滤波器的数字化。调用格式为： Bz,Az=impinvar(B,A,Fs)其中：调用参数为模拟滤波器传递函数分母和分子多

38、项式的系数；为抽样频率；其默认值为。 返回参数为用脉冲响应不变法设计得到的数字滤波器分母和分子多项式系数。5、把（2），（3），（4）合为一步的数字滤波器设计函数butter（N,wc,ftype）B,A=butter(N,Wc, ftype)其中：调用参数为阶数和截止频率； 返回参数设计得到的数字滤波器分母和分子多项式的系数；ftype用来控制设计得到的滤波器是低通、高通、带通还是带阻滤波器。B,A=butter（N,wc）设计低通或带通数字滤波器系数。 B,A=butter（N,wc,high）设计高通数字滤波器系数。 B,A=butter（N,wc,stop）设计带阻数字滤波器系数。【注

39、】（5）采用的是双线性变换函数bilinear，如果用脉冲响应不变法就得分开来作。由此可见：通常情况之下，有了（1）和（5）两类函数，数字滤波器设计问题就解决了。二、IIR滤波器设计1、低通巴特沃斯模拟滤波器设计低通巴特沃斯模拟滤波器的系统函数一般形式如下：极点：由此可见，低通巴特沃斯滤波器的系统函数完全由阶数和截止频率决定。而和是由滤波器设计指标决定的，其计算公式如下：，取：所设计的滤波器的通带指标刚好满足，阻带指标富裕。取：所设计的滤波器的阻带指标刚好满足，通带指标富裕。MATLAB工具箱buttord，butter就是根据以上关系式编写的，因此设计的时候不需要再记忆和使用这些公式了。例1

40、：设计一个低通巴特沃斯模拟滤波器指标如下：通带截止频率：，通带最大衰减：通带截止频率：，阻带最小衰减：解：MATLAB程序如下：%低通巴特沃斯低通滤波器设计clear;close allfp=3400;fs=4000;Ap=3;As=40; %输入滤波器指标N,fc=buttord(fp,fs,Ap,As,'s') %计算阶数N和3dB截止频率fcB,A=butter(N,fc,'s'); %设计低通巴特沃斯模拟滤波器hf,f=freqs(B,A,1024); %计算模拟滤波器频率响应，freqs为工具函数plot(f,20*log10(abs(hf)/abs(

41、hf(1)grid;xlabel('f/Hz');ylabel('幅度(dB)');axis(0,4000,-40,5)line(0,4000,-3,-3);line(0,4000,-3,-3)line(3400,3400,-90,5)2、模拟低通滤波器转换成数字低通滤波器模拟滤波器离散化的基本方法有脉冲响应不变法和双线性变换法。1、脉冲响应不变法函数impinvar 对只有单阶极点的阶模拟滤波器，用部分分式展开为：则脉冲响应不变法取：MATLAB工具箱函数impinvar用来实现这个转换。2、双线性变换函数bilinear双线性变换原理是用替换中的值，得到。M

42、ATLAB工具箱函数bilinear用来实现这个转换。例2：已知一模拟滤波器的系统函数为，分别用脉冲响应不变法和双线性变换法将转换成数字滤波器系统函数，并图示和的幅频响应曲线。分别取采样频率和观察脉冲响应不变法中存在的频率混叠失真和双线性变换法中存在的非线性频率失真。解：MATLAB程序如下：%用脉冲响应不变法和双线性变换法将模拟滤波器离散化clear;close allb=1000;a=1,1000;w=0:1000*2*pi; %设定模拟频率hf,w=freqs(b,a,w); %计算模拟滤波器的频响函数subplot(2,3,1) %画出模拟滤波器的幅频特性plot(w/2/pi,abs

43、(hf);grid;title('(a)')xlabel('f(Hz)');ylabel('幅度');Fs0=1000,500;for m=1:2 Fs=Fs0(m) %T=0.001s及T=0.002s d,c=impinvar(b,a,Fs) %用impinvar函数实现离散化 wd=0:512*pi/512 %设定数字归一化频率 hw1=freqz(d,c,wd) % 计算数字滤波器频响函数subplot(2,3,2);plot(wd/pi,abs(hw1)/abs(hw1(1);hold on; %画出数字滤波器幅频特性endgrid;x

44、label('f/Hz');text(0.25,0.88,'T=0.002s');text(0.12,0.54,'T=0.001s');for m=1:2 Fs=Fs0(m) %T=0.001s及T=0.002s f,e=bilinear(b,a,Fs) %用bilinear函数实现离散化 wd=0:512*pi/512 %设定数字归一化频率 hw2=freqz(f,e,wd) % 计算数字滤波器频响函数subplot(2,3,3);plot(wd/pi,abs(hw2)/abs(hw2(1);hold on; %画出数字滤波器幅频特性endgr

45、id;xlabel('f/Hz');text(0.5,0.74,'T=0.002s');text(0.12,0.34,'T=0.001s');三、实验内容1、运行例题一及例题二程序；显示运行结果；根据运行结果说明脉冲响应不变法和双线性变换法各自的特点。2、用双线性法设计一个Butterworth低通数字滤波器，频率在处衰减，在处衰减，采样频率为，确定系统函数。并图示和的幅频响应曲线。【提示】数字低通双线性变换与脉冲不变法程序：wp=500*2*pi;ws=150*2*pi;Ap=3;As=15;Fs=2000;n,Wn=buttord(wp,ws

46、,Ap,As,'s');b,a=butter(n,Wn,'s');bn1,an1=bilinear(b,a,2000); %双线性变换H1,w=freqz(bn1,an1);%采用脉冲响应不变法进行离散化处理bn2,an2=impinvar(b,a,2000); %脉冲响应不变法H2,w=freqz(bn2,an2)实验六 FIR数字滤波器的设计 IIR数字滤波器的设计简单方便，特别是采用双线性变换法设计的数字滤波器不存在频谱混叠问题，效果较好。但是IIR数字滤波器有一个明显的缺点，就是其相位特性一般是非线性的。如果滤波器在有效传输频带内的相位特性不是线性的，将

47、造成有用信号的传输失真。在有些实际应用场合，例如数据传输，图像处理等，对滤波器的线性相位特性要求颇为严格，所以在这些场合中IIR滤波器一般是不能胜任的。FIR滤波器最大的优点是容易设计成线性相位特性，其幅度特性可以随意设计，而且不存在稳定性问题。设FIR滤波器的单位脉冲响应为，则其系统函数为：上式表明，FIR滤波器的系统函数为的多项式，而IIR滤波器的系统函数为的有理分式形式。因此，FIR滤波器在S平面上找不到与之响应的模拟系统函数，也就是说，FIR滤波器的设计不能借用模拟滤波器设计的一套成熟理论。MATLAB信号处理工具箱提供的FIR数字滤波器的设计方法有以下两种：窗函数法，等波纹最佳一致逼

48、近法。本次实验主要介绍窗函数法。一、MATLAB提供的窗函数窗函数法设计滤波器时，阶数选择如下表所示：表1 六种窗函数特性的表窗函数旁瓣峰值近似过渡带精确过渡带最小阻带衰减矩形窗巴特利特汉宁海明布拉克曼凯泽窗Matlab提供了相应的子程序来实现本节中的窗函数，如下：wd=boxcar(N) %数组中返回点矩形窗函数wd=triang(N) %数组中返回点Bartlett（三角）窗函数wd=hanning(N) %数组中返回点汉宁窗函数wd=hamming(N) %数组中返回点海明窗函数wd=blackman(N) %数组中返回点布拉克曼窗函数wd=kaiser(N,beta) %数组中返回给定

49、beta值时点凯泽窗函数【注】以上函数的输入变元一般只要窗函数的长度就够了，只有凯泽窗还需要规定值。输出变元就是中心值归一化为1的窗函数序列，它是列向量。利用语句plot就可以画出窗函数序列的形状。二、窗函数法设计FIR DF的步骤1、设计步骤（1）根据对过渡带宽及阻带衰减要求，选择窗函数的类型并估计窗口长度（或阶数）窗函数类型可根据其阻带最小衰减的条件独立选择，因为其长度对阻带最小衰减没有影响。在确定窗函数之后，就可根据过渡带宽小于给定指标的条件，确定所拟用的窗函数长度。（2）根据和求出相应的理想滤波器单位脉冲响应。在一般情况下，是不能用封闭公式表示的，需要采用数值方法。从到采样点，用IDF

50、T来求出。（3）计算FIR数字滤波器单位脉冲响应。由MATLAB得出的是一个无限长序列，因而不可用FIR滤波器实现。所以要选择合适的窗函数来截取的适当长度（即阶数），以保证实现要求的阻带衰减；在matlab中用点乘命令表示为：（4）验证技术指标是否满足要求。为了计算数字滤波器在频域中的特性，可以用freqz子程序。例1：分别用矩形窗和Hamming窗设计线性相位的FIR数字低通滤波器。要求通带截止频率，单位脉冲响应的长度。绘出及其幅频响应特性曲线。解：【建模】 用窗函数法设计FIR滤波器时，先求出相应的理想低通滤波器单位脉冲响应，在根据阻带最小衰减选择合适的窗函数，最后得到FIR数字滤波器单位脉冲响应。本题中，所以线性相位理想低通滤波器的单位脉冲响应为：为了满足线性相位FIR滤波器条件，要求【程序】%窗函数法设计FIR低通滤波器clear;close allN=21;wc=pi/4n=0:N-1;r=(N-1)/2;hdn=sin(wc*(n-r)/pi./(n-r); %计算理想低通单位脉冲响应 if rem(N,2