南京郵電大學(xué)實(shí)驗(yàn)報(bào)告實(shí)驗(yàn)名稱：采樣、系統(tǒng)性質(zhì)及濾波系統(tǒng)頻率響應(yīng)和樣本處理算法實(shí)現(xiàn)加窗和離散傅氏變換數(shù)字濾波器設(shè)計(jì)課程名稱:數(shù)字信號(hào)處理姓名：學(xué)號(hào)開(kāi)課時(shí)間2011/2012學(xué)年，第2學(xué)期實(shí)驗(yàn)一．實(shí)驗(yàn)名稱：采樣、系統(tǒng)性質(zhì)及濾波二．實(shí)驗(yàn)?zāi)康暮腿蝿?wù)，實(shí)驗(yàn)內(nèi)容：一、觀察采樣引起的混疊。設(shè)模擬信號(hào)為x(t)=cos(5兀t)+4sin(2兀t)-sin(3兀t),t的單位為毫秒(ms)。設(shè)采樣頻率為3kHz，確定與x(t)混疊的采樣重建信號(hào)x(t)。a畫出x(t)和x(t)在0<t<6(ms)范圍內(nèi)的連續(xù)波形。(因數(shù)字計(jì)算機(jī)無(wú)法真正畫出a連續(xù)波形，可用較密的離散點(diǎn)的連線來(lái)近似。)分別用""和"x"在兩信號(hào)波形上標(biāo)記出3kHz采樣點(diǎn)。兩信號(hào)波形是否相同？采樣后的兩序列是否相同？二、 判別離散時(shí)間系統(tǒng)的時(shí)不變性。(來(lái)源：p105例3.2.2)設(shè)輸入序列為x(n)，系統(tǒng)y(n)=x(2n)實(shí)現(xiàn)對(duì)x(n)的抽取。2兀設(shè)x(n)=sin(而n),n=1,2,...,500。取延遲量D(例如D=30)。記x(n)=x(n-D)，畫出x(n)、x(n)的序列波形。DD編程求出系統(tǒng)對(duì)x(n)的響應(yīng)y(n)以及對(duì)x(n)的響應(yīng)y(n)DD畫出y(n-D)、y(n)的波形。D該系統(tǒng)是否為時(shí)不變的？三、 利用卷積計(jì)算信號(hào)通過(guò)FIR濾波器的輸出，并觀察輸出信號(hào)的input-on暫態(tài)、input-of暫態(tài)和穩(wěn)態(tài)階段。(來(lái)源：p144例4.1.8)考慮兩個(gè)濾波器，h1(n)=] 0 其宀，h2=5[1,—5,0,-10‘5'-1];[10<x<24輸入x(n)為周期方波，第一個(gè)周期內(nèi)x(n)= ”，/“。[0 25<x<491..分別畫出x(n)通過(guò)兩個(gè)濾波器的輸出y(n)、y(n)的波形(0<n<199)，并與書12上p144例4.1.8的兩幅圖比較是否一致。2.計(jì)算圖中穩(wěn)態(tài)部分的響應(yīng)值。三．實(shí)驗(yàn)過(guò)程與結(jié)果1、實(shí)驗(yàn)程序：%=============%problem1%=============clear%estimatex(t)andxa(t)withamuchhighersamplingfreq.'fs1'time_period=6;%unit:msfs1=50;%unit:kHzT1=1/fs1;%unit:msn1=0:fix(time_period/T1);x=cos(5*pi*n1*T1)+4*sin(2*pi*n1*T1).*sin(3*pi*n1*T1);xa=cos(pi*n1*T1);%obtainx(nT)andxa(nT)withgivensamplingfreq.'fs'fs=3;T=1/fs;n=0:fix(time_period/T);x_sample=cos(5*pi*n*T)+4*sin(2*pi*n*T).*sin(3*pi*n*T);xa_sample=cos(pi*n*T);figure,plot(n1*T1,x,'r',n1*T1,xa,'b',n*T,x_sample,'ro'),holdon,stem(n*T,xa_sample,'b:x')legend('x(t)','xa(t)','x(nT)','xa(nT)'),xlabel('t(ms)')%=============%problem2%=============clear%plotx(n)andx(n-D)D=30;N=500;n=1:N;x=sin(2*pi/100*n);forn=1:N+D,if(n-D)<=0,xD(n)=0;elsexD(n)=x(n-D);endendfigure,subplot(2,1,1),plot(1:N,x,'r:',1:length(xD),xD,'b'),legend('x(n)','xD(n)'),,xlabel('n')%ploty(n)andyD(n)forn=1:fix(N/2)y(n)=x(2*n);endforn=1:length(y)+D,if(n-D)<=0,y_delay(n)=0;elsey_delay(n)=y(n-D);endendforn=1:fix(length(xD)/2)yD(n)=xD(2*n);endsubplot(2,1,2),plot(1:length(y),y,'r:',1:length(y_delay),y_delay,'r.:',1:length(yD),yD,'b.'),legend('y(n)','y(n-D)','yD(n)'),,xlabel('n')axis([0530-11])%=============%problem3%=============clearhl=0.25*0.75.人(0:14);h2=1/5*[1-510-105-1];N=200;n=0:N-1;x1=[ones(1,25)zeros(1,25)];%oneperiodof'x(n)'x=[x1x1x1x1];y1=conv(x,h1);y2=conv(x,h2);figure,subplot(2,1,1),,plot(n,x,'r:',n,y1(1:N),'b'),axis([0200-0.52.5]),gridon,legend('input','output'),xlabel('n')subplot(2,1,2),plot(n,x,'r:',n,y2(1:N),'b'),axis([0200-1.52.5]),gridon,legend('input','output'),xlabel('n')[4恥[4恥「51覓(3)寸RQUitIEjI-i囲LYiirlustiI<-^3Quk^ipILtinJilp□S?R凈型嗖電題圮］_P廚I□口四?實(shí)驗(yàn)分析a).兩信號(hào)x(t)和x(t)波形是否相同？采樣后的兩序列x(nT)和x(nT)是否相同？a a反映了什么現(xiàn)象？答：兩信號(hào)波形不相同。采樣的兩序列相同。反映了采樣重構(gòu)時(shí)發(fā)生混淆現(xiàn)象。.參考程序第一段的語(yǔ)句：x=cos(5*pi*nl*Tl)+4*sin(2*pi*nl*Tl).*sin(3*pi*nl*Tl)中用計(jì)算符“*”代替“?*”，結(jié)果如何？如果進(jìn)一步將參與計(jì)算的兩正弦矢量sin(2*pi*nl*Tl)和sin(3*pi*nl*Tl)分別進(jìn)行轉(zhuǎn)置(提示：矢量y的轉(zhuǎn)置為y',再進(jìn)行“*”，結(jié)果又如何？答：前者程序報(bào)錯(cuò)。后者計(jì)算正確。.改用綠色畫出x(t)。a答：將程序中的figure,plot(nl*Tl,x,'r',nl*Tl,xa,'b',n*T,x_sample,'ro')改為：figure,plot(nl*Tl,x,'r',nl*Tl,xa,'g',n*T,x_sample,'ro'),a)?該系統(tǒng)是否為時(shí)不變的？答：不是.b).自上而下用四個(gè)小窗口分別畫出x(n)和y(n),y(nD),y(n)的波形D0(3).a).兩個(gè)濾波器的DCgain分別是多少？答：第一個(gè)：DCgain=l-(0.75)15沁0.987第二個(gè)：DCgain=0b).響應(yīng)的input-on暫態(tài)、input-off暫態(tài)和穩(wěn)態(tài)段范圍分別是多少？答：第一個(gè)周期：input-on暫態(tài)：0Wn<14input-off暫態(tài)：49<n<63穩(wěn)態(tài):14<n<49(4)實(shí)驗(yàn)心得通過(guò)實(shí)驗(yàn)一，我鞏固了采樣重構(gòu)、系統(tǒng)性質(zhì)和濾波這些章節(jié)所學(xué)習(xí)的知識(shí)。在實(shí)驗(yàn)中，了解并加強(qiáng)了對(duì)Matlab的學(xué)習(xí)，提高了編程能力。同時(shí)，從Matlab編程所得程序的結(jié)果中，更加直觀的反映了前幾章節(jié)的重點(diǎn)，理解更加深刻。實(shí)驗(yàn)二一.實(shí)驗(yàn)名稱：系統(tǒng)頻率響應(yīng)和樣本處理算法實(shí)現(xiàn)二實(shí)驗(yàn)?zāi)康暮腿蝿?wù)，實(shí)驗(yàn)內(nèi)容：一、觀察序列頻譜，觀察信號(hào)通過(guò)系統(tǒng)后波形與頻譜的變化兀n 兀n已知輸入信號(hào)x(n)=5x(n)+2x(n)，其中x(n)=cos()，x(n)=cos()，121824n=0,…N-1,N可取5000點(diǎn)。畫出x(n),x(n),x(n)的前100點(diǎn)波形12畫出x(n)的DTFT頻譜X(w)(0 ＜兀)由于計(jì)算機(jī)無(wú)法畫出連續(xù)頻譜，所以可在0＜w＜兀內(nèi)均勻取足夠密的點(diǎn)數(shù)，如兀M=5000個(gè)頻率點(diǎn)w=k(k=0,…M-1),求出這些頻率點(diǎn)上的頻譜值kMX(w)=￡'x(n)e-wk(k=0,…M-1)，并畫出IX(w)I隨w變化的曲線。TOC\o"1-5"\h\zk k kn=0某LTI系統(tǒng)h(n)=u(n)-u(n-8)，畫出系統(tǒng)的幅度頻響IH(w)I求系統(tǒng)對(duì)x(n)的響應(yīng)y(n)(可以自己編程也可利用卷積函數(shù)conv(-))畫出y(n)的\o"CurrentDocument"兀 兀 兀 兀 兀 兀波形，并與5IH(石)lcos(￡n+argH(石))+21H(丁)Icos(n+argH(丁))的波形比較\o"CurrentDocument"8 8 8 4 4 4(各畫100點(diǎn))；畫出y(n)的幅度譜IY(w)|，并與IX(w)-H(w)1比較。你從中觀察到什么？二、系統(tǒng)函數(shù)H(z)= ，根據(jù)正準(zhǔn)型結(jié)構(gòu)(canonicalform)編寫樣本處理算法。1-0.8z-1內(nèi)部狀態(tài)的初始值設(shè)為零，輸入信號(hào)x(n)采用逐個(gè)樣本手動(dòng)輸入的方式，求輸出信號(hào)y(n)。三．實(shí)驗(yàn)過(guò)程與結(jié)果1.實(shí)驗(yàn)代碼：% %Programforthefirstproblem% colordefnoneN=5000; %setlengthofthetimeseriesx1=cos(pi/8*(0:N-1));x2=cos(pi/4*(0:N-1));x=5*x1+2*x2;m=100; figure,plot(1:m,x1(1:m),1:m,x2(1:m),'g',1:m,x(1:m),'b.-'),legend('x1','x2','x')%---plotX(w)M=5000; %setsamplesnumberinfrequencydomainw=pi/M*(0:M-1);%discretedfrequencyj=sqrt(-1);Xw=zeros(1,M);fork=1:M,Xw(k)=sum(x*(exp(-j*w(k)*(0:N-1)')));endfigure,plot(w/pi,abs(Xw)),xlabel('\omega(\pi)'),ylabel('|X(\omega)|')%---plotH(w)h=ones(1,8);fork=1:M,Hw(k)=sum(h*(exp(-j*w(k)*(0:7)')));endfigure,plot(w/pi,abs(Hw)),xlabel('\omega(\pi)'),ylabel('|H(\omega)|')%---plotY(w)y=conv(x,h);%calculatetheresponseH1=sum(h*(exp(-j*pi/8*(0:7)')));H2=sum(h*(exp(-j*pi/4*(0:7)')));m=100;yy=5*abs(H1)*cos(pi/8*(0:m-1)+angle(H1))+2*abs(H2)*cos(pi/4*(0:m-1)+angle(H2));figure,plot(y(1:m),'r.-'),holdon,plot(yy,'g'),legend('y','yy')Ny=length(y);fork=1:M,Yw(k)=sum(y*(exp(-j*w(k)*(0:Ny-1)')));endfigure,plot(w/pi,abs(Yw),'r.-'),holdon,plot(w/pi,abs(Xw.*Hw)),xlabel('\omega(\pi)'),legend('|Y(\omega)|','|X(\omega)|*|H(\omega)|')% %Programforthesecondproblem% w1=0;fori=1:10,x=input('inputx=');w0=0.8*w1+x;y=w0-0.5*w1;w1=w0;yend2.運(yùn)行結(jié)果(1)

