MATLAB实验指导书

1、数字信号处理实验指导书（基于matlab）华北电力大学二六年十一月前 言matlab 是一套功能强大的工程计算及数据处理软件，广泛应用于工业，电子，医疗和建筑等众多领域。它是一种面向对象的，交互式程序设计语言，其结构完整，又具有优良的可移植性。它在矩阵运算，数字信号处理方面有强大的功能。另外，matlab提供了方便的绘图功能，便于用户直观地输出处理结果。本课程实验要求学生掌握用matlab语言仿真数字信号处理算法的基本方法，加深对教学内容的理解。说明：本书给出的为本课程基本实验，教师可根据学生情况安排内容延伸实验。目 录实验一 序列的产生和运算1实验二 因果性数字系统的时域实现5实验三 离散傅

2、里叶变换（dft）及其快速算法（fft）10实验四 fft的典型应用17实验五 fir滤波器设计23实验六 iir滤波器设计29附 录32 实验一 序列的产生和运算一、实验目的：1、熟悉matlab环境，掌握基本编程方法。2、熟悉matlab中序列产生和运算的基本函数。二、实验内容：1、matlab入门（1）了解matlab 的工作窗口进入matlab环境，选择view中的各个选项，可以打开命令窗口command window、变量浏览器workspace、路径浏览器current directory，如下图示：图11 matlab 工作窗口matlab命令窗口，是键入指令的地方也是matla

3、b显示计算结果的地方。在符号之后键入“1+2+3”，或者“x=1+2+3”，显示结果有什么不同？如果在上述的例子结尾加上“;”，则计算结果不会在窗口自动显示，要显示计算结果须键入该变量x。进行以上操作后打开变量浏览器，观察里面的变量名和变量的值。在matlab环境选择file 再选择new，即进入程序编辑/调试窗口，如下图示。存档时必须以.m名称储存。要执行 m-file 可以在命令窗口下直接键入该文件名；或是选择debug下的run m-file来执行m-file。如果要修改 m-file 可以选择功能表上的open m-file ，即可搜寻要修改的 m-file，修改后再存档。尝试编写程序

4、文件，完成以上操作。图12 程序编辑/调试窗口练习使用路径浏览器打开特定目录下的m-file。m-file还可以用来定义函数，然后储存起来，就可以和那些内建的函数（如sin, cos,log等）一样的自由使用。举例来说，我们可以定义一函数cirarea来计算圆的面积，以下的 m-file: cirarea.m就是定义这个函数的：% m-file function, cirarea.m % calculate the area of a circle with raduis r % r can be a scalar or an array function c=cirarea(r) c=pi*

5、r.2;注意，m-file定义的函数语法上有一些规定： 第一行指令以function这个字做为起头，接着是输出的变量，等号，函数名称，输入的变量是放在括号之内。function out1=userfun(in1)，这行的out1是输出的变量，userfun是函数名称，in1是输入的变量。function out1, out2=serfun(in1, in2) 如果输出变量 out1,out2和输入变数(in1, in2)不只一个时，则在输出变量部份须加上 。 上述的输入变量是在调用函数时输入的，而输出的变量即是函数的返回值。 函数名称的取法的规定与一般变量相同。 在定义函数之前，最好加上注解行

6、来说明这个函数的特色及如何使用，这样，使用指令如help cirarea，该函数的注解行会出现在指令视窗。尝试编写函数文件，并在你编写的程序文件中应用此函数。matlab会将绘图结果展示在另一个视窗,称为matlab figure windows，如下图示。如果你看不到此视窗，可以进入windows再选择figure。在图形窗口中，可以利用edit菜单中的选项来改变显示效果以及拷贝图形，尝试这些操作。图13 matlab 图形视窗（2）寻求帮助在matlab系统中相关的线上(on-line)求助方式有三种： 利用help指令，如果你已知要找的主题为何的话，直接键入help 。所以即使身旁没有使

