70华北电力大学Matla实验指导书

数字信号处理实验指导书（基于MATLAB）尚秋峰宋文妙华北电力大学二○○六年十一月前言MATLAB是一套功能强大的工程计算及数据处理软件，广泛应用于工业，电子，医疗和建筑等众多领域。它是一种面向对象的，交互式程序设计语言，其结构完整，又具有优良的可移植性。它在矩阵运算，数字信号处理方面有强大的功能。另外，MATLAB提供了方便的绘图功能，便于用户直观地输出处理结果。本课程实验要求学生掌握用MATLAB语言仿真数字信号处理算法的基本方法，加深对教学内容的理解。说明：本书给出的为本课程基本实验，教师可根据学生情况安排内容延伸实验。目录实验一序列的产生和运算...........................................1实验二因果性数字系统的时域实现...................................5实验三离散傅里叶变换（DFT）及其快速算法（FFT）...................10实验四FFT的典型应用.............................................17实验五FIR滤波器设计.............................................23实验六IIR滤波器设计.............................................29附录............................................................321实验一序列的产生和运算一、实验目的：1、熟悉MATLAB环境，掌握基本编程方法。2、熟悉MATLAB中序列产生和运算的基本函数。二、实验内容：1、MATLAB入门（1）了解MATLAB的工作窗口进入MATLAB环境，选择View中的各个选项，可以打开命令窗口CommandWindow、变量浏览器Workspace、路径浏览器CurrentDirectory，如下图示：图1－1Matlab工作窗口MATLAB命令窗口，是键入指令的地方也是MATLAB显示计算结果的地方。在符号之后键入“1+2+3”，或者“x=1+2+3”，显示结果有什么不同？如果在上述的例子结尾加上“;”，则计算结果不会在窗口自动显示，要显示计算结果须键入该变量x。进行以上操作后打开变量浏览器，观察里面的变量名和变量的值。在MATLAB环境选择File再选择New，即进入程序编辑/调试窗口，如下图示。存档时必须以.m名称储存。要执行M-file可以在命令窗口下直接键入该文件名；或是选择Debug2下的RunM-file来执行M-file。如果要修改M-file可以选择功能表上的OpenM-file，即可搜寻要修改的M-file，修改后再存档。尝试编写程序文件，完成以上操作。图1－2程序编辑/调试窗口练习使用路径浏览器打开特定目录下的M-file。M-file还可以用来定义函数，然后储存起来，就可以和那些内建的函数（如sin,cos,log等）一样的自由使用。举例来说，我们可以定义一函数cirarea来计算圆的面积，以下的M-file:cirarea.m就是定义这个函数的：%M-filefunction,cirarea.m%Calculatetheareaofacirclewithraduisr%rcanbeascalaroranarrayfunctionc=cirarea(r)c=pi*r.^2;注意，M-file定义的函数语法上有一些规定：第一行指令以function这个字做为起头，接着是输出的变量，等号，函数名称，输入的变量是放在括号之内。functionout1=userfun(in1)，这行的out1是输出的变量，userfun是函数名称，in1是输入的变量。function[out1,out2]=serfun(in1,in2)如果输出变量[out1,out2]和输入变数(in1,in2)不只一个时，则在输出变量部份须加上[]。上述的输入变量是在调用函数时输入的，而输出的变量即是函数的返回值。函数名称的取法的规定与一般变量相同。在定义函数之前，最好加上注解行来说明这个函数的特色及如何使用，这样，使用指令如helpcirarea，该函数的注解行会出现在指令视窗。3尝试编写函数文件，并在你编写的程序文件中应用此函数。MATLAB会将绘图结果展示在另一个视窗，称为MATLABFigureWindows，如下图示。如果你看不到此视窗，可以进入Windows再选择Figure。在图形窗口中，可以利用Edit菜单中的选项来改变显示效果以及拷贝图形，尝试这些操作。图1－3MATLAB图形视窗（2）寻求帮助在MATLAB系统中相关的线上(on-line)求助方式有三种：利用help指令，如果你已知要找的主题为何的话，直接键入helptopic。所以即使身旁没有使用手册，也可以使用help指令查询不熟悉的指令或是题材之用法，例如helpsqrt,helptopic。利用lookfor指令，它可以从你键入的关键字(key-word)（即使这个关键字并不是MATLAB的指令）列出所有相关的主题，例如lookforcosine,lookforsine。利用指令视窗的功能选单中的Help，从中选取TableofContents（目录）或是Index（索引）。练习以上寻求帮助的方法。2、波形产生学习附录中的有关知识，完成下列练习。练习一：建立一个MATLAB函数impseq.m，用来生成迟延的单位脉冲序列。其输入参数为：序列起始位置ns，序列终止位置nf，脉冲位置np。练习二：编写MATlAB程序，以产生下图所示的波形。并将序列绘制出来。4练习三：设x＝[2345]，求将它延拓6个周期所得的序列。练习四：设x＝[2345]，实现以下运算并绘制波形练习五：编写程序实现下列函数，并绘制波形以t以t=0.01（n=0:N-1）进行取样，N=64三、思考题：算术运算符*和．*之间的区别是什么？)(8.0)8sin(5)4sin(2)(twtttx()(1,)wnrandnN5实验二因果性数字系统的时域实现一、实验目的编制实现下列差分方程的程序:y(n)=∑bkx(n-k)+∑aky(n-k)二、实验内容1、编制nonrec.m函数文件,实现FIR滤波y(n)=h(n)*x(n)。这里给定h(n)=R8(n)，x(n)=nR16(n)，求y(n)。nonrec.m函数文件:functiony=nonrec(x,h)x=[x,zeros(1,length(h)-1)];%补零w=zeros(1,length(h));fori=1:length(x)forj=length(h):-1:2w(j)=w(j-1);%得到每一延时单元输出变量endw(1)=x(i);y(i)=w*h.’;%h与w对应相乘end主程序文件:x=0:15;h=ones(1,8);y=nonrec(x,h);n=0:22;stem(n,y);图2－1卷积运算实现FIR滤波6分析：线性卷积y(n)=x(n)*h(n)的长度为16+8-1=23，可利用y(n)=∑h(m)x(n-m)直接计算得n(n+1)/2,n≤7y(n)=4(2n-7),8≤n≤15(n+8)(23-n)/2,16≤n≤22即y=[0136101521283644526068768492847554422915],与曲线相符。2、编制rec.m函数文件,实现纯递归IIR滤波y(n)=x(n)+∑aky(n-k)。这里给定a1=2rcosw0,a2=-r2,r=0.95,w0=π/8,求单位抽样响应h(n).rec.m函数文件:functiony=rec(x,a,n)x=[x,zeros(1,n-length(x))];%补零到所需长度sum=0;w=zeros(1,length(a));fori=1:ny(i)=sum+x(i);%递归forj=length(a):-1:2w(j)=w(j-1);%延时endw(1)=y(i);sum=w*a';end主程序文件:x=[1];a=[2*0.95*cos(pi/8),-0.95^2];h=rec(x,a,75);%取h(n)的长度为75点n=0:74;stem(n,h);图2－2纯递规IIR滤波分析计算:7由题意,a1=2*0.95*cos(π/8),a2=-0.952,所以,得到系统函数H(z)=1/[1-1.9cos(π/8)z-1+0.952z-2],做逆Z变换得h(n)=0.95ncos(πn/8)+ctg(π/8)*0.95nsin(πn/8),利用MATLAB直接画h(n),即,使用下列语句n=0:74;h=0.95.^n.*cos(pi.*n./8)+cot(pi/8).*(0.95.^n).*sin(pi.*n./8);stem(n,h);得到如图2－3的曲线图2－3此曲线与用rec函数所画曲线完全一致,可见,用MATLAB编制的FIR滤波程序与理论计算所得结果是相同的.3、用nonrec.m和rec.m函数编制dfilter.m函数文件,构成完整的一般IIR滤波程序,并完成下列信号滤波:x(n)=cos(2πn/5)R64(n)这里给定系统函数H(z)=(1-2z-1+z-2)/(1-0.5z-1+0.5z-2)求y(n).dfilter.m函数文件:functiony=dfilter(x,b,a,n)y1=nonrec(x,b);y=rec(y1,a,n);主程序文件:n=0:63;x=cos(2*pi/5*n);b=[1,-2,1];%由H(z)得到系数a,b8a=[0.5,-0.5];y=dfilter(x,b,a,64);%取y(n)的长度为64点stem(n,y);图2－44．用help查看内部函数filter.m,了解调用格式,重做3,并与你编写的dfilter.m结果进行比较.用help可以看到内部函数为Y=FILTER(B,A,X),且有“ThefilterisaDirectFormIITransposedimplementationofthestandarddifferenceequation”:a(1)*y(n)=b(1)*x(n)+b(2)*x(n-1)+...+b(nb+1)*x(n-nb)-a(2)*y(n-1)-...-a(na+1)*y(n-na)因此,调用内部函数filter时,要对原系数a做适当变化:a1(1)=1,a1(i)=-a(i-1),i1.n=0:63;x=cos(2*pi/5*n);b=[1,-2,1];a=[0.5,-0.5];y=dfilter(x,b,a,64);%调用自己编的dfilter函数a1=[1,-0.5,0.5];%a变为a1,用于调用内部函数filtery1=filter(b,a1,x);subplot(2,1,1);stem(n,y);subplot(2,1,2);stem(n,y1);9图2－5三、思考题1、内部函数filter.m的调用格式是什么?与你编写的dfilter.m调用格式是否一致?差异在何处?2、实验中如何合理地选取单位抽样响应h(n)的点数?实验内容及要求中第三项的物理概念是什么?给定的H(z)是具有什么性质的滤波器?3、为实现各实验内容中的滤波器，你编写程序时采用的什么结构？4、试利用Matlab中的卷积函数conv完成实验内容1,并与本书的例程进行比较。10实验三离散傅里叶变换（DFT）及其快速算法（FFT）一、实验目的1、学习DFT和FFT的原理及编程实现方法。2、熟悉MATLAB编程的特点。二、实验内容1、用三种不同的DFT程序计算x(n)=R8(n)的傅里叶变换X(ej)，并比较三种程序计算机运行时间。（1）用forloop语句的M函数文件dft1.m，用循环变量逐点计算X(k)；（2）编写用MATLAB矩阵运算的M函数文件dft2.m,完成下列