EalvE4iLbin■-IewiHTwh lijulnIaL. w□ it'i魚鬲4□目口(2)inputx=2y=2inputx=4y=4.6000inputx=5y=6.6800inputx=6y=8.8440inputx=7y=11.0752inputx=8y=13.3602inputx=8y=14.6881inputx=6y=13.7505inputx=4y=12.0004inputx=2y=9.6003四?實(shí)驗(yàn)分析(1).a).濾波器h(n)是什么頻響類型的濾波器?答：低通b).你從以上實(shí)驗(yàn)中觀察到什么？與課本Chapter6的什么重要結(jié)論相吻合？希望同學(xué)們通過(guò)此實(shí)驗(yàn)加深對(duì)理論教學(xué)內(nèi)容的理解。1.Y(z)=H(z)X(z)2.對(duì)正弦信號(hào)：Hcos(?n)T|H)|cos(?n+argH))0000sin(?n)T|H)1sin(?n+argH))0000c).觀察實(shí)際仿真中的現(xiàn)象，指出有限長(zhǎng)正弦輸入引起的input-ontransients,input-offtransients響應(yīng)段。input-ontransients段：0~8實(shí)驗(yàn)三一．實(shí)驗(yàn)名稱：加窗和離散傅氏變換二．實(shí)驗(yàn)?zāi)康暮腿蝿?wù)，實(shí)驗(yàn)內(nèi)容：1．觀察窗函數(shù)的影響。信號(hào)為x(t)二cos(2吋t)+cos(2吋t)+cos(2吋t)，f=2KHz，f=2.5KHz，f=31 2 3 1 2 3KHz，采樣頻率f二10KHz。sa)寫出x(n)(—g<n<+a)的頻譜X(①)；分別畫出窗長(zhǎng)度L=10，L=20，L=40，L=100的矩形窗頻譜和Hamming窗頻譜。觀察L的變化對(duì)窗譜的主瓣寬度、旁瓣密集度、相對(duì)旁瓣水平的影響。時(shí)域采樣點(diǎn)數(shù)分別取L=10，L=20，L=40，L=100，分別畫出x(n)加矩形窗及加Hamming窗時(shí)DTFT頻譜X(①)；你能否從頻譜上分辨出信號(hào)的三個(gè)頻率分量？若能分辨出，它們的位置和相對(duì)大小是否準(zhǔn)確？L的大小對(duì)頻率的物理分辨率(physicalfrequencyresolution)有何影響？2．理解頻率的物理分辨率和計(jì)算分辨率的區(qū)別信號(hào)同上，加矩形窗。時(shí)域采樣點(diǎn)數(shù)分別取L=10，L=20，L=40，L=100。畫出以上各種時(shí)長(zhǎng)情況下，頻域采樣點(diǎn)數(shù)分別為N=32，N=64時(shí)的DFT(在同一個(gè)圖上用虛線畫出相應(yīng)的DTFT頻譜，用于比較)。離散頻譜DFT和連續(xù)頻譜DTFT有什么關(guān)系？L一定的情況下，能否通過(guò)增加N改善頻率的物理分辨率？N的作用是什么？附：參考程序（注：由一個(gè)主程序和三個(gè)函數(shù)程序組成。其中函數(shù)程序名要與函數(shù)名相同）三．實(shí)驗(yàn)過(guò)程與結(jié)果1.實(shí)驗(yàn)代碼：%================%Problem1%================% Plotwaveformandspectrumofthewindowsfigure,L=10;windows_spectrum(L,[221222])L=20;windows_spectrum(L,[223224])figure,L=40;windows_spectrum(L,[221222])L=100;windows_spectrum(L,[223224])% plotthesignal'sspectrum L=100;x=cos(0.4*pi*(0:(L-1)))+cos(0.5*pi*(0:(L-1)))+cos(0.6*pi*(0:(L-1)));L=10;figure,windowed_spectrum(x(1:L),[221222])L=20;windowed_spectrum(x(1:L),[223224])L=40;figure,windowed_spectrum(x(1:L),[221222])L=100;windowed_spectrum(x(1:L),[223224])%================%Problem2%================L=10;figure,windowed_dft_spectrum(x(1:L),[221222])L=20;windowed_dft_spectrum(x(1:L),[223224])L=40;figure,windowed_dft_spectrum(x(1:L),[221222])L=100;windowed_dft_spectrum(x(1:L),[223224])子程序1:function[]=windowed_dft_spectrum(x,a);L=length(x);W1=abs(fft(x,1024));omega=(0:1023)*2/1024;N=32;kk=2/N*(0:N-1);ifN>=L,dft=abs(fft(x,N));elsedft=zeros(N,1);fork=1:N,dft(k)=sum(x(:)'*exp(-j*2*pi/N*(k-1)*(0:L-1)'));endendsubplot(a(1),a(2),a(3)),plot(omega,W1,':'),holdon,stem(kk,abs(dft),'r.'),xlabel('\omega(\pi)'),axis([02050]),title(['L=',num2str(L),',N=',num2str(N)])N=64;kk=2/N*(0:N-1);ifN>=L,dft=abs(fft(x,N));elsedft=zeros(N,1);fork=1:N,dft(k)=sum(x(:)'*exp(-j*2*pi/N*(k-1)*(0:L-1)'));endendsubplot(a(4),a(5),a(6)),plot(omega,W1,':'),holdon,stem(kk,abs(dft),'r.')xlabel('\omega(\pi)'),axis([02050]),title(['L=',num2str(L),',N=',num2str(N)])子程序2：functionwindowed_spectrum(x,a)L=length(x);w=hamming(L);W1=abs(fft(x,1024));W2=abs(fft(x.*w',1024));omega=(0:1023)*2/1024;subplot(a(1),a(2),a(3)),plot(omega,W1),xlabel('\omega(\pi)'),axis([02050])title(['X(w)withrectangularwindow,L=',num2str(L)]),gridonsubplot(a(4),a(5),a(6)),plot(omega,W2),xlabel('\omega(\pi)'),axis([02050])title(['X(w)withHammingwindow,L=',num2str(L)]),gridon子程序3：functionwindows_spectrum(L,a)rect=rectwin(L);%rectangularwindowhamm=hamming(L);%Hammingwindoww=2*pi/1024*(-511:512);%discretedfrequencyW1=fft(rect,1024);W1=abs(fftshift(W1));W2=fft(hamm,1024);W2=abs(fftshift(W2));subplot(a(1),a(2),a(3)),holdon,plot(1:L,rect,'bo',1:L,hamm,'r*'),stem(1:L,rect),axis([1L-0.52]),xlabel('n'),title(['L=',num2str(L),'(waveform)']),legend('rectangular','Hamming')subplot(a(4),a(5),a(6)),plot(w/pi,abs(W1),w/pi,abs(W2),'r'),xlabel('\omega(\pi)'),title(['L=',num2str(L),'(Magnitudespectrum)'])legend('rectangular','Hamming')2.實(shí)驗(yàn)結(jié)果：i*斤QUl亡］SCLn3*詁*1YtwIrw+r: gktrfIninliLfn￡9^:;*＜丸倉(cāng)％口廚匚可L-1Qjrarfonrij L-1QIMagnUdespedniri|^0IKiJngjjUr*Hinrnrq■ ■ ■ ■門■8A Wt^ngULiir HinrnnqL=2SI

Ln3*詁*1KtwIrw+r: gktrfLi-ainliLf2曰日If□nL=1DFJ=JJLn3*詁*1KtwIrw+r: gktrfLi-ainliLf2曰日If□nL=1DFJ=JJL=1D5Q ! ! ! ! 4J]3]-2i]L如曰？四?實(shí)驗(yàn)分析(1)a).觀察窗長(zhǎng)L的變化對(duì)窗函數(shù)頻譜的主瓣寬度、旁瓣密集度、相對(duì)旁瓣水平的影響。答：L越長(zhǎng)，主瓣寬度越窄、旁瓣寬度越密集、相對(duì)旁瓣水平影響越小。.你能否從信號(hào)頻譜上分辨出信號(hào)的三個(gè)頻率分量？若能分辨出，它們的位置和相對(duì)大小是否準(zhǔn)確？能，不準(zhǔn)確.矩形窗與漢明窗矩形窗有旁瓣，主瓣寬度是同等條件下漢明窗的二分之一，主瓣高度高于漢明窗。A 2兀A?=c ，其中，矩形窗c取1,漢明窗c取2。wL(2)a).離散頻譜DFT和連續(xù)頻譜DTFT有什么關(guān)系？DTFT定義：X(?)=藝x(n)e-=藝x(n)e-LLn=0 n二一8X(?)=Ex(n)e-j?n(DTFToflenth-Lsignal)n=0DFT(N一pointDFToflength一Lsignal)的定義：X(?)=X(k)=Ex(n)e-j?k=Ex(n)e-jnn(9.2.8)kn=0 n=0(即相當(dāng)于在Z平面單位圓上取N等分點(diǎn)，在這些等分點(diǎn)上計(jì)算DTFT, Fig.9.2.3)N-pointDFTX(k)的周期為：Nb).本實(shí)驗(yàn)實(shí)際上是對(duì)確定性信號(hào)序列的譜分析。須通過(guò)實(shí)驗(yàn)加深理解time-windowing和frequencysampling對(duì)譜分析結(jié)果的影響，而譜分析分辨率可以用physical/computationalfrequencyresolution來(lái)衡量。回答以下問(wèn)題：L的大小對(duì)頻率的物理分辨率(physicalfrequencyresolution)有何影響？L一定的情況下，能否通過(guò)增加N改善頻率的物理分辨率？N的作用是什么？答：L越大，物理分辨率越高L一定時(shí)，不能通過(guò)增加N改善物理分辨率，N的作用是改善計(jì)算分辨率。實(shí)驗(yàn)四一．實(shí)驗(yàn)名稱：數(shù)字濾波器設(shè)計(jì)二．實(shí)驗(yàn)?zāi)康暮腿蝿?wù)，實(shí)驗(yàn)內(nèi)容1.窗口法設(shè)計(jì)FIR數(shù)字濾波器(來(lái)源：p571習(xí)題10.18)(a)用矩形窗設(shè)計(jì)長(zhǎng)度分別為N=11、41、81、121的低通FIR濾波器，要求截止頻率為①二0.3兀。畫出濾波器的單位沖激響應(yīng)h(n)和幅度頻響IH(w)1曲線。c問(wèn)題：理想濾波器的頻響是怎樣的？當(dāng)N增大時(shí)，F(xiàn)IR濾波器在w附近的最大紋波幅度是否降低？其余紋波的幅度是否減c?。考y波的密度怎樣變化？過(guò)渡帶寬度怎樣變化？(如有必要可增大N值觀察)。在N=ll時(shí)，畫出濾波器的相頻曲線。它是否是線性的？(b)用漢明窗再次設(shè)計(jì)同樣的濾波器。用漢明窗設(shè)計(jì)出的濾波器與用矩形窗相比有什么特點(diǎn)？2．以Butterworth模擬低通濾波器為原型，設(shè)計(jì)IIR數(shù)字濾波器。截止頻率w二0?3兀。設(shè)計(jì)11階IIR數(shù)字低通濾波器，畫出幅頻、相頻曲線。c截止頻率w二0?3兀。設(shè)計(jì)11階IIR數(shù)字高通濾波器，畫出幅頻、相頻曲線。c問(wèn)題：所設(shè)計(jì)的IIR濾波器與FIR濾波器的頻率特性有何區(qū)別？三．實(shí)驗(yàn)過(guò)程與結(jié)果1.實(shí)驗(yàn)代碼：%================================%DesignFIRDFwithwindowmethod%===============================wc=0.3*pi;% rectangularwindowwin=1;p=1; %observemagnituderesponseN=11;figure,FIR_design(N,wc,win,[221222],p)N=41;FIR_design(N,wc,win,[223224],p)N=81;figure,FIR_design(N,wc,win,[221222],p)N=121;FIR_design(N,wc,win,[223224],p)p=2; %observephaseresponseN=11;figure,FIR_design(N,wc,win,[111111],p)% hammingwindowwin=2;p=1;N=11;figure,FIR_design(N,wc,win,[221222],p)N=41;FIR_design(N,wc,win,[223224],p)N=81;figure,FIR_design(N,wc,win,[221222],p)N=121;FIR_design(N,wc,win,[223224],p)%=================================%DesignIIRDF%==================================N=11;Wn=wc/pi;[b,a]=BUTTER(N,Wn); %lowpass[H,w]=freqz(b,a,8192);figure,subplot(2,1,1),plot(w/pi,abs(H)),xlabel('\omega(\pi)'),ylabel('|H(\omega)|'),title(['order=',num2str(N),',LowpassIIRDF:magnituderesponse']),gridsubplot(2,1,2),plot(w/pi,angle(H)/pi),xlabel('\omega(\pi)'),ylabel('ang(H(\omega))(\pi)'),title('phaseresponse')N=11;[b,a]=BUTTER(N,Wn,'high');%highpass[H,w]=freqz(b,a,8192);figure,subplot(2,1,1),plot(w/pi,abs(H)),xlabel('\omega(\pi)'),ylabel('|H(\omega)|'),title(['order=',num2str(N),',HIghpassIIRDF:magnituderesponse']),gridsubplot(2,1,2),plot(w/pi,angle(H)/pi),xlabel('\omega(\pi)'),ylabel('ang(H(\omega))(\pi)'),title('phaseresponse')子程序1:functionFIR_design(N,wc,win,a,p)%LPFIRdesignwithwindowmethodn=0:N-1;M=(N-1)/2;warningoffMATLAB:divideByZeroh=sin(wc*(n-M))./(pi*(n-M));h(M+1)=wc/pi;ifwin==2,h=h(:).*hamming(N);end[Hw]=freqz(h,1,16400);ifp==1,subplot(a(1),a(2),a(3)),stem(0:N-1,h,'.'),xlabel('n'),ylabel('h(n)'),title(['N=',num2str(N),',h(n)'])subplot(a(4),a(5),a(6)),plot(w/pi,abs(H)),xlabel('\omega(\pi)'),ylabel('|H(\omega)|'),title(['N=',num2str(N),',|H(\omega)|']),gridelsesubplot(a(1),a(2),a(3)),plot(w/pi,angle(H)/pi,w/pi,abs(H),'r:'),legend('ang(H(\omega)','|H(\omega)|')xlabel('\omega(\pi)'),ylabel('ang(H(\omega