数字信号处理实验讲义

1、数字信号处理讲义 数字信号处理实验指导书   曲阜师范大学物理工程学院信号与信息处理教研室2012.0369 前 言数字信号处理是是电子信息工程专业的一门专业基础课。本课程主要研究如何对信号进行分析、变换、综合、估计与识别等加工处理的基本理论和方法，数字信号处理实验是验证、巩固和补充课堂讲授的理论知识的必要环节。通过实验，使学生巩固所学基本理论，掌握最基本的数字信号处理的理论和方法，提高综合运用所学知识，提高计算机编程的能力。进一步加强学生独立分析问题、解决问题的能力、综合设计及创新能力的培养，同时注意培养学生实事求是、严肃认真的科学作风和良好的实验习惯

2、，为今后的工作打下良好的基础。数字信号处理实验指导书针对每个实验介绍了MATLAB语言数字信号处理工具箱中的相应函数，举例并附有相应的程序。为配合课堂理论教学，实验内容安排仍从认识性和验证性入手，逐步增加设计性和工程应用性内容，达到训练实验技能和积累工程实际应用经验之目的。数字信号处理实验成绩占本门课程总成绩的20%。实 验 要 求在实验过程中，要求学生做到：1、预习实验指导书有关部分，认真做好实验内容的准备工作，就实验可能出现的情况提前作出思考和分析，需要计算的参数提前完成计算工作，并认真写出预习报告。2、仔细观察实验过程中图形随参数的变化，记录图形变化的主要情况，作出必要说明和分析。3、认

3、真书写实验报告并在规定的时间内把实验报告交给辅导教师。实验报告包括实验目的和要求，实验情况及其分析。对需要编程的实验，写出程序设计说明，给出源程序框图和清单以及实验结论、相关图表，并做必要的分析总结。4、遵守机房纪律，服从辅导教师指挥，爱护实验设备，注意卫生。5、实验课程不迟到。如有事不能出席，所缺实验一般不补。实验验收分为两个部分。第一部分是上机操作，包括检查程序运行和即时提问。第二部分是提交书面的实验报告。每个实验都应当在规定的时间内完成并检查通过，过期视为未完成该实验，扣该实验操作成绩。为避免期末集中检查方式产生的诸多不良问题，希望同学们抓紧时间，合理安排，认真完成。目 录实验一 MAT

4、LAB基本命令1实验二 离散时间信号分析4实验三 离散系统分析17实验四 离散傅立叶变换30实验五 快速傅里叶变换（FFT）及其应用37实验六 IIR滤波器设计42实验七 FIR滤波器设计55实验八 窗函数及其在谱分析中的作用59实验一 MATLAB基本命令一、实验目的：掌握Matlab基本命令的使用，以便在以后的实验和进一步的学习中使用Matlab命令。实验任务：1、 熟悉Matlab 语言环境，学习基本命令的使用，并要求掌握。2、 熟悉Matlab中的信号处理Demo使用，并了解其中的程序结构和命令。MATLAB代表Matrix Laboratory。 MATLAB 最早是为方便矩阵存取而

5、编写的一种软件，现已发展成为一个具有高性能数值计算和可视化功能的科学计算环境。MATLAB集成了数值分析、矩阵计算、信号处理和图形等众多功能，具有简单易用的环境，问题的提出和解答只须以数学方式表达和描述，不需要大量原始而传统的编程过程。 MATLAB 一些简单命令：1. Demo命令：MATLAB例程演示初次学习MATLAB可以使用demo命令来熟悉MATLAB语言环境。你可以在MATLAB的窗口中直接敲demo，它将带你学习使用MATLAB。2. Help命令：在线帮助。语法格式为Help topic topic 为要取得的帮助主题。例如：help fft 将显示fft( )函数的说明。3.

6、 编写MATLAB程序：MATLAB 的程序以 .M 结尾，你可以使用MATLAB窗口中File-àNewàM-file 的菜单命令来建立一个.m空白文件，编写完成后，你保存一下，然后可以在MATLAB窗口中敲此文件名执行该文件。 函数编写 script编写4. 常用的MATLAB命令归类：通用命令：Help 在线帮助Demo 例程演示Path 路径设置Who 列出当前变量What 列出当前路径下的.m文件Clear 清除所有变量Clc 清除当前命令窗口命令Clf 清除当前图形Quit 退出MATLAB运算符：加法： +减法： 乘法： *点乘： .*指数： . 除法： 点除

