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

《华北电力大学Matla实验指导书.doc》由会员分享，可在线阅读，更多相关《华北电力大学Matla实验指导书.doc（57页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

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命令窗口，是键入指令的地方也

3、是MATLAB显示计算结果的地方。在符号之后键入“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)

5、 c=pi*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函数impseq.m，用来生成迟延的单位脉冲序列。其输入参数为：序列起始

8、位置ns，序列终止位置nf，脉冲位置np。练习二：编写MATlAB程序，以产生下图所示的波形。并将序列绘制出来。 练习三：设x2 3 4 5，求将它延拓6个周期所得的序列。练习四：设x2 3 4 5，实现以下运算并绘制波形练习五： 编写程序实现下列函数，并绘制波形以 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

9、(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: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(

10、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 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,zer

11、os(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); %延时 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

12、.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); 得到如图23的曲线图23此曲线与用rec函数所画曲线完全一致, 可见, 用MATLAB编制的FIR滤波程序与理论计算所得结果是相同的.3、用nonrec.m和rec.m函数编制dfilter.m 函数文件, 构成完整的一般IIR滤波程序, 并完成下列信号

13、滤波: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)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, 并与你编写的dfi

14、lter.m结果进行比较.用help可以看到内部函数为Y = FILTER(B,A,X), 且有 “The filter is a Direct Form II Transposed implementation of the standard 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

15、: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);图25三、 思考题1、 内部函数filter.m的调用格式是什么?与你编写的dfilter.m调用格式是否一致?差异在何处?2、 实验中如何合理地选取单位抽样响应h(n)的点数?实验内容及要求中第三项的物理概念是什么?给定的H(z)是

16、具有什么性质的滤波器?3、 为实现各实验内容中的滤波器，你编写程序时采用的什么结构？4、 试利用Matlab中的卷积函数conv完成实验内容1,并与本书的例程进行比较。实验三 离散傅里叶变换（DFT）及其快速算法（FFT）一、实验目的1、学习DFT和FFT的原理及编程实现方法。2、熟悉MATLAB编程的特点。二、实验内容1、用三种不同的DFT程序计算x(n)=R8(n)的傅里叶变换X (ejw)，并比较三种程序计算机运行时间。（1）用for loop 语句的M函数文件dft1.m，用循环变量逐点计算X(k)；（2）编写用MATLAB矩阵运算的M函数文件dft2.m, 完成下列矩阵运算；（3）调

17、用FFT库函数，直接计算X(k)；（4）分别利用上述三种不同方式编写的DFT程序计算序列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);

18、n=0:N-1;k=0:N-1;w=exp(-j*2*pi/N);nk=n*k;wnk=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,

19、b); grid;title(Magnitude part);xlabel(frequency 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=（2）x=ones(1,8),zeros(1,248);t=cputime;Am2,pha2=dft2(x);t2=cputime-tn=0:(length(x)-1);w=(2*

20、pi/length(x)*n;figure(2)subplot(2,1,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=（3）x=ones(1,8),zeros(1,248);t=cputime;

21、Am3,pha3=dft3(x);t3=cputime-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

22、DFT2记录运行时间：t3=从以上运行结果可以看出，调用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*

23、t);Xk2=fft(x2,128); %做128点傅里叶变换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

24、;w=(2*pi/128)*n;figure(3)plot(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)的关系，并说

25、明为什么抽取后能分辨出两个频率分量？有无任何条件的限制？实验四 FFT的典型应用一、实验目的1、学习FFT的应用于卷积运算和谱分析时的编程实现方法。2、熟悉MATLAB编程的特点。二、实验内容1、利用FFT实现两序列的卷积运算，并研究FFT点数与混叠的关系。给定x(n)=R16(n)，h(n)=R8(n)，用FFT和IFFT分别求线性卷积和混叠结果输出，并用函数stem(n,y)画出相应图形。（一）实验例程源程序如下：N=16;x=0:N-1;h=ones(1,8);线性卷积：Xk1=fft(x,23); %做23点fftHk1=fft(h,23);Yk1=Xk1.*Hk1;y1=ifft(Y

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

27、性。（2） 采集数据长度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);Xk1=fft (xa);Am1=abs(Xk1);n=0:(length(xa)-1);w=(2*pi/length(xa)*n;fig

28、ure(1)plot(w,Am1,b); grid;title(Magnitude part);xlabel(frequency in radians);ylabel(|X(exp(jw)|);x=xa(1:16),zeros(1,240); %补零到256Xk2=fft(x); %做256点DFTAm2=abs(Xk2);n=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)|)

29、;N0=256; %采集256点n0=0: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);Xk3=fft(xa0); Am3=abs(Xk3);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

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

31、)的FFT变换X(k);(2)将（1）中的 x(n)以补零方式加长到长度为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(1,8); %h为单位抽样响应H=fft(h); for i=0:14 if

32、 (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),则可得到如下曲线. 与重叠保留法所得结果相同.图36 直接卷积结果三、思考题1、什么是高密度频谱

33、？什么是高分辨率频谱？实验内容1中为区分三个频率分量至少应采集几个样本？2、实验内容2中，分别求X1(ejww)和x2(ejww)与原X(ejww)的关系，并说明为什么抽取后能分辨出两个频率分量？有无任何条件的限制？实验五 FIR滤波器设计一、实验目的1、学习FIR滤波器的设计方法。2、熟悉MATLAB编程的特点。二、实验内容1、利用窗函数法设计FIR滤波器 设计具有下列指标的低通数字滤波器（线性相位）选择合适的窗函数，并画出滤波器的频率响应（dB）、窗函数、理想和实际的冲激响应。实验参考例程：%=FIR by windowsclear;wp=0.2*pi;ws=0.3*pi;tr_width

34、=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);figure(1);subplot(2,2,1);stem(n,hd);subplot(2,2,2);stem(n,w_ham);subplot(2

35、,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);function db,mag,pha,grd,w=freqz_m(b,a)H,w=freqz(b,a,1000,whole);H=(H(1:501

36、);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（）函数实现题目要求，并和例程比较运行结果。练习四：分析并运行下列程序，说明该程序实现的是何种功能的滤波器？%=FIR by windowsclear;wp1=0.35*pi;ws1=0.2*pi;wp1=0.65*pi;ws1=0.8*pi;As=60;tr_width=min(wp1-ws1),(ws2

37、-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/delta_w+1:wp2/delta_w+1); %Actual passband RippleAs=-round(max(db(ws2/delta_w+1:501); %Min Stopband

38、 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）、理想和实际的冲激响应。实验参考例程：%=FIR by sampling in frequecy Clear N=33;alpha=(N-1)/2;L=0:N-1;wL=(2*pi/N)*L;T

39、1=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);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,

40、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.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,

41、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,XTickMode,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