7、用手册，也可以使用help指令查询不熟悉的指令或是题材之用法，例如help sqrt, help topic。 利用lookfor指令，它可以从你键入的关键字(key-word)（即始这个关键字并不是matlab的指令）列出所有相关的主题，例如lookfor cosine, lookfor sine。 利用指令视窗的功能选单中的help，从中选取table of contents（目录）或是index（索引）。 练习以上寻求帮助的方法。2、波形产生学习附录中的有关知识，完成下列练习。练习一：建立一个matlab函数impseqm，用来生成迟延的单位脉冲序列。其输入参数为:序列起始位置ns，序列

8、终止位置nf，脉冲位置np。function x,n=impseq(np,ns,nf)n=ns:nf;x=(n-n0)=0;n=0:5;x=impseq(3,0,5);subplot(2,2,1);stem(n,x);xlabel(n);ylabel(xl(n);练习二：编写matlab程序，以产生下图所示的波形。并将序列绘制出来。 n=0:29;x1=2.5*(stepseq(0,0,9)-stepseq(6,0,9);x2=-2.5*(stepseq(6,0,9)-stepseq(10,0,9);x=x1+x2;x3=x,x,x;subplot(2,2,1);stem(n,x3);xlab

9、el(n);ylabel(x(n);n=0:29;x1=2.5*(stepseq(0,0,9)-stepseq(3,0,9);x2=-2.5*(stepseq(3,0,9)-stepseq(10,0,9);x=x1+x2;x3=x,x,x;subplot(2,2,2);stem(n,x3);xlabel(n);ylabel(x(n);练习三：设x2 3 4 5，求将它延拓6个周期所得的序列。x=2,3,4,5;y=x,x,x,x,x,x;subplot(2,2,1);stem(y);xlabel(n);ylabel(x(n);练习四：设x2 3 4 5，实现以下运算并绘制波形练习五： 编写程序

10、实现下列函数，并绘制波形以 t以t=0.01（n=0:n-1）进行取样，n=64三、思考题：算术运算符*和*之间的区别是什么？ 标量和矢量实验二 因果性数字系统的时域实现一、 实验目的编制实现下列差分方程的程序: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函数文件:function y=nonrec(x,h)x=x,zeros(1,length(h)-1); %补零w=zeros(1,length(h);for i=1

11、:length(x) for j=length(h):-1:2w(j)=w(j-1); %得到每一延时单元输出变量 end w(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);图21 卷积运算实现fir滤波分析：线性卷积y(n)=x(n)*h(n)的长度为16+8-1=23，可利用y(n)=h(m)x(n-m)直接计算得 n(n+1)/2, n7y(n)= 4(2n-7), 8n15 (n+8)(23-n)/2, 16n22即 y= 0 1 3 6 10 15 21 28

12、 36 44 52 60 68 76 84 92 84 75 54 42 29 15 ,与曲线相符。2、编制rec.m函数文件,实现纯递归iir滤波 y(n)=x(n)+aky(n-k)。 这里给定a1=2rcosw0, a2=-r2, r=0.95, w0=/8, 求单位抽样响应h(n).rec.m函数文件:function y=rec(x,a,n)x=x,zeros(1,n-length(x); %补零到所需长度sum=0;w=zeros(1,length(a);for i=1:n y(i)=sum+x(i); % 递归 for j=length(a):-1:2 w(j)=w(j-1);

13、%延时 end w(1)=y(i); sum=w*a;end主程序文件:x=1;a=2*0.95*cos(pi/8),-0.952;h=rec(x,a,75); %取h(n)的长度为75点 n=0:74;stem(n,h);图22 纯递规iir滤波分析计算:由题意, 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

14、(pi.*n./8)+cot(pi/8).*(0.95.n).*sin(pi.*n./8);stem(n,h); 得到如图23的曲线图23此曲线与用rec函数所画曲线完全一致, 可见, 用matlab编制的fir滤波程序与理论计算所得结果是相同的.3、用nonrec.m和rec.m函数编制dfilter.m 函数文件, 构成完整的一般iir滤波程序, 并完成下列信号滤波:x(n)=cos(2n/5)r64(n) 这里给定系统函数 h(z)=(1-2z-1+z-2)/(1-0.5z-1+0.5z-2) 求y(n).dfilter.m函数文件:function y=dfilter(x,b,a,n)

15、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,ba=0.5,-0.5;y=dfilter(x,b,a,64); %取y(n)的长度为64点stem(n,y);图244用help查看内部函数filter.m, 了解调用格式, 重做3, 并与你编写的dfilter.m结果进行比较.用help可以看到内部函数为y = filter(b,a,x), 且有 “the filter is a direct form ii transposed implementation of the st

16、andard difference equation”: 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

17、=filter(b,a1,x);subplot(2,1,1);stem(n,y);subplot(2,1,2);stem(n,y1);图25三、 思考题1、 内部函数filter.m的调用格式是什么?与你编写的dfilter.m调用格式是否一致?差异在何处? y = filter(b,a,x) 有差异 y=dfilter(x,b,a,n) 一个定义了y（n）长度 ，一个没有定义 但对系数a做适当变化2、 实验中如何合理地选取单位抽样响应h(n)的点数?实验内容及要求中第三项的物理概念是什么?给定的h(z)是具有什么性质的滤波器?3、 为实现各实验内容中的滤波器，你编写程序时采用的什么结构？4、

18、 试利用matlab中的卷积函数conv完成实验内容1,并与本书的例程进行比较。实验三 离散傅里叶变换（dft）及其快速算法（fft）一、实验目的1、学习dft和fft的原理及编程实现方法。2、熟悉matlab编程的特点。二、实验内容1、用三种不同的dft程序计算x(n)=r8(n)的傅里叶变换x (ejww),并比较三种程序计算机运行时间。（1）用for loop 语句的m函数文件dft1.m，用循环变量逐点计算x(k)；（2）编写用matlab矩阵运算的m函数文件dft2.m, 完成下列矩阵运算；（3）调用fft库函数，直接计算x(k)；（4）分别利用上述三种不同方式编写的dft程序计算序

19、列x(n)的傅立叶变换x(ejw)，并画出相应的幅频和相频特性，再比较各个程序的计算机运行时间。（一）实验例程m函数文件如下：dft1.m:functionam,pha=dft1(x)n=length(x);w=exp(-j*2*pi/n);for k=1:n sum=0; for n=1:n sum=sum+x(n)*w(k-1)*(n-1); end am(k)=abs(sum); pha(k)=angle(sum);enddft2.m:functionam,pha=dft2(x)n=length(x);n=0:n-1;k=0:n-1;w=exp(-j*2*pi/n);nk=n*k;wnk

20、=w.(nk);xk=x*wnk;am=abs(xk); pha=angle(xk);dft3.m:functionam,pha=dft3(x)xk=fft(x);am=abs(xk);pha=angle(xk);源程序及运行结果：（1）x=ones(1,8),zeros(1,248);t=cputime;am1,pha1=dft1(x);t1=cputime-tn=0:(length(x)-1);w=(2*pi/length(x)*n;figure(1)subplot(2,1,1), plot(w,am1,b); grid;title(magnitude part);xlabel(frequ

21、ency in radians);ylabel(|x(exp(jw)|);subplot(2,1,2), plot(w,pha1,r); grid;title(phase part);xlabel(frequency in radians);ylabel(argxexp(jw)/radians);图31 dft1记录运行时间：t1 = 0.1719（2）x=ones(1,8),zeros(1,248);t=cputime;am2,pha2=dft2(x);t2=cputime-tn=0:(length(x)-1);w=(2*pi/length(x)*n;figure(2)subplot(2,1

22、,1), plot(w,am2,b); grid;title(magnitude part);xlabel(frequency in radians);ylabel(|x(exp(jw)|);subplot(2,1,2), plot(w,pha2,r); grid;title(phase part);xlabel(frequency in radians);ylabel(argxexp(jw)/radians);图32 dft2 记录运行时间：t2 =0.1406（3）x=ones(1,8),zeros(1,248);t=cputime;am3,pha3=dft3(x);t3=cputime-

23、t;n=0:(length(x)-1);w=(2*pi/length(x)*n;figure(3)subplot(2,1,1), plot(w,am3,b); grid;title(magnitude part);xlabel(frequency in radians);ylabel(|x(exp(jw)|);subplot(2,1,2), plot(w,pha3,r); grid;title(phase part);xlabel(frequency in radians);ylabel(argxexp(jw)/radians);图32 dft2记录运行时间：t3=0从以上运行结果可以看出，调

24、用fft库函数直接计算x(k)速度最快，矩阵运算次之，用循环变量逐点计算运行速度最慢。因此fft算法大大提高了dft的实用性。（二）练习改变输入信号类型和点数，重做例程中的内容，并整理总结实验结果,给出结论。2、实现序列的内插和抽取所对应的傅里叶变换。给定序列x(n)=cos(/36n)+cos(1.5/36n)r128(n)，做n=128点的傅里叶变换，并求x1(n)=x(4n) 和 对应的傅里叶变换(n=128点)。（一）实验例程源程序如下：n=32; %抽取32点t=0:n-1;x2=cos(pi/9*t)+cos(1.5*pi/9*t);xk2=fft(x2,128); %做128点傅

25、里叶变换am2=abs(xk2);n=0:127;w=(2*pi/128)*n;figure(2)plot(w,am2,b); title(magnitude part);xlabel(frequency in radians);ylabel(|x(exp(jw)|);x3=zeros(1,128); %初始化数组为零n=32;for i=0:n-1 %n=4i时赋值 x3(1,4*i+1)=cos(pi/36*i)+cos(1.5*pi/36*i);endxk3=fft(x3); %做128点fftam3=abs(xk3);n=0:127;w=(2*pi/128)*n;figure(3)pl

26、ot(w,am3,b);title(magnitude part);xlabel(frequency in radians);ylabel(|x(exp(jw)|);抽取后的离散傅里叶变换：图34 抽取的离散傅里叶变换内插后的离散傅里叶变换：图33 内插的离散傅里叶变换（二）练习仿照例程编写给定序列x(n)=cos(/36n)+cos(1.5/36n)r128(n)的128点傅里叶变换程序，比较这三个计算结果得到的幅频特性图，分析其差别产生的原因？。三、思考题实验内容4中，分别求x1(ejww)和x2(ejww)与原x(ejww)的关系，并说明为什么抽取后能分辨出两个频率分量？有无任何条件的限

27、制？实验四 fft的典型应用一、实验目的1、学习fft的应用于卷积运算和谱分析时的编程实现方法。2、熟悉matlab编程的特点。二、实验内容1、利用fft实现两序列的卷积运算，并研究fft点数与混叠的关系。给定x(n)=nr16(n)，h(n)=r8(n)，用fft和ifft分别求线性卷积和混叠结果输出，并用函数stem(n,y)画出相应图形。（一）实验例程源程序如下：57 n=16;x=0:n-1;h=ones(1,8);线性卷积：xk1=fft(x,23); %做23点ffthk1=fft(h,23);yk1=xk1.*hk1;y1=ifft(yk1);n=0:22;figure(1)st

28、em(n,y1)运行结果：线性卷积：图31 线性卷积（二）练习改写上面例程做18点fft,由程序运行结果总结前多少（5）个点为混叠结果，后多少（14）个点是线性卷积的结果。如果改做30点fft，又会出现什么结果？依据实验结果总结：要用dft来做线性卷积，dft的点数必须满足什么要求？l=n1+n2-12、研究高密度频谱与高分辨率频谱（一）实验例程设有连续信号xa(t)=cos(26.5103t)+cos(27103t)+cos(29103t)以采样频率fs=32khz对信号xa(t)采样,分析下列三种情况的幅频特性。（1） 采集数据长度n=16点，做n=16点的dft，并画出幅频特性。（2）

29、采集数据长度n=16点，补零到256点，做n=256点的dft，并画出幅频特性。（3） 采集数据长度n=256点，做n=256点的dft，并画出幅频特性。观察三幅不同频率特性图，分析和比较它们的特点以及形成的原因。源程序如下：fs=32000; %采样频率n=16; %采集16点n=0:n-1; %做16点dftxa=cos(2*pi*6.5*103*n/fs)+cos(2*pi*7*103*n/fs)+cos(2*pi*9*103*n/fs);am1,pha1=fft (xa);n=0:(length(xa)-1); w=(2*pi/length(xa)*n;figure(1)plot(w,

30、am1,b); grid;title(magnitude part);xlabel(frequency in radians);ylabel(|x(exp(jw)|);x=xa(1:16),zeros(1,240); %补零到256am2,pha2=fft(x); %做256点dftn=0:(length(x)-1);w=(2*pi/length(x)*n;figure(2)plot(w,am2,b); grid;title(magnitude part);xlabel(frequency in radians);ylabel(|x(exp(jw)|);n0=256; %采集256点n0=0:

31、n0-1; %做256点dftxa0=cos(2*pi*6.5*103*n0/fs)+cos(2*pi*7*103*n0/fs)+cos(2*pi*9*103*n0/fs);am3,pha3=fft(xa0);n=0:(length(xa0)-1);w=(2*pi/length(xa0)*n;figure(3)plot(w,am3,b); grid;title(magnitude part);xlabel(frequency in radians); ylabel(|x(exp(jw)|); 运行结果：n=16点dft幅频特性图32 n=16点的dft高密度频谱：图33 补零到256点的dft

32、高分辨率频谱：图34 256点的dft观察运行结果可知，n=16点时的dft，由于点数太少，很难反映出频谱的细节特征，只能分辨出两个频率分量。将16点dft补零到256点，做256点dft,由结果可以看出，频率间隔缩小（由2p/16减小为2 p/256），得到了比较平滑的连续曲线，此为高密度频谱。它是由16点经过补零所得到的，并没有增加任何信息，频率分辨率并未提高，仍然只能看出两个频率分量。第三次做的是256点的dft，采集点数大大增加，分辨率提高，可以清楚看到三根谱线的位置，是为高分辨率频谱。（二）练习练习一：设(1)取 时，求x(n)的fft变换x(k);(2)将（1）中的 x(n)以补零

33、方式加长到长度为100，求 x(k);（3）取 ，求 x(k)要求：用三个图形窗口；各图形窗口包括两个子图，分别显示x(n)和x(k)的幅值；比较三种情况的频谱成分；说明高密度谱（补0）和高分辨率谱（3）的不同。频率分辨率提高 能够清楚地看到谱线的位置尝试n=120,200等，有什么现象？并解释。点数越多 分辨率越高，可以更清楚的看到谱线的位置。原因：点数越多 得到的信息越多 分辨率提高练习二：模拟信号以t=0.01n（n=0:n-1）进行取样，n=64，试比较有无噪声时的信号谱。尝试增加噪声的幅值。3.重叠保留法做线性卷积x=0:134; %x为长序列y=;h=ones(1,8),zeros

34、(1,8); %h为单位抽样响应h=fft(h); for i=0:14 if (i=0) xk=zeros(1,7),0:8; %x序列的前面补零m-1点(m-1=7) else xk=x(9*i-6:9*i+9); %截取xk(n)中的n点(n=16) end xk=fft(xk); yk=h.*xk; yk=ifft(yk); %利用fft和ifft做循环卷积 y=y,yk(8:16); %截取yk(n)中的未混叠部分并衔接endstem(x,y); 图35 重叠保留法做线性卷积利用matlab内部函数filter对此结果进行验证, 输入语句filter(h,1,x),则可得到如下曲线.

35、 与重叠保留法所得结果相同.图36 直接卷积结果三、思考题1、什么是高密度频谱？什么是高分辨率频谱？实验内容1中为区分三个频率分量至少应采集几个样本？ 2、实验内容2中，分别求x1(ejww)和x2(ejww)与原x(ejww)的关系，并说明为什么抽取后能分辨出两个频率分量？有无任何条件的限制？实验五 fir滤波器设计一、实验目的1、学习fir滤波器的设计方法。2、熟悉matlab编程的特点。二、实验内容1、利用窗函数法设计fir滤波器 设计具有下列指标的低通数字滤波器（线性相位）wp通带截止频率 rp通带衰减 ws阻带截止频率 as阻带衰减 选择合适的窗函数，并画出滤波器的频率响应（db）、

36、窗函数、理想和实际的冲激响应。实验参考例程：%=fir by windowsclear;wp=0.2*pi;ws=0.3*pi;tr_width=ws-wp;n=ceil(6.6*pi/tr_width)+1;n=0:n-1;wc=(ws+wp)/2;hd=ideal_lp(wc,n);w_ham=(hamming(n);h=hd.*w_ham;db,mag,pha,grd,w=freqz_m(h,1);delta_w=2*pi/1000; 取1000点rp=-(min(db（1:wp/delta_w+1); 间隔1dbas=-round(max(db(ws/delta_w+1:501);fi

37、gure(1);subplot(2,2,1);stem(n,hd);subplot(2,2,2);stem(n,w_ham);subplot(2,2,3);plot(w/pi,pha);% stem(n,h);subplot(2,2,4);plot(w/pi,db);figure(2);h,w=freqz(h,1,1000,whole);h=(h(1:501);w=(w(1:501);freqzplot(h,w);各函数文件：function hd=ideal_lp(wc,m)alpha=(m-1)/2;n=0:m-1;m=n-alpha+eps;hd=sin(wc*m)./(pi*m);fu

38、nction db,mag,pha,grd,w=freqz_m(b,a)h,w=freqz(b,a,1000,whole);h=(h(1:501);w=(w(1:501);mag=abs(h);db=20*log10(mag+eps)/max(mag);pha=angle(h);grd=grpdelay(b,a,w);grd=grd;练习一：请给例程加上注释。练习二：选用其它窗函数，修改程序并和例程比较运行结果。练习三：选用fir1（）函数实现题目要求，并和例程比较运行结果。？练习四：分析并运行下列程序，说明该程序实现的是何种功能的滤波器？理想带通滤波器 用布莱克曼窗函数设计 p236%=fi

39、r by windowsclear;wp1=0.35*pi;ws1=0.2*pi;wp2=0.65*pi;ws2=0.8*pi;as=60;tr_width=min(wp1-ws1),(ws2-wp2);m=ceil(11*pi/tr_width)+1;n=0:m-1;wc1=(ws1+wp1)/2;wc2=(ws2+wp2)/2;hd=ideal_lp(wc2,m)-ideal_lp(wc1,m);w_bla=(blackman(m);h=hd.*w_bla;db,mag,pha,grd,w=freqz_m(h,1);delta_w=2*pi/1000;rp=-(min(db(wp1/del

40、ta_w+1:wp2/delta_w+1); %actual passband rippleas=-round(max(db(ws2/delta_w+1:501); %min stopband attenuationfigure(1);subplot(2,2,1);stem(n,hd);subplot(2,2,2);stem(n,w_bla);subplot(2,2,3);plot(w/pi,pha);% stem(n,h);subplot(2,2,4);plot(w/pi,db);2、利用频率取样法设计fir滤波器（采用频域抽样方法设计具有下列指标的高通滤波器： 画出滤波器的频率响应（db）

41、、理想和实际的冲激响应。实验参考例程：%=fir by sampling in frequecy clear n=33;alpha=(n-1)/2;l=0:n-1;wl=(2*pi/n)*l;t1=0.1095;t2=0.5980;hrs=zeros(1,11),t1,t2,ones(1,8),t2,t1,zeros(1,10);hdr=0,0,1,1;wdl=0,0.6,0.8,1;k1=0:floor(n-1)/2);k2=floor(n-1)/2)+1:n-1;angh=-alpha*(2*pi)/n*k1,alpha*(2*pi)/n*(n-k2);h=hrs.*exp(j*angh)

42、;h=real(ifft(h,n);db,mag,pha,grd,w=freqz_m(h,1);hr,ww,a,l=hr_type1(h);%plot%figure(1)subplot(2,2,1);plot(wl(1:17)/pi,hrs(1:17),o,wdl,hdr);axis(0,1,-0.1,1.1);title(highpass:n=33,t1=0.1095,t2=0.598)ylabel(hr(k)set(gca,xtickmode,manual,xtick,0;.6;.8;1)set(gca,xticklabelmode,manual,xticklabels,0; .6; 0.

43、8; 1)set(gca,xtickmode,manual,ytick,0;0.109;0.59;1) ;gridsubplot(2,2,2);stem(l,h);axis(-1,m,-0.4,0.4)title(impulseresponse);xlabel(n);ylabel(h(n);subplot(2,2,3);plot(ww/pi,hr,wl(1:17)/pi,hrs(1:17),o);axis(0,1,-0.1,1.1);title(amplitude response)xlabel(frequency in pi units);ylabel(hr(w)set(gca,xtickm

44、ode,manual,xtick,0;.6;.8;1)set(gca,xticklabelmode,manual,xticklabels,0;.6;.8;1)set(gca,ytickmode,manual,ytick,0;0.109;0.59;1);grid subplot(2,2,4);plot(w/pi,db);axis(0,1,-100,10);gridtitle(magnitude response);xlabel(frequency in pi units);ylabel(decibels);set(gca,xtickmode,manual,xtick,0;.6;.8;1)set(

45、gca,xticklabelmode,manual,ytick,0; .6; .8; 1)set(gca,yticklabelmode,manual,yticklabels,0;.6;.8;1)各函数文件：function hr,w,a,l=hr_type1(h);%computes amplitude response of type-1 lp fir filter%=%hr=amplitude response%w=500 frequencies between 0 pi over which hr is computed%a=type-1 lp filter coefficients%l

46、=order of he% h=type-1 lp filter impulse response% m=length(h); l=(m-1)/2; a=h(l+1) 2*h(l:-1:1); n=0:l; w=0:500*pi/500; hr=cos(w*n)*a function hr,w,b,l=hr_type2(h);%computes amplitude response of type-2 lp fir filter%=%hr=amplitude response%w=500 frequencies between 0 pi over which hr is computed%b=

47、type-2 lp filter coefficients%l=order of hr% h=type-2 lp filter impulse response% n=length(h); l=n/2; b= 2*h(l:-1:1); n=1:l;n=n-0.5 w=0:500*pi/500; hr=cos(w*n)*b function hr,w,c,l=hr_type3(h);%computes amplitude response of type-3 lp fir filter%=%hr=amplitude response%w=frequencies between 0 pi over

48、 which hr is computed%c=type-3 lp filter coefficients%l=order of hr% h=type-3 lp filter impulse response% n=length(h); l=(n-1)/2 b=2*h(l+1:-1:1); n=0:l; w=0:500*pi/500; hr=cos(w*n)*c;function hr,w,c,l=hr_type4(h);%computes amplitude response of type-4 lp fir filter%=%hr=amplitude response%w=500 frequencies between 0 pi over which hr is computed%c=type-4 lp filter coefficients%l=order of hr% h=type-4 lp filter impulse response% n=length(h); l=(n)/2 d=2*h(l:-1:1); n=1:l;n=n-0.5; w=0:500*pi/