7、： .关系运算符：小于 < 大于 > 小于等于 <= 大于等于 >= 等于 = 不等于 =基本矩阵zeros 零阵 ones 1阵 信号运算abs取模值angle相角real实部imag虚部conv卷积fft快速傅立叶变换ifft 快速傅立叶反变换图形函数figure生成一图框axis设置坐标轴hold保持原图text在图上标记文字gtext通过鼠标点击方式在图上标记文字plot画图subplot 子图title添加图名xlabel给X轴增加文本标注ylabel给Y轴增加文本标注grid设置网格线stem离散序列数据杆状图bar条形图polar极坐标实验二 离散时间信号

8、分析实验目的：通过Matlab实现基本信号的表示和运算。理解信号的构成和性质。实验任务：1 学习使用基本信号在Matlab 中的基本表示方法。2 了解本实验中的扩展函数的结构和组成，实现自己设计的扩展函数的编程： （1）写一新的卷积函数conv_m，它可求出带下标的序列卷积。（2）修改给出的evenodd函数，使他能接受任意序列并把它分解成共轭对称分量和共轭反对称分量。3 掌握复杂信号的表示方法；并对复杂信号进行编程实现。 （1）产生下列序列并绘出离散图：（a） x(n)=2(n+2)-(n-4) -5n5（b） x(n)=nu(n)-u(n-10)+10e-0.3(n-10)u(n-10)-

9、u(n-2) 0n20（c） x(n)=cos (0.04n)+0.2w(n) 0n50 其中w(n)是均值为零、方差为1的白噪声序列（d） （n）=x1 x1 x1 x1,其中x1=5 4 3 2 1 (e) 已知x(n)=u(n)-u(n-10),要求将它分解成奇偶序列。(2) 用Matlab产生一矩形脉冲串，连续时间为2s,脉冲个数为10，并绘制其波形图。4 掌握复序列x(n)的实部、虚部、幅值和相位离散图，并编程实现。 产生复信号： x(n)=e(-0.1+j0.3)n -10n10并画出复序列x(n)的实部、虚部、幅值和相位离散图。一离散信号在数字信号处理（DSP）中，所有的信号都是

10、离散（时间）信号，因此首先应解决在MATLAB中如何表示离散信号。设一模拟信号经A/D变换后，得到序列信号x(n)=x(n)= , x (-1), x ( 0 ) , x (1 ),由于MATLAB对下标的约定为从1开始递增，因此要表示x(n),一般应采用两个矢量，如:n=-3,-2,-1,0,1,2,3,4,5y=1,-1,3,2,0,4,5,2,1这表示了一个含9个采样点的矢量：y(n)=x(-3), x (-2) , x (-1) , x (0) ,x (1 ) , x (5) 。通常情况下，序列值从x(0)开始，因此一个N点序列x(n)=x (0) ,x (1) , ,x (N-1)可

11、简单地表示成y=x (0), x (1) , , x (N-1)这时y的下标为1N。二基本信号表示 1单位取样序列这一函数可利用MATLAB的zeros函数实现：x=zeros(1,N);x(1)=1 还可以借助于关系操作符实现：n=1：N;x=n=1;如要产生 则可采用MATLAB实现：n=n1：n2;x=(n-n0)=0;2单位阶跃序列这一函数可利用MATLAB的ones函数实现：x=ones(1,N);还借助于关系操作符“>=”来实现。如要产生在n1n0n2上的单位阶跃序列则可采用MATLAB实现：n=n1:n2;x=(n-n0)>=0;3实指数序列x(n)=an 采用MAT

12、LAB实现：n=0:N-1;x=a.n;4.复指数序列采用MATBLAB实现：n=0:N-1;x=exp(lu+j*w0)*n);5正（余）弦序列采用MATLAB实现：n=0:N-1;x=cos(w0*n+Q);6随机序列MATLAB中提供了两类（伪）随机信号：rand(1,N)产生0，1上均匀分布的随机矢量；randn(1,N)产生均值为0，方差为1的高斯随机序列，也就是白噪声序列。其它分布的随机数可通过上述随机数的变换而产生。7周期序列x(n)=x(n+N) 例如，设x1表示x序列中一个周期的序列，要产生P个周期的x序列，可以把它重复4次，用MALAB实现：x=x1,x1,x1,.,x1;

13、但是高明的方法是利用MATLAB的强有力的下标能力。先产生一个包含P行x(n)值的矩阵，然后用结构（：）把它的P行串接起来成为一个长行，不过这种结构只能用于列向，所以我们往往还要用矩阵转置把它扩展成行向。>>xtilde=x1*ones(1,P); %P列x1>>xtilde=xtilde（：）； %长的列向量>>xtilde=xtilde; 长的行向量三 序列操作1信号加x(n)=x1(n)+x2(n)采用MATLAB实现：x=x1+x2;注意：当x1和x2序列的长度不等或其位置不对应时，信号相加就不是这么简单。首先应使x1、x2具有相等的长度，然后两者对

14、齐，最后进行相加。2信号乘x(n)=x1(n)*x2(n)这是一种样本对样本的相乘，也即点乘运算，在MATLAB中可采用.*（数组乘法）来实现，但两序列x1、x2也应经过处理。3改变比例y(n)=ax(n)=ax(n)采用MATLAB实现：y=alpha*x;4移位y(n)=x(n-k)5折叠y(n)=x(-n)它将序列x(n)在n=0处倒转，在MATLAB中可直接用fliplr函数实现。6取样和它不同于信号加。采用MATLAB实现：y=sum(x(n1:n2);7取样积采用MATLAB实现：y=prod(x(n1:n2);8卷积和 y(n)=x(n)*h(n)=可直接采用MATLAB中的函数

15、conv，即y=conv(x,h)；它默认序列从n=0开始。但如果序列是从一负值开始，即如x(n):nxbnnxeh(n):nhbnnhe其中nxb<0或nhb<0，或两者同时为负时，就不能直接采用conv函数。通过卷积结果序列 y(n):nybnnye且 nyb=nxb+nhb , nye=nxe+nhe这样我们就可以构成一新的卷积函数conv_m，它可求出带下标的序列卷积。9信号能量采用MATLAB实现：Px=sum(abs(x).2); Px=sum(x.*conj(x);10奇偶综合任何一个序列x(n)都可分解成偶对称部分xe(n)和奇对称部分x0(n)，即x(n)=xe(

16、n)+x0(n)其中xe(n)=x(n)+x(-n)x0(n)=x(n)-x(-n)这样我们可以设计一函数evenodd，完成将任一给定序列x(n) 分解成xe(n)和 x0(n)。四MATLAB信号生成函数1square（）方波发生器。调用格式为：square(T) 返回一周期为2的方波，采样时刻由向量T指定。square(T,DUTY) 产生一给定占空比的方波，占空比DUTY是信号为正值的比例。2 sawtooth（）锯齿波和方波发生器。调用格式为：sawtooth(T) 返回一周期为2的锯齿波，采样时刻由向量T指定。sawtooth(T,WIDTH) 产生一三角波，WIDTH指定最大值出

17、现的地方，其取值在0到1之间，1对应2。3 diric（）周期sinc函数（Dirichlet函数）。调用格式：Y=diric(X,N) 返回一大小与X相同的矩阵，其元素为Dirichlet函数。N必须为正整数，该函数将0到2等间隔分成N等分。4 pulstran （）通过对连续函数或脉冲原型进行采样而得到脉冲序列的发生器5 rectplus（）非周期矩形波发生器6 tripuls（）非周期三角脉冲发生器五MATLAB七个基本的扩展函数1单位取样序列impseq.mfunction x,n=impseq(n0,n1,n2)%Generates x(n)=delta(n-n0);n1<=n

18、,n0<-n2%-%x,n=impseq(n0,n1,n2)%if(n0<n1)|(n0>n2)|(n1>n2) error(arguments must satisfy n1<=n0<=n2)end n=n1:n2;x=(n-n0)=0;2单位阶跃序列stepseq.mfunctionx,n=stepseq(n0,n1,n2)%Generates x(n)=u(n0);n1<=n0<=n2%-%x,n=stepseq(n0,n1,n2)%if(n0<n1)|(n0>n2)|(n1>n2) error(arguments mus

19、t satisfy n1<=n0<=n2)end n=n1:n2;x=(n-n0)>=0;3信号加sigadd.mfunctiony,n=sigadd(x1,n1,x2,n2)%implements y(n)=x1(n)+x2(n)%-%y,n=sigadd(x1,n1,x2,n2)%y=sum sequence over n, which includes n1 and n2%x1=first sequence over n1%x2=second sequence over n2(n2 can be different from n1)%n=min(min(n1),min(

20、n2):max(max(n1),max(n2);y1=zeros(1,length(n);y2=y1;y1(find(n>=min(n1)&(n<max(n1)= =1)=x1;y2(find(n>=min(n2)&(n<=max(n2)= =1)=x2;y=y1+y24信号乘sigmult.mfunctiony,n=sigmult(x1,n1,x2,n2)%implements y(n)=x1(n)*x2(n)%-%y,n=sigmult(x1,n1,x2,n2)%y=product sequence over n ,which includes n1

21、 and n2%x1=first sequence over n1%x2=second sequence over n2(n2 can be different from n1)%n=min(min(n1),min(n2):max(max(n1),max(n2);y1=zeros(1,length(n);y2=y1;y1(find (n>=min(n1)&(n<=max(n1)= =1)=x1;y2(find(n>=min(n2)&(n<=max(n2)= =1)=x2;y=y1.*y25移位sigshift.m functiony,n=sigshift

22、(x,m,n0)%implements y(n)=x(n-n0)%-%y,n=sigshift(x,m,n0)%n=m+n0;y=x;6折叠sigfold.m functiony,n=sigfold(x,n)%implements y(n)=x(-n)%-%y,n=sigsfold(x,n)%y=fliplr(x);n=-fliplr(n);7 奇偶综合evenodd.mfunction xe,xo,m=evenodd(x,n)%Real signal decomposition into even and odd parts%-%xe,xo,m=evenodd(x,n)%if any (im

23、ag(x)=0) error(x is not a real sequence)end m=-fliplr(n);m1=min(m,n);m2=max(m,n);m=m1:m2;nm=n(1)-m(1);n1=1:length(n);x1=zeros(1,length(m);x1(n1+nm)=x;x=x1;xe=0.5*(x+fliplr(x);xo=0.5*(x-fliplr(x);【例2-1】实验任务3.（1）产生下列序列并绘出离散图：x(n)=2(n+2)-(n-4) -5n5解： 利用MATLAB及信号处理工具箱函数，再加上前面构造的几个函数如sigadd,sigmult等，可很容易

24、编写出可直接执行的MATLAB程序ex31.m:%Example 3.1%a) x(n)=2*delta(n+2)-delta(n-4), -5<=n<=5%figure(1);clfn=-5:5;x=2* impseq (-2,-5,5)-impseq(4,-5,5);subplot(2,2,1);stem(n,x);title (Sequence in Example 3.1a);【例2-2】 设线性时不变（LTI）系统的单位冲激响应为 h(n)=(0.9)nu(n)输入序列为 x(n)=u(n)-u(n-10) 求系统输出y(n)。解： 系统输出y(n)为输入x(n)与系统单

25、位冲激响应h(n)的卷积，可直接采用conv-m函数求出输出序列。MATLAB程序为ex22.m:%Example2-2%x(n)=u(n)-u(n-10);h(n)=(0.9)n*u(n) n=0:04%y(n)=conv(x,h)%n=-5:50;u1=stepseq(0,-5,50);u2=stepseq(10,-5,50);%input x(n)x=u1-u2;%impulse response h(n)h=(0.9).n).*u1; figure(1)subplot (3,1,1);stem(n,x);axis(-5,50,0,2)title(Input Sequence)ylabe

26、l(x(n)subplot(3,1,2);stem(n,h);axis(-5,50,0,2)ylabel(h(n)%output responsey,ny=conv_m(x,n,h,n);subplot(3,1,3);stem(ny,y);axis(-5,50,0,8)gtext(Impulse Response)xlabel(n),ylabel(y(n)执行结果分析所得的序列的离散图.实验报告要求1. 实验报告包括目的、要求、内容、步骤、结果、总结，形成完整实验报告。实验步骤不是书本内容的复制，而是自己结合实验内容进行探索的过程。2. 提交打印版实验报告，A4纸张打印，并附源程序清单。3.

27、提交的报告和源程序中标注清楚姓名学号专业等基本信息，禁止抄袭。实验三 离散系统分析实验目的：通过用Matlab实现信号的分析和表示方法，掌握离散系统的分析方法。实验任务：1. 了解系统的Matlab描述和转换2. 掌握分析线性时不变系统的方法，并编程实现。给定系统H(z)=-0.2z/(z2+0.8),(1) 求出H(z)的幅频响应和相频响应；(2) 绘制极零点图；(3) 求出并绘出该系统的单位抽样响应；(4) 令x(n)=u(n)，求出并绘出系统的单位阶跃响应y(n)。MATLAB系统的描述和转换MATLAB的信号处理工具提供了几种线性时不变系统的模型，我们可以灵活的选择合适的模型以满足不同

28、应用场合的需要。信号处理工具箱提供的离散时间系统模型常用于数字滤波器，这些模型包括传递函数、零极点增益形式、状态空间形式、部分分式、二阶级联形式、格状结构和卷积矩阵等。1 传递函数形式在MATLAB信号处理工具箱中，数字滤波器的传递函数具有以下形式：这种形式称为DSP（Digital Signal Processing）,它和控制工具箱的z传递函数模型形式不同。这一模型用z传递函数分子和分母多项式的系数构成两个向量来确定，即 num=b(1) b(2)b(nb+1)den=a(1) a(2)a(na+1)式中nb和na分别为分子num和den的阶数在MATLAB中，要定义一个DSP形式的传递函

29、数可利用FILT。调用格式： sys=filt(num,den) sys=filt(num,den,Ts)2 零极点增益形式 在MATLAB信号处理工具箱中，离散时间系统（数字滤波器）的零极点形式具有以下形式：这一模型用零点向量和极点向量及增益来表示，即 q(z)=q(1) ,q(2),q(n)p(z)=p(1) ,p(2),p(n)MATLAB函数poly和roots可实现传递函数形式和零极点增益形式的转换。3 状态空间形式离散系统的状态空间表达式为： x(n+1)=Ax(n)+Bu(n) y(n)=Cx(n)+Du(n)式中，u为输入，y为输出,x为状态向量。在MATLAB信号处理工具箱中

30、，用矩阵A、B、C、D表示数字滤波器模型（和控制工具箱一样）MATLAB信号处理工具箱还提供传递函数、零极点和状态空间三种模型之间的相互转换函数，它们是zp2ft、zp2ss 、tf2ss、 tf2zp、ss2tf和 ss2zp4 部分分式离散时间系统（数字滤波器）的传递函数H(z)也可以用下面部分分式展开来描述：（1）若H(z)部包含重极点（2）若H(z)包含多重极点，如p(j)有Sr个重极点，H(z)部分分式展开式中包含相应的形式为：MATLAB信号处理工具箱中，函数RESIDUEZ()用来把传递函数形式转变为部分分式展开式。调用格式为：r,p,k=residuez(b,a) 其中，b,a

31、分别为z传递函数的分子、分母；p为传递函数极点向量；r为相应极点的留数向量；k为传递函数商向量。若原传递函数为真分数，此项为空向量。Residuez也可用于将部分分式展开式转换为传递函数形式，调用格式为：b,a=residuez(r,p,k) 可用函数residuez将H(z)展成部分分式的形式，然后利用公式直接求H(z)的z反变换。5 二阶级联形式任何传递函数可写为二阶因子级联形式，即即H(z)可以化成多个二阶因式相乘，L为二阶因式的个数。这样一个离散系统可用一个L×6的数组SOS表示，每一行表示一个二阶因子，分子三个系数和分母三个系数。 这种形式用于滤波器的二阶节级联结构中。MA

32、TLAB中提供了一般模型和SOS形式的相互转换函数zp2sos、sos2zp、sos2os、sos2tf、6 格状结构滤波器的实现通常采用格状结构。对于N阶全零点和全极点滤波器采用格型结构，可以用多项式a(n),n=1,2,N来描述，也可以通过格型结构反射系数k(n)，n=1,2,.,N来描述。 MATLAB信号处理工具箱提供了滤波器格状结构和转换函数latcfilt、poly2rc、rc2poly、latc2tf和tf2latc。7 卷积矩阵若两个向量x和h，满足下面等式的矩阵C称为卷积矩阵。C·x=h*x （式36）当滤波器的单位采样序列为h(n)，滤波器输入为x(n)，滤波器输

33、出信号为： y(n)=h*x=C·x可见，卷积矩阵是数字滤波器的另一种描述形式。MATLAB信号处理工具箱提供了计算卷积矩阵函数convmtx。调用格式为： C=convmtx(h,n)式中，h为m×1列向量；n为x列向量维数；c为(m+n-1)×n阶卷积矩阵。【例31】已知滤波器的传递函数 试求该滤波器的零极点和状态空间表达式和部分分式展开式。解：用MATLAB编程如下%MATLAB Program 3-1%input parametersnum=2 3;den=1 0.4 1;f=filt(num,den)disp(Convert to zero-pole-g

34、ain)z,p,k=tf2zp(num,den)disp(Convert to state space)A,B,C,D=tf2ss(num,den)transform partial-fraction expansionr,p1,k1=residuez(num,den)MATLAB系统分析MATLAB信号处理工具箱提供了许多工具函数来分析滤波器（系统）的特性，包括时域分析，对任意输入响应、脉冲响应等；频域分析，幅值响应、相位响应、零极点位置；其他特性，群延迟等。滤波器时域和频域分析是设计各类滤波器、评价滤波器性能和滤波器应用基础。这里简单介绍一下有关函数。一时间响应分析工具filter,fil

35、tic,fftfilt和impulse1 Filter（）函数filter用于实现IIR和FIR滤波器对数据滤波，常用来计算滤波器对输入的响应，函数调用格式为：y=filter(b,a,x)其中，b，a分别是滤波器系统函数H(z)的分子和分母多项式系数向量；x为滤波器输入向量或矩阵； 函数filter可用于全极点、全零点、零极点系统。2 Filtic（）函数filtic用于从系统过去值y和x求系统状态初始值，调用格式：z=filtic(b,a,y,x)z=filtic(b,a,y)其中，y=y(n-1),y(n-2),.,y(-na)为输出y的过去向量；x=x(n-1),x(n-2),.,x(

36、-nb)为输出x的过去向量;b,a分别为分子、分母多项式系数向量；z为系统初始状态。如果，n<=-1时，x(n)=0,则filtic函数中不用指定x。【例32】求解 其中 ，初始条件为y(-1)=4和y(-2)=10。解：用MATLAB计算全响应：n=0:7; x=(1/4).n; xic=1,2;b=1;a=1,1.5,0.5;format longy1=filter(b,a,x,xic);Matlab提供了filtic函数计算对于上例：Y=4,10;xic=filter(b,a,Y)3Fftfilt（）函数fftfilt利用效率高的基于FFT重叠相加法实现对数据滤波，该函数只适用于F

37、IR滤波器，函数调用格式为：y=fftfilt(b,x)y=fftfilt(b,x,n)式中，b为FIR滤波器系数向量；x为输入数据；n为FFT长度，省时，函数选择最佳的FFT长度；y为滤波器输出。Y=fftfilt(b,x)等价于y=filter(b,1,x)。【例33】：FIR低通滤波器截止频率为200Hz，采样频率Fs=1000Hz，对信号x(t)=sin(2f1t)+sin(2f2t)滤波，f1=50Hz,f2=250Hz,求滤波器输出。%Matlab Program 3.2%Create Orignal signal dataN=1000;Fs=1000;n=0：N-1;t=n/.F

38、s;x=sin(2*pi*50*t)+sin(2*pi*240*t);%Create a FIR filter with lawpassb=fir1(40,200/500);%Filtering the datayfft=fftfilt(b,x,256);Outputn1=81:241;t1=t(n1);x1=x(n1);y1=yfft(n1);subplot(221);plot(t1,x1);tilte(Original Signal);subplot(222);plot(t1,y1);title(Signal after the filter);grid;执行上面这个程序，分析程序结果。4

39、Impz（）impz用于产生数字系统的脉冲响应。基本调用格式为：h,t=impz(b,a,n,Fs)b,a分别为系统传递函数分子和分母多项式系数向量；n为采样点数；Fs为采样周期，缺省值1；h为单位脉冲响应向量；t为和h对应的时间向量。当函数输出缺省时，绘制系统脉冲响应图；当n缺省时，函数自动选择n值。在Matlab中，h=impz(b,a,n)等价于下面语句：x=1 zeros(1,n);h=filter(b,a,x);因此，不少场合下，用后一种方法求系统的脉冲响应更为灵活。5 Conv（）线性时不变系统的输入输出关系还可通过单位脉冲响应h(n)表示：y(n)=x(n)*h(n)=所以可直接

40、采用MATLAB中的函数conv或者用可求出带下标的序列卷积的函数conv_m来求系统对输入的响应。二 频率响应1 Freqs()函数freqs用于求模拟滤波器的频率响应。调用格式为h=freqs(b,a,)； h,=freqs(b,a,n); freqs(b,a)其中，b,a分别为模拟滤波器传递函数分子和分母多项式向量；n为频率点数(整数),缺省200；h为频率响应,复数；频率向量，实数。函数缺省时，绘制模拟滤波器的幅频响应图和相频响应图。模拟滤波器的系统函数形式为输出h为模拟滤波器的频率响应H(ej)。2Freqz（）函数freqz用于求数字滤波器的频率响应。调用格式为h,=freqz(b

41、,a,n); h,f=freqz(b,a,n,Fs);h,=freqz(b,a,n,whole); h,f=freqz(b,a,n,whole,Fs);h=freqz(b,a,f,Fs); freqz(b,a);说明：（1）函数输入：b,a分别为数字滤波器传递函数分子和分母多项式向量；n为频率点数(整数),最好为2的幂，缺省512；Fs为采样频率，Hz；f为给定的频率矢量；whole 表示返回的值为包含z整个单位圆频率矢量，即02；缺省时，仅为包含z平面上半单位圆（0）之间等间距N个点频率矢量。函数返回：h为频率响应,复数；为n点频率向量（弧度），返回范围与输入参量whole有关；f为n点频率

42、向量(Hz)。函数缺省时，绘制模拟滤波器的幅频响应图和相频响应图。（2）该函数适用于和式31形式的数字滤波器。（3）该函数执行补零的FFT算法，即h=fft(b,n)./fft(a,n)（4）Matlab提供两个函数ABS和ANGLE由频率响应H(ej)求幅频响应| H(ej)|和相频响应argH(ej).（5）函数freqz输出的频率向量在02之间，为了获得一个滤波器真正的相频特性图，要对相位角进行修正。为此Matlab信号处理工具箱提供一个函数UNWRAP.3 UNWRAP.函数UNWARP展开由函数freqz产生的频率.调用格式为：P=unwrap()三 零极点图系统零极点位置决定了系统

43、稳定性和性能，因此，考察系统零极点位置是分析系统特性重要方面之一。Matlab信号处理工具箱提供绘制线性离散时间系统零极点图的工具函数zplane.调用格式为： zplane(z,p) zplane(b,a)其中，z,p分别为系统零、极点向量；b,a为系统z传递函数分子、分母多项式系数向量。四 群延迟 由信号传输不失真条件，滤波器相频特性应该是一条经过原点直线，即 ()= td td为常数 （式38）但一般实际系统不满足上述条件，衡量实际滤波器相位平均延迟用群延迟。群延迟定义为信号通过系统的平均延迟对频率的函数，即滤波器相频特性图上切线的负斜率 （式39）Matlab信号处理工具箱提供计算群延

44、迟函数GRPDELAY，调用格式为：gd,=grpdelay(b,a,n); gd,f=grpdelay(b,a,n,Fs);gd,=grpdelay(b,a,n,whole);gd,f=grpdelay(b,a,n,wholeFs);其中，gd为群延迟，其他各项意义同函数freqz。函数输出项缺省时绘制群延迟图。【例34】：编一函数实现：求系统的绝对的和相对的幅度响应，相位响应及群延迟响应。解：函数freqz_m是对freqz的修正，它给出了系统的绝对的和相对的幅度响应，相位响应及群延迟响应。程序如下：function db,mag,pha,grd,w=freqz_m(b,a)% db,ma

45、g,pha,grd,w=freqz_m(b,a)%db=0pi区间内的相对振幅mag=0pi区间内的绝对振幅pha=0pi区间内的相位响应grd=0pi区间内的群延迟w=0pi区间内的501个频率样本向量b=H(z)的分子多项式的系数（对FIR：b=h）a=H(z)的分母多项式的系数(对FIR:a=1)H,w=freqz(b,a,1000,whole);H=(H(1:1:501);w=(w(1:1:501);mag=abs(H);db=20*log10(mag+eps)/max(mag);pha=angle(H); grd=grpdelay(b,a,w);【例35】 已知滤波器的传递函数 H(

46、z)=0.15/(1-0.8z-1)输入信号为x(t)=2sin(0.05t)+w(t)，w(t)为随机信号，幅值为0.2。试绘出滤波器的输出信号波形。用两种方法求滤波器的输出，MATLAB编写程序如下：%MATLAB Program 3-3b=0.15;a=1 0.8;%create input of filter n=0:100;x=2*sin(0.05*pi*n)+0.2*randn(1,101);%Calculate response according to convolutionimp=1,zeros(100,1);h=impz(b,a,101);yc=conv(h,x);yc=y

47、c(1:101);%Calculate response according to matlab functiony1=filter(b,a,x);plot(n,x,r,n,y1,b,n,y,m);xlabel(n);ylabel(x y yc);title(Time response);grid;分析程序执行结果。实验报告要求1. 实验报告包括目的、要求、内容、步骤、结果、总结，形成完整实验报告。实验步骤不是书本内容的复制，而是自己结合实验内容进行探索的过程。2. 提交打印版实验报告，A4纸张打印，并附源程序清单。3. 提交的报告和源程序中标注清楚姓名学号专业等基本信息，禁止抄袭。实验四 离

48、散傅立叶变换实验目的：理解离散傅立叶变换及其性质，与线性卷积和线性移位相比较掌握循环卷积和循环移位的实现方法。实验任务：1. 根据给出的函数，编程实现某一序列DFT, IDFT。已知x1为5点序列，x2为14点序列，各作15点的DFT，然后两个DFT相乘，再求乘积的IDFT，设结果为y(n)，分析y(n)哪些点对应于x1*x2应该得到的点？2. 编程实现某一序列的循环移位和循环卷积。已知x1=(n+1)R4(n)，h(n)=R4(n-2)，试求x1和h(n)的循环卷积，并作图。3. 线性卷积与循环卷积的关系已知x1=1 2 2，x2=1 2 3 4，试计算并分析结果（1）x1x2；(2)x1x

49、2；(3) x1x2；(4)x1*x2.4DFT的奇偶虚实性已知序列x(n)=10(0.8)n，序列长度N=21,绘出x(n)傅立叶变换的实部和虚部并验证其奇偶性。4. 执行exam41例程序，分析其结果。【例 41】为说明高密度频谱和高分辩频谱之间的差异，考虑 x(n)=cos(0.48n)+cos(0.52n) （1）取x(n)(0n10)时，求x(n)的DFT X(k)； （2）将（1）中的x(n)以补零方式使x(n)加长到0n100,求X(k) （3）取x(n)(0n100),求X(k).解 MATLAB程序为 exam41.m:%MATLAB Program 4-1%part 1%s

50、pecttum based on the first 10 amples of x(n)n=0:1:10;x=cos(0.48*pi*n)+cos(0.52*pi*n);figure(1)n1=0:1:9;y1=x(1:1:10);subplot(2,1,1);stem(n1,y1);title('signal x(n),0<=n<=9');xlabel('n')%axis(0,10,-4,4);Y1=fft(y1);magY1=abs(Y1(1:1:6);k1=0:1:5;w1=2*pi/10*k1;Subplot(2,1,2);stem(w1/p

51、i,magY1);title('Samples of DTET Magnitude');Xlabel('frequency in pi units')axis(0,1,0,10)%part 2%High density spectrum (100 samples )based on the first 10 samples of x(n)figure(2)n3=0:1:99;y3=x(1:1:10)zeros(1,90);subplot(2,1,1);stem(n3,y3);title('signal x(n),0<=n<=9+90 zeros');xlabel('n')axis(0,100,-2.5,2.5);Y3=fft(y3);magY3=abs(Y3(1:1:51);k3=0:1:50;w3=2*pi/100*k3;subplot(2,1,2);stem(w3/pi,magY3