實(shí)驗(yàn)1 利用DFT分析信號頻譜

實(shí)驗(yàn)1利用DFT分析信號頻譜一、實(shí)驗(yàn)?zāi)康蘑偌由顚FT原理的理解②應(yīng)用DFT分析信號的頻譜③深刻理解利用DFT分析信號頻譜的原理，分析實(shí)現(xiàn)過程中出現(xiàn)的現(xiàn)象及解決方法二、實(shí)驗(yàn)設(shè)備與環(huán)境計(jì)算機(jī)、MATLAB軟件環(huán)境。三、實(shí)驗(yàn)基礎(chǔ)理論1.DFT與DTFT的關(guān)系有限長序列的離散時(shí)間傅里葉變換在頻率區(qū)間的N個(gè)等間隔分布的點(diǎn)上的N個(gè)取樣值可以由下式表示：由上式可知，序列的N點(diǎn)DFT，實(shí)際上就是序列的DTFT在N個(gè)等間隔頻率點(diǎn)上的樣本。2.利用DFT求DTFT方法1：由恢復(fù)出的方法如圖2.1所示：DTFTIDFTDTFTIDFT圖2.1由N點(diǎn)DFT恢復(fù)頻譜DTFT的流程由圖2.1所示流程可知：其中為內(nèi)插函數(shù)方法2：然而在MATALAB計(jì)算中，上述插值運(yùn)算不見得是最好的辦法。由于DFT是DTFT的取樣值，其相鄰兩個(gè)頻率樣本點(diǎn)的間距為，所以如果我們增加數(shù)據(jù)的長度N，使得到的DFT譜線就更加精細(xì)，其包絡(luò)就越接近DTFT的結(jié)果，這樣就可以利用DFT來近似計(jì)算DTFT。如果沒有更多的數(shù)據(jù)，可以通過補(bǔ)零來增加數(shù)據(jù)長度。3.利用DFT分析連續(xù)時(shí)間信號的頻譜采用計(jì)算機(jī)分析連續(xù)時(shí)間信號的頻譜，第一步就是吧連續(xù)時(shí)間信號離散化，這里需要進(jìn)行兩個(gè)操作：采樣，截?cái)唷τ谶B續(xù)時(shí)間非周期信號，按采樣間隔T進(jìn)行采樣，截取長度為M，那么對進(jìn)行進(jìn)行N點(diǎn)頻域采樣，得到因此，可以將利用DFT分析聯(lián)系非周期信號頻譜的步驟歸納如下：（1）確定時(shí)域采樣間隔，得到離散序列；（2）確定截取長度M，得到M點(diǎn)離散序列，這里為窗函數(shù)。（3）確定頻域采樣點(diǎn)數(shù)N，要求。（4）利用FFT計(jì)算離散序列的N點(diǎn)DFT，得到。（5）根據(jù)式（2-6）由計(jì)算采樣點(diǎn)的近似值。采用上述方法計(jì)算的頻譜，需要注意如下三個(gè)問題：（1）頻譜混疊。如果不滿足采樣定理的條件，頻譜會(huì)出現(xiàn)混疊誤差。對于頻譜無限寬的信號，應(yīng)考慮覆蓋大部分主要頻率分量的范圍。（2）柵欄效應(yīng)和頻譜分辨率。使用DFT計(jì)算頻譜，得到的結(jié)果只是N個(gè)頻譜樣本值，樣本值之間的頻譜是未知的，像通過一個(gè)柵欄觀察頻譜，成“柵欄效應(yīng)”。頻譜分辨率與記錄長度成反比，要提高頻譜分辨率，就要增加記錄時(shí)間。（3）頻譜泄露。對信號截?cái)鄷?huì)把窗函數(shù)的頻譜引入信號頻譜，造成頻譜泄露。解決這個(gè)問題的主要辦法是采用瓣小的窗函數(shù)，頻譜泄露和窗函數(shù)均會(huì)引起誤差。因此，要合理選取采樣間隔和截取長度，必要時(shí)還需考慮增加適當(dāng)?shù)拇病τ谶B續(xù)時(shí)間周期信號，我們采用計(jì)算機(jī)進(jìn)行計(jì)算時(shí)，也總是要進(jìn)行截?cái)?，序列總是有限長的，仍然可以采用上述方法近似計(jì)算。4.可能用到的MATLAB函數(shù)與代碼實(shí)驗(yàn)中DFT運(yùn)算可采用MATLAB中提供的函數(shù)FFT來實(shí)現(xiàn)。DTFT可以利用MATLAB矩陣運(yùn)算的方法進(jìn)行計(jì)算四、實(shí)驗(yàn)內(nèi)容1.已知，完成以下要求：（1）計(jì)算其DTFT，并畫出區(qū)間的波形。程序代碼：n=0:3;x=[2,-1,1,1];w=-pi:0.01*pi:pi;M=x*exp(-j*n'*w);subplot(211);plot(w,abs(M));xlabel('\Omega/\pi');title('Magnitude');axistightsubplot(212);plot(w,angle(M)/pi);xlabel('\Omega/\pi');title('Magnitude');axistight程序結(jié)果：（2）計(jì)算4點(diǎn)DFT，并把結(jié)果顯示在（1）所畫的圖形中程序代碼：n=0:3;x=[2,-1,1,1];w=-pi:0.01*pi:pi;M=x*exp(-j*n'*w);subplot(211);plot(w,abs(M));xlabel('\Omega/\pi');title('Magnitude');axistightsubplot(212);plot(w,angle(M)/pi);xlabel('\Omega/\pi');title('Magnitude');axistightx=[2,-1,1,1];y=fft(x);subplot(211)holdonstem(0:3,abs(y)/pi,'fill')title('Magnitude')subplot(212)holdonstem(0:3,angle(y)/pi,'fill')ylabel('/\pi')title('Phase')（3）對補(bǔ)零，計(jì)算64點(diǎn)DFT，并顯示結(jié)果程序代碼：x=[2-111];x=[x,zeros(1,60)];y=fft(x,64);subplot(211);stem(0:63,abs(y),'fill')xlabel('\Omega/\pi')title('Magnitude')subplot(212);stem(0:63,angle(y)/pi,'fill')xlabel('\Omega/\pi')ylabel('/\pi')title('Phase')程序結(jié)果：（4）根據(jù)實(shí)驗(yàn)結(jié)果，分析是否可以由DFT計(jì)算DTFT，如果可以，如何實(shí)現(xiàn)。根據(jù)以上圖形可得我們可以通過由DFT計(jì)算DTFT?？梢栽谌〉命c(diǎn)數(shù)后面補(bǔ)零再做DFT，則可以很好的逼近DTFT的結(jié)果。2.考察序列（1）時(shí)，用DFT估計(jì)的頻譜；將補(bǔ)零加長到長度為100點(diǎn)序列用DFT估計(jì)的頻譜。要求畫出相應(yīng)波形。程序代碼：n=0:10;x=cos(0.48*pi*n)+cos(0.52*pi*n);y=fft(x);subplot(211);stem(0:10,abs(y),'fill')title('Magnitude')subplot(212);stem(0:10,angle(y)/pi,'fill')ylabel('/\pi')title('Phase')程序結(jié)果：n=0:10;x=cos(0.48*pi*n)+cos(0.52*pi*n);x=[x,zeros(1,89)];y=fft(x);subplot(211);stem(0:99,abs(y),'fill')title('Magnitude')subplot(212);stem(0:99,angle(y)/pi,'fill')ylabel('/\pi')title('Phase')（2）時(shí)，用DFT估計(jì)的頻譜，并畫出波形。程序代碼：n=0:100;x=cos(0.48*pi*n)+cos(0.52*pi*n);y=fft(x);subplot(211);stem(0:100,abs(y),'fill')title('Magnitude')subplot(212);stem(0:100,angle(y)/pi,'fill')ylabel('/\pi')title('Phase')（3）根據(jù)實(shí)驗(yàn)結(jié)果，分析怎樣提高頻譜分辨率。答:用樣點(diǎn)加密和補(bǔ)零充位提高頻譜分辨率。

3.已知信號，其中，，。從的表達(dá)式可以看出，它包含三個(gè)頻率的正弦波，但是，從其時(shí)域波形看來，似乎是一個(gè)正弦信號，利用DFT做頻譜分析，確定合適的參數(shù)，是得到的頻譜的頻率分辨率符合需要。程序代碼：n=0:39;x1=0.15*sin(2*pi*n)-0.1*sin(6*pi*n);X1=fft(x1,40);subplot(311)stem(n,abs(X1),'fill');x2=0.15*sin(2*pi*(1/6)*n)+sin(4*pi*(1/6)*n);X2=fft(x2,40);subplot(312)stem(n,abs(X2),'fill');x3=0.15*sin(2*pi*0.1*n)+sin(4*pi*0.1*n)-0.1*sin(6*pi*0.1*n);X3=fft(x3,40);subplot(313)stem(n,abs(X3),'fill');利用DFT近似分析連續(xù)時(shí)間信號的頻譜（幅度譜）。分析采用不同的采樣間隔和截取的長度進(jìn)行計(jì)算的結(jié)果，并最終確定合適的參數(shù)。實(shí)驗(yàn)結(jié)果：采樣區(qū)間可以為0到100，采樣間隔為2程序代碼：n=-6:0.1:6;x=exp(-0.1*n).*heaviside(n);y=fft(x);subplot(211);stem(-60:60,abs(y),'fill')title('Magnitude')subplot(212);stem(-60:60,angle(y)/pi,'fill')ylabel('/\pi')title('Phase')n=-6:0.2:6;x=exp(-0.1*n).*heaviside(n);y=fft(x);subplot(211);stem(-30:30,abs(y),'fill')title('Magnitude')subplot(212);stem(-30:30,angle(y)/pi,'fill')ylabel('/\pi')title('Phase')五、收獲和體會(huì)通過本次試驗(yàn)我加深了對DFT原理的理解，熟悉了如何應(yīng)用DFT分析信號頻譜的方法，并對試驗(yàn)中出現(xiàn)的問題進(jìn)行解決。應(yīng)用DFT分析信號的頻譜，深刻理解利用DFT分析信號頻譜的原理，分析實(shí)現(xiàn)過程中出現(xiàn)的現(xiàn)象及解決方法。實(shí)驗(yàn)的主要內(nèi)容是DFT對離散信號的分析可以用來近似估計(jì)連續(xù)信號的頻譜，因?yàn)橛邢揲L序列的DFT仍是有限長序列，故其特別適合數(shù)字系統(tǒng)。對于無限長的序列我們可以使用窗函數(shù)截取為有限長序列，之后再用DFT進(jìn)行頻譜分析。在利用DFT分析連續(xù)信號頻譜是，會(huì)出現(xiàn)混疊現(xiàn)象，減少時(shí)域抽樣間隔可以克服這一誤差；頻譜泄露，選用合適的窗函數(shù)可以改變這種誤差；柵欄現(xiàn)象，補(bǔ)零增加序列的長度可以改善這種誤差，而補(bǔ)零增長序列長度時(shí)，同時(shí)也減小的頻譜分析時(shí)的普賢間隔，從而所計(jì)算出的頻譜將會(huì)顯示更多的細(xì)節(jié)。最后，感謝老師的耐心講解和解答問題，希望在大三、大四有多門課程可以上機(jī)實(shí)驗(yàn)。實(shí)驗(yàn)2利用FFT計(jì)算線性卷積一、實(shí)驗(yàn)?zāi)康?.掌握利用FFT計(jì)算線性卷積的原理及具體實(shí)現(xiàn)方法。2.加深理解重疊相加法和重疊保留法。3.考察利用FFT計(jì)算線性卷積各種方法的適用范圍。二、實(shí)驗(yàn)設(shè)備與環(huán)境計(jì)算機(jī)、MATLAB軟件環(huán)境。三、實(shí)驗(yàn)基礎(chǔ)理論1.線性卷積與圓周卷積設(shè)為L點(diǎn)序列，為M點(diǎn)序列，和的線性卷積為的長度為。和的N點(diǎn)圓周卷積為圓周卷積與線性卷積相等而不產(chǎn)生交疊的必要條件為圓周卷積定理：根據(jù)DFT的性質(zhì)，和的N點(diǎn)圓周卷積的DFT等于它們DFT的乘積2.快速卷積快速卷積算法用圓周卷積實(shí)現(xiàn)線性卷積，根據(jù)圓周卷積定理利用FFT算法實(shí)現(xiàn)圓周卷積?？梢詫⒖焖倬矸e的步驟歸納如下：（1）為了使線性卷積可以利用圓周卷積來計(jì)算，必須選擇；同時(shí)為了能使用基2-FFT完成卷積運(yùn)算，要求。采用補(bǔ)零的辦法是和的長度均為N。（2）計(jì)算x(n)和h(n)的N點(diǎn)FFT（3）組成卷積（4）利用IFFT計(jì)算IDFT，得到線性卷積3.分段卷積我們考察單位取樣響應(yīng)為的線性系統(tǒng)，輸入為，輸出為，則當(dāng)輸入序列極長時(shí)，如果要等全部集齊時(shí)再開始進(jìn)行卷積，會(huì)使輸出相對輸入有較大的延時(shí)，再者如果序列太長，需要大量的存儲(chǔ)單元。為此，我們把分段，分別求出每段的卷積，合在一起其到最后的總輸出。這種方法稱為分段卷積。分段卷積可細(xì)分為重疊相加法和重疊保留法。重疊保留法：設(shè)的長度為，的長度為M。我們把序列分成多段N點(diǎn)序列，每段與前一段重疊M-1個(gè)樣本。由于第一段沒有前一段保留信號，為了修正，我們在第一個(gè)輸入段前面填充M-1個(gè)零。計(jì)算每一段與的圓周卷積，則其每段卷積結(jié)果的前M-1個(gè)樣本不等于線性卷積值，不是正確的樣本值。所以我們將每段卷積結(jié)果的前M-1個(gè)樣本舍去，只保留后面的N-M+1個(gè)正確輸出樣本，把這些輸出樣本合起來得到總的輸出。利用FFT實(shí)現(xiàn)重疊保留法的步驟如下：（1）在前面填充M-1個(gè)零，擴(kuò)大以后的序列為（2）將分為若干N點(diǎn)子段，設(shè)為每一段的有效數(shù)據(jù)長度，則第i段的數(shù)據(jù)為（3）計(jì)算每一段與的N點(diǎn)圓周卷積，利用FFT計(jì)算圓周卷積：（4）舍去每一段卷積結(jié)果的前M-1個(gè)樣本，連接剩下樣本，得到卷積結(jié)果。重疊相加法：設(shè)長度為M，將信號分解成長為L的子段，建議L選擇與M數(shù)量級相同。以表示每段信號，則每一段卷積的長度為L+M-1，所以在做求和時(shí)，相鄰兩段序列有M-1個(gè)樣本重疊，即前一段的最后M-1個(gè)樣本和下一段的前M-1個(gè)序列重疊，這個(gè)重疊部分相加，再與不重疊部分共同組成輸出。利用FFT實(shí)現(xiàn)重疊保留法的步驟如下：（1）將分為若干L點(diǎn)子段。（2）計(jì)算每一段與的卷積，根據(jù)快速卷積算法利用FFT計(jì)算卷積。（3）將各段（包括重疊部分）相加，得到輸出4.可能用到的MATLAB函數(shù)實(shí)驗(yàn)中FFT運(yùn)算可采用MATLAB中提供的函數(shù)fft來實(shí)現(xiàn)。MATLAB提供了函數(shù)conv來計(jì)算線性卷積，實(shí)驗(yàn)中可以將編程計(jì)算的結(jié)果和conv函數(shù)的計(jì)算結(jié)果相比較，以驗(yàn)證結(jié)果的正確性，conv函數(shù)用法如下：x和h為要進(jìn)行卷積運(yùn)算的兩個(gè)序列，y為卷積結(jié)果統(tǒng)計(jì)程序運(yùn)行時(shí)間可以利用MATLAB提供的tic和toc兩個(gè)命令，具體用法如下：tic……需要運(yùn)行的程序代碼toc即在需要統(tǒng)計(jì)運(yùn)行時(shí)間的程序代碼前加上tic命令，之后加上toc命令，此時(shí)會(huì)在命令窗口中顯示該程序的運(yùn)行時(shí)間。四、實(shí)驗(yàn)內(nèi)容及結(jié)果分析假設(shè)要計(jì)算序列x(n)=u(n)-u(n-L)，0≤n≤L和的線性卷積，完成以下實(shí)驗(yàn)內(nèi)容：1.設(shè)L=M，根據(jù)線性卷積表達(dá)式和快速卷積原理，分別編程實(shí)現(xiàn)計(jì)算兩個(gè)序列線性卷積的方法，比較當(dāng)序列長度分別為8，16，32，64，256，512，1024時(shí)，兩種方法計(jì)算線性卷積所用的時(shí)間。（1）當(dāng)長度為8時(shí):程序代碼：L=8;M=8;n1=0:L-1;x=ones(1，L);n2=0:M-1;h=cos(0.2*pi*n2);y=conv(x，h);n3=0:L+M-2;stem(n3，y)實(shí)驗(yàn)結(jié)果：程序代碼：L=8;M=8;n1=0:L-1;x=ones(1，L);n2=0:M-1;h=cos(0.2*pi*n2);ticy=conv(x，h);toc實(shí)驗(yàn)結(jié)果：Elapsedtimeis0.000062seconds.（2）當(dāng)長度為16時(shí):程序代碼：L=16;M=16;n1=0:L-1;x=ones(1,L);n2=0:M-1;h=cos(0.2*pi*n2);y=conv(x,h);n3=0:L+M-2;stem(n3,y)實(shí)驗(yàn)結(jié)果：程序代碼：clearclcL=16;M=16;n1=0:L-1;x=ones(1,L);n2=0:M-1;h=cos(0.2*pi*n2);ticy=conv(x,h);toc實(shí)驗(yàn)結(jié)果：Elapsedtimeis0.000077seconds.（3）當(dāng)長度為32時(shí)：程序代碼：L=32;M=32;n1=0:L-1;x=ones(1,L);n2=0:M-1;h=cos(0.2*pi*n2);y=conv(x,h);n3=0:L+M-2;stem(n3,y)實(shí)驗(yàn)結(jié)果：程序代碼：clearclcL=32;M=32;n1=0:L-1;x=ones(1,L);n2=0:M-1;h=cos(0.2*pi*n2);ticy=conv(x,h);toc實(shí)驗(yàn)結(jié)果：Elapsedtimeis0.000083seconds.（4）當(dāng)長度為64時(shí)：程序代碼：L=64;M=64;n1=0:L-1;x=ones(1，L);n2=0:M-1;h=cos(0.2*pi*n2);y=conv(x,h);n3=0:L+M-2;stem(n3,y)實(shí)驗(yàn)結(jié)果：程序代碼：clearclcL=64;M=64;n1=0:L-1;x=ones(1,L)；n2=0:M-1;h=cos(0.2*pi*n2);ticy=conv(x,h)；toc實(shí)驗(yàn)結(jié)果：Elapsedtimeis0.000082seconds.（5）當(dāng)長度為256時(shí)：程序代碼：L=256;M=256;n1=0:L-1;x=ones(1,L);n2=0:M-1;h=cos(0.2*pi*n2);y=conv(x,h);n3=0:L+M-2;stem(n3,y)實(shí)驗(yàn)結(jié)果：程序代碼：clearclcL=256;M=256;n1=0:L-1；x=ones(1,L);n2=0:M-1;h=cos(0.2*pi*n2);ticy=conv(x,h);toc實(shí)驗(yàn)結(jié)果：Elapsedtimeis8.800034seconds.（6）當(dāng)長度為512時(shí)：程序代碼：L=512;M=512;n1=0:L-1;x=ones(1,L);n2=0:M-1;h=cos(0.2*pi*n2);y=conv(x,h);n3=0:L+M-2;stem(n3,y)實(shí)驗(yàn)圖形：程序代碼：clearclcL=512;M=512;n1=0:L-1;x=ones(1,L);n2=0:M-1;h=cos(0.2*pi*n2);ticy=conv(x,h);toc實(shí)驗(yàn)結(jié)果：Elapsedtimeis0.000221seconds.（7）長度為1024：程序代碼：L=1024;M=1024;n1=0:L-1;x=ones(1,L);n2=0:M-1;h=cos(0.2*pi*n2);y=conv(x,h);n3=0:L+M-2;stem(n3,y)實(shí)驗(yàn)圖形：程序代碼：clcclearL=1024;M=1024;n1=0:L-1;x=ones(1,L);n2=0:M-1;h=cos(0.2*pi*n2);ticy=conv(x,h);toc實(shí)驗(yàn)結(jié)果：Elapsedtimeis0.000562seconds.快速卷積：程序代碼：i=input('Pleaseinputthenumberofdegree:');L=i;M=i;x=[ones(1,L),zeros(1,M-1)];X=fft(x);n2=0:M-1;h=[cos(0.2*pi*n2),zeros(1,L-1)];H=fft(h);Y=X.*H;y=ifft(Y);n3=0:L+M-2;stem(n3,y)①實(shí)驗(yàn)圖形：（1）當(dāng)長度為度為8：（2）當(dāng)長度為16：（3）當(dāng)長度為32：（4）當(dāng)長度為64：（5）當(dāng)長度為256：（6）當(dāng)長度為512：（7）當(dāng)長度為1024：程序代碼：clccleari=input('Pleaseinputthenumberofdegree:');L=i;M=i;x=[ones(1,L),zeros(1,M-1)];X=fft(x);n2=0:M-1;h=[cos(0.2*pi*n2),zeros(1,L-1)];H=fft(h);ticY=X.*H;toc實(shí)驗(yàn)結(jié)果：（1）長度為8：Pleaseinputthenumberofdegree:8Elapsedtimeis0.000008seconds.（2）長度為16：Pleaseinputthenumberofdegree:16Elapsedtimeis0.000009seconds.（3）長度為32：Pleaseinputthenumberofdegree:32Elapsedtimeis0.000012seconds.（4）長度為64：Pleaseinputthenumberofdegree:64Elapsedtimeis0.000008seconds.（5）長度為256：Pleaseinputthenumberofdegree:256Elapsedtimeis0.000015seconds.（6）長度為512：Pleaseinputthenumberofdegree:512Elapsedtimeis0.000027seconds.（7）長度為1024：Pleaseinputthenumberofdegree:1024Elapsedtimeis0.000037seconds.結(jié)論：快速卷積所用的時(shí)間要比線性卷積所用的時(shí)間少2.當(dāng)L=2048且M=256時(shí)，比較直接計(jì)算線性卷積和快速卷積所需的時(shí)間，進(jìn)一步考察當(dāng)L=4096且M=256時(shí)兩種算法所需的時(shí)間。線性卷積代碼：clcclearL=input('PleaseinputtheL:');M=input('PleaseinputtheM:');n1=0:L-1;x=ones(1,L);n2=0:M-1;h=cos(0.2*pi*n2);ticy=conv(x,h);tocn3=0:L+M-2;stem(n3,y)實(shí)驗(yàn)結(jié)果：①L=2048且M=256PleaseinputtheL:2048PleaseinputtheM:256Elapsedtimeis0.011759seconds.②L=4096且M=256PleaseinputtheL:4096PleaseinputtheM:256Elapsedtimeis0.000718seconds.快速卷積代碼：clcclearL=input('PleaseinputtheL:');M=input('PleaseinputtheM:');x=[ones(1,L),zeros(1,M-1)];X=fft(x);n2=0:M-1;h=[cos(0.2*pi*n2),zeros(1,L-1)];H=fft(h);tic Y=X.*H;tocy=ifft(Y);n3=0:L+M-2;stem(n3,y)實(shí)驗(yàn)結(jié)果：①L=2048且M=256PleaseinputtheL:2048PleaseinputtheM:256Elapsedtimeis0.000103seconds.②L=4096且M=256PleaseinputtheL:4096PleaseinputtheM:256Elapsedtimeis0.000059seconds.3.編程實(shí)現(xiàn)利用重疊相加法計(jì)算兩個(gè)序列的線性卷積，考察L=2048且M=256時(shí)計(jì)算線性卷積的時(shí)間，與第二題結(jié)果進(jìn)行比較。程序代碼：clcclearL=2048;M=256;x=ones(1,L);n=0:M-1;h=cos(0.2*pi*n);N=L+M-1;y=zeros(1,N);subLen=M;subNum=L/M;subLenY=2*M-1;ticH=fft(h,2*M-1);forcount=1:subNumx_temp=x(((count-1)*subLen+1):(count*subLen));X=fft(x_temp,2*M-1);Y=X.*H;y_temp=ifft(Y,2*M-1);y(((count-1)*subLen+1):(count*subLen+M-1))=y(((count-1)*subLen+1):(count*subLen+M-1))+y_temp;endtoc實(shí)驗(yàn)結(jié)果：Elapsedtimeis0.001721seconds.4.編程實(shí)現(xiàn)利用重疊保留法計(jì)算兩個(gè)序列的線性卷積，考察L=2048且M=256時(shí)計(jì)算線性卷積的時(shí)間，與第二題結(jié)果進(jìn)行比較。程序代碼：clcclearL=2048;M=256;x=ones(1,L);n=0:M-1;h=cos(0.2*pi*n);N=L+M-1;useLen=M;subLen=useLen+M-1;subNum=round(L/useLen);x_bar=[zeros(1,M-1),x];ticH=fft(h,subLen);forcount=1:subNumx_temp=x_bar(((count-1)*useLen+1):(count-1)*useLen+subLen);X=fft(x_temp,subLen);Y=X.*H;y_temp=ifft(Y,subLen);y(((count-1)*useLen+1):(count*useLen))=y_temp(M:subLen);endx_temp=x_bar((count*useLen+1):N);X=fft(x_temp,subLen);Y=X.*H;y_temp=ifft(Y,subLen);y((count*useLen+1):N)=y_temp(M:2*M-2);toc實(shí)驗(yàn)結(jié)果：Elapsedtimeis0.000906seconds.實(shí)驗(yàn)結(jié)果分析：1.當(dāng)L和M的數(shù)量級相當(dāng)，且L=M比較大時(shí)，圓周卷積的時(shí)間比線性卷積有明顯的優(yōu)勢，當(dāng)M大于32時(shí)，優(yōu)勢以及相當(dāng)明顯。 2.當(dāng)L和M相差很大時(shí)，直接進(jìn)行圓周卷積的時(shí)間則明顯增加。 3.重疊相加法和重疊保留法將L分為n個(gè)和M同樣數(shù)量級的子段，可以顯著的提高運(yùn)算效率。五、實(shí)驗(yàn)心得與體會(huì)通過本次實(shí)驗(yàn)，我熟練掌握利用FFT計(jì)算線性卷積的原理及具體實(shí)現(xiàn)方法，加深理解重疊相加法和重疊保留法。利用FFT計(jì)算線性卷積各種方法的適用范圍。并且通過自己實(shí)際的編程實(shí)現(xiàn)了重疊相加法和重疊保留法兩種方法，進(jìn)一步加深了對兩種方法的理解及運(yùn)用，并且對課本上的重疊保留和重疊相加法的理解更加的透徹，但實(shí)現(xiàn)編程時(shí)編程時(shí)遇到了很多問題，后來經(jīng)過請教老師和同學(xué)得到了解決，所以實(shí)驗(yàn)前一定要對算法原理有充分的熟悉才能在實(shí)驗(yàn)中比較流暢地寫出算法。最后，感謝老師的耐心講解和解答問題，希望在大三、大四有多門課程可以上機(jī)實(shí)驗(yàn)。實(shí)驗(yàn)三IIR數(shù)字濾波器設(shè)計(jì)一、實(shí)驗(yàn)?zāi)康?.掌握利用脈沖響應(yīng)不變法和雙線性變換法設(shè)計(jì)IIR數(shù)字濾波器的原理及具體方法。2.加深理解數(shù)字濾波器和模擬濾波器之間的技術(shù)指標(biāo)轉(zhuǎn)化。3.掌握脈沖響應(yīng)不變法和雙線性變換法設(shè)計(jì)IIR數(shù)字濾波器的優(yōu)缺點(diǎn)及適用范圍。二、實(shí)驗(yàn)設(shè)備與環(huán)境計(jì)算機(jī)、MATLAB軟件環(huán)境。三、實(shí)驗(yàn)基礎(chǔ)理論1、基本原理從時(shí)域響應(yīng)出發(fā)，使數(shù)字濾波器的單位脈沖響應(yīng)h(n)模仿模擬濾波器的單位沖激響應(yīng)，h(n)等于的取樣值。2、變換方法思路：(1)將進(jìn)行部分分式展開(2)對進(jìn)行拉式變換(3)對時(shí)域采樣得到h(n)(4)對h(n)進(jìn)行z變換3、設(shè)計(jì)步驟(1)確定數(shù)字濾波器的性能指標(biāo)。(2)將數(shù)字濾波器頻率指標(biāo)轉(zhuǎn)換成響應(yīng)的模擬濾波器頻率指標(biāo)(3)根據(jù)指標(biāo)，，和設(shè)計(jì)模擬濾波器。(4)將展成部分分式形式。(5)把模擬極點(diǎn)轉(zhuǎn)換成數(shù)字極點(diǎn)，得到數(shù)字濾波器?？梢娭罤(z)間的變換關(guān)系為：方法1：利用residue函數(shù)和residuez函數(shù)實(shí)現(xiàn)脈沖響應(yīng)不變變換法，實(shí)用方法如下：[r,p,k]=residue(b,a)[b,a]=residue(r,p,k)實(shí)現(xiàn)多項(xiàng)式形式和部分分式形式之間的轉(zhuǎn)換[r,p,k]=residuez(b,a)[b,a]=residuez(r,p,k)實(shí)現(xiàn)多項(xiàng)式形式和部分分式形式之間的轉(zhuǎn)換方法2：MATLAB中提供了impinvar函數(shù)采用脈沖響應(yīng)不變法實(shí)現(xiàn)模擬濾波器到數(shù)字濾波器的變換，其使用如下：[bz,az]=impinvar(b,a,fs)采用脈沖響應(yīng)不變法將模擬濾波器系統(tǒng)函數(shù)的系數(shù)向量b和a變換成為數(shù)字濾波器系統(tǒng)函數(shù)的系數(shù)向量bz和az，fs為采樣頻率（默認(rèn)為1）。[bz,az]=impinvar(b,a)采樣頻率默認(rèn)為1的情況下，采用脈沖響應(yīng)不變法將模擬濾波器變換為數(shù)字濾波器。四、實(shí)驗(yàn)內(nèi)容及結(jié)果分析設(shè)采樣頻率為：，設(shè)計(jì)數(shù)字低通濾波器，滿足如下指標(biāo)：通帶截止頻率：，通帶波動(dòng)：阻帶截止頻率：，阻帶衰減：要求分別設(shè)計(jì)巴特沃斯、切比雪夫Ⅰ型、切比雪夫Ⅱ型和橢圓模擬原型濾波器，并分別結(jié)合脈沖響應(yīng)不變法和雙線性變換法進(jìn)行設(shè)計(jì)。結(jié)合實(shí)驗(yàn)結(jié)果，分別討論采用上述方法設(shè)計(jì)的數(shù)字濾波器是否都能滿足給定指標(biāo)要求，分析脈沖響應(yīng)不變法和雙線性變換法設(shè)計(jì)IIR數(shù)字濾波器的優(yōu)缺點(diǎn)和適用范圍1、巴特沃斯型濾波器：程序代碼：wp=1000*2*pi;ws=1500*2*pi;Rp=1;As=15;N=ceil((log10((10^(Rp/10)-1)/(10^(As/10)-1)))/(2*log10(wp/ws)));%階數(shù)OmegaC=wp/((10^(Rp/10)-1)^(1/(2*N)));[b,a]=butter(N,OmegaC,'s');fs=10000;[bz,az]=impinvar(b,a,fs);[H,w]=freqz(bz,az);subplot(221);plot(w/pi,abs(H));gridon;xlabel('\omega(\pi)');ylabel('|H(e^j^\omega)|');subplot(222);plot(w/pi,20*log10(abs(H)/max(abs(H))));gridon;xlabel('\omega(\pi)');ylabel('|H(e^j^\omega)|(dB)');subplot(223);plot(w/pi,angle(H)/pi);gridon;xlabel('\omega(\pi)');ylabel('PhaseofH(e^j^\omega)(\pi)');grd=grpdelay(bz,az,w);subplot(224);plot(w/pi,grd);gridon;xlabel('\omega(\pi)');ylabel('ImpulseResponse');程序結(jié)果：將數(shù)字化的代碼修改為：[bz,az]=bilinear(b,a,fs);用雙線性變換法改進(jìn)后的結(jié)果：2、切比雪夫Ⅰ型：程序代碼：Rp=1;As=15;wp=1000*2*pi;ws=1500*2*pi;E=((10^(Rp/10))-1)^0.5;A=10^(As/20);N=ceil((acosh(((A^2-1)^0.5)/E))/(acosh(ws/wp)));[b,a]=cheby1(N,Rp,wp,'s');fs=10000;[bz,az]=impinvar(b,a,fs);[H,w]=freqz(bz,az);subplot(221);plot(w/pi,abs(H));gridon;xlabel('\omega(\pi)');ylabel('|H(e^j^\omega)|');subplot(222);plot(w/pi,20*log10(abs(H)/max(abs(H))));gridon;xlabel('\omega(\pi)');ylabel('|H(e^j^\omega)|(dB)');subplot(223);plot(w/pi,angle(H)/pi);gridon;xlabel('\omega(\pi)');ylabel('PhaseofH(e^j^\omega)(\pi)');grd=grpdelay(bz,az,w);subplot(224);plot(w/pi,grd);gridon;xlabel('\omega(\pi)');ylabel('ImpulseResponse');程序結(jié)果：將數(shù)字化的代碼修改為：[bz,az]=bilinear(b,a,fs);用雙線性變換法改進(jìn)后的結(jié)果：切比雪夫II型低通濾波器程序代碼：Rp=1;As=15;wp=1000*2*pi;ws=1500*2*pi;E=((10^(Rp/10))-1)^0.5;A=10^(As/20);N=ceil((acosh(((A^2-1)^0.5)/E))/(acosh(ws/wp)));fs=10000;[b,a]=cheby2(N,As,ws,'s');%求切比雪夫II型濾波器[bz,az]=impinvar(b,a,fs);%將模擬濾波器數(shù)字化[H,w]=freqz(bz,az); %求頻率響應(yīng)subplot(221);plot(w/pi,abs(H));gridon;xlabel('\omega(\pi)');ylabel('|H(e^j^\omega)|');subplot(222);plot(w/pi,20*log10(abs(H)/max(abs(H))));gridon;xlabel('\omega(\pi)');ylabel('|H(e^j^\omega)|(dB)');subplot(223);plot(w/pi,angle(H)/pi);gridon;xlabel('\omega(\pi)');ylabel('PhaseofH(e^j^\omega)(\pi)');grd=grpdelay(bz,az,w);subplot(224);plot(w/pi,grd);gridon;xlabel('\omega(\pi)');ylabel('ImpulseResponse');程序結(jié)果：同理用雙線性變換法改進(jìn)后的結(jié)果：橢圓模擬原型濾波器程序代碼：Rp=1; %通帶波動(dòng)As=15; %阻帶波動(dòng)wp=1000*2*pi; %通帶截止頻率ws=1500*2*pi; %阻帶截止頻率E=((10^(Rp/10))-1)^0.5;%通帶波動(dòng)系數(shù)A=10^(As/20); %阻帶衰減參數(shù)k=wp/ws; %增益k1=E/((A^2-1)^0.5);N=ceil(ellipke(k)*ellipke((1-k1^2)^0.5)/(ellipke(k1)*ellipke((1-k^2)^0.5))); fs=10000;[b,a]=ellip(N,Rp,As,wp,'s');%橢圓濾波器[bz,az]=impinvar(b,a,fs);%將模擬濾波器數(shù)字化[H,w]=freqz(bz,az);%求頻率響應(yīng)subplot(221);plot(w/pi,abs(H));gridon;xlabel('\omega(\pi)');ylabel('|H(e^j^\omega)|');subplot(222);plot(w/pi,20*log10(abs(H)/max(abs(H))));gridon;xlabel('\omega(\pi)');ylabel('|H(e^j^\omega)|(dB)');subplot(223);plot(w/pi,angle(H)/pi);gridon;xlabel('\omega(\pi)');ylabel('PhaseofH(e^j^\omega)(\pi)');grd=grpdelay(bz,az,w);subplot(224);plot(w/pi,grd); gridon;xlabel('\omega(\pi)');ylabel('ImpulseResponse');程序結(jié)果：同理用雙線性變換法改進(jìn)后的結(jié)果：結(jié)合實(shí)驗(yàn)結(jié)果，分別討論采用上述設(shè)計(jì)的數(shù)字濾波器是否都能滿足給定指標(biāo)要求，分析脈沖響應(yīng)不變法和雙線性變換法設(shè)計(jì)IIR數(shù)字濾波器。巴特沃斯低通濾波器的通帶滿足指標(biāo)，在1kHz的時(shí)候衰減1dB，滿足通帶波動(dòng)；阻帶滿足指標(biāo)，在1.5kHz滿足衰減15dB；。切比雪夫I型低通濾波器在通帶滿足指標(biāo)，在0到1kHz的時(shí)候通帶波動(dòng)1dB；阻帶滿足指標(biāo)，在1.5kHz滿足衰減21.5dB。切比雪夫II型低通濾波器通帶和阻帶都不滿足指標(biāo)。通帶波動(dòng)6.6dB，阻帶衰減8.8dB。橢圓低通濾波器通帶不滿足指標(biāo)，通帶波動(dòng)3.5dB，阻帶不滿足指標(biāo)，阻帶衰減12.5dB。用脈沖響應(yīng)不變法設(shè)計(jì)IIR數(shù)字濾波器的優(yōu)缺點(diǎn)及適用范圍：優(yōu)點(diǎn)：①頻率坐標(biāo)的變換是線性的：；②模擬濾波器時(shí)域的優(yōu)點(diǎn)會(huì)在數(shù)字濾波器中保留下來。比如若模擬濾波器具有良好的階躍響應(yīng)特性，如小的上升時(shí)間和低的過沖鋒值等，這些特性會(huì)在數(shù)字濾波器中保留下來。缺點(diǎn)：因?yàn)槟M濾波器都不能是完全帶限的，所以會(huì)引起混疊。②頻率響應(yīng)與T成反比。采樣頻率很高時(shí)，可能出現(xiàn)所不希望得到的高增益。脈沖響應(yīng)不變法設(shè)計(jì)IIR數(shù)字濾波器因?yàn)樗淖畲笕秉c(diǎn)是由頻率響應(yīng)的混疊效應(yīng)。只有當(dāng)模擬濾波器的頻響是限帶于折疊頻率以內(nèi)時(shí)，才能使數(shù)字濾波器的頻率響應(yīng)在折疊頻率以內(nèi)，重現(xiàn)模擬濾波器的頻率響應(yīng)而不產(chǎn)生混疊失真。所以，此法只適用于限帶的模擬濾波器，而對于對于高通和帶阻濾波器，無論取樣周期為多少，無法滿足這一條件。五、實(shí)驗(yàn)心得與體會(huì)通過本次實(shí)驗(yàn)，我熟練掌握利用脈沖響應(yīng)不變法設(shè)計(jì)IIR數(shù)字濾波器的原理及具體方法和脈沖響應(yīng)不變法設(shè)計(jì)IIR數(shù)字濾波器的優(yōu)缺點(diǎn)及適用范圍，加深了數(shù)字濾波器和模擬濾波器之間的技術(shù)指標(biāo)轉(zhuǎn)化理解，并且對書本上IIR數(shù)字濾波器的講述有了更加深刻的了解。設(shè)計(jì)IIR數(shù)字濾波器的基本原理是從時(shí)域響應(yīng)出發(fā)，使數(shù)字濾波器的單位脈沖響應(yīng)h(n)模仿模擬濾波器的單位沖擊響應(yīng)。實(shí)驗(yàn)初期遇到各個(gè)濾波器指標(biāo)的轉(zhuǎn)化問題，對各個(gè)域中指標(biāo)的區(qū)別不是很清晰，后來經(jīng)過查閱課本慢慢糾正，得到了最后的濾波器，對我進(jìn)一步學(xué)習(xí)濾波器設(shè)計(jì)有很大的提升。通過實(shí)驗(yàn)，我學(xué)習(xí)了利用脈沖響應(yīng)不變法設(shè)計(jì)IIR數(shù)字濾波器的基本原理為：從時(shí)域響應(yīng)出發(fā)，使數(shù)字濾波器的單位脈沖響應(yīng)h(n)模仿模擬濾波器的單位沖激響應(yīng)ha(t)，h(n)等于ha脈沖響應(yīng)不變法設(shè)計(jì)IIR數(shù)字濾波器的優(yōu)點(diǎn)是一個(gè)線性相位的模擬濾波器通過脈沖響應(yīng)不變法得到的仍然是一個(gè)一個(gè)線性相位的數(shù)字濾波器。它的最大缺點(diǎn)是由頻率響應(yīng)的混疊效應(yīng)。所以，此法只適用于限帶的模擬濾波器。通過學(xué)習(xí)實(shí)驗(yàn)教程上的實(shí)例，自己進(jìn)行一個(gè)知識的運(yùn)用就很快就完成了實(shí)驗(yàn)，只是對于切比雪夫II型低通濾波器和橢圓低通濾波器的原型模擬濾波器設(shè)計(jì)出來的頻率響應(yīng)曲線不能理解，自己一度認(rèn)為是錯(cuò)的，后來向老師請教得到，這個(gè)就是這樣的，只是這兩個(gè)不能使用脈沖響應(yīng)不變法來設(shè)計(jì)IIR數(shù)字濾波器。這就是在分析實(shí)驗(yàn)結(jié)果的時(shí)候就開始懂得了，為什么會(huì)出現(xiàn)這種現(xiàn)象。這是因?yàn)槊}沖響應(yīng)不變法它的主要缺點(diǎn)就是實(shí)際上沒有一個(gè)是帶限的理想模擬濾波器，這樣就一定會(huì)產(chǎn)生折疊效應(yīng)。最后，感謝老師的耐心講解和解答問題，希望在大三、大四有多門課程可以上機(jī)實(shí)驗(yàn)。實(shí)驗(yàn)四窗函數(shù)法設(shè)計(jì)FIR數(shù)字濾波器一、實(shí)驗(yàn)?zāi)康恼莆沾昂瘮?shù)法設(shè)計(jì)FIR數(shù)字濾波器的原理及具體方法。二、實(shí)驗(yàn)設(shè)備與環(huán)境計(jì)算機(jī)、MATLAB軟件環(huán)境。三、實(shí)驗(yàn)原理1、基本原理 窗函數(shù)設(shè)計(jì)法的基本思想為，首先選擇一個(gè)適當(dāng)?shù)睦硐氲臑V波器，然后用窗函數(shù)截取它的單位脈沖響應(yīng)，得到線性相位和因果的FIR濾波器。這種方法的重點(diǎn)是選擇一個(gè)合適的窗函數(shù)和理想濾波器，使設(shè)計(jì)的濾波器的單位脈沖響應(yīng)逼近理想濾波器的單位脈沖響應(yīng)。2、設(shè)計(jì)步驟 （1）給定理想濾波器的頻率響應(yīng)，在通帶上具有單位增益和線性相位，在阻帶上具有零響應(yīng)。一個(gè)帶寬為的低通濾波器由下式給定：其中為采樣延遲，其作用是為了得到一個(gè)因果系統(tǒng)。（2）確定這個(gè)濾波器的單位脈沖響應(yīng)為了得到一個(gè)長度為N的因果的線性相位FIR濾波器，我們令（3）用窗函數(shù)截取得到所設(shè)計(jì)FIR數(shù)字濾波器：3、窗函數(shù)的選擇常用的窗函數(shù)有矩形（Rectangular）窗，漢寧（Hanning）窗，海明（Hamming）窗、布萊克曼（Blackman）窗、凱瑟（Kaiser）窗等。MATLAB提供了一些函數(shù)用于產(chǎn)生窗函數(shù)。如表7-1所示。表7-1MATLAB中產(chǎn)生窗函數(shù)的命令MATLAB函數(shù)窗函數(shù)MATLAB函數(shù)窗函數(shù)Boxcar矩形窗函數(shù)Blackman布萊克曼窗Hanning漢寧窗函數(shù)Kaiser凱瑟窗函數(shù)Hamming海明窗在設(shè)計(jì)過程中我們需要根據(jù)給定的濾波器技術(shù)指標(biāo)，選擇濾波器長N和窗函數(shù)。表7-2列出了常用窗函數(shù)的一些特性，可供設(shè)計(jì)參考。表7-2常用窗函數(shù)的特性窗函數(shù)窗函數(shù)頻率特性加窗后濾波器指標(biāo)旁瓣峰值dB主瓣寬度過渡帶寬最小阻帶衰減dB矩形窗-134π/N1.8π/N-21漢寧窗-318π/N6.2π/N-44海明窗-418π/N6.6π/N-53布萊克曼窗-5712π/N11π/N-74凱瑟窗是一種廣泛在實(shí)際中廣泛應(yīng)用的窗函數(shù)，它由下式給定：其中是修正的零階貝塞爾函數(shù)，參數(shù)控制最小阻帶衰減，這種窗函數(shù)對于相同的N可以提供不同的過渡帶寬。由于貝塞爾函數(shù)比較復(fù)雜，這種窗函數(shù)的設(shè)計(jì)方程很難推導(dǎo)，然而幸運(yùn)的是，有一些經(jīng)驗(yàn)設(shè)計(jì)方程可以直接使用。已知給定的指標(biāo)，濾波器長度N和凱瑟窗參數(shù)可以按如下凱瑟窗方程給出過渡帶帶寬：四、實(shí)驗(yàn)內(nèi)容及結(jié)果分析1.設(shè)計(jì)一個(gè)數(shù)字低通FIR濾波器，其技術(shù)指標(biāo)如下：，，分別采用矩形窗，漢寧窗，海明窗，布萊克曼窗，凱瑟窗設(shè)計(jì)該濾波器，結(jié)合實(shí)驗(yàn)結(jié)果，分別討論采用上述方法設(shè)計(jì)的數(shù)字濾波器是否滿足給定的指標(biāo)要求。矩形窗：程序代碼：wp=0.2*pi;wst=0.3*pi;tr_width=wst-wp;N=ceil(1.8*pi/tr_width)+1;N n=0:N-1;wc=wp+wst;wc=wc/2;alpha=(N-1)/2;hd=(wc/pi)*sinc((wc/pi)*(n-alpha));w_boxcar=boxcar(N)';h=hd.*w_boxcar;subplot(221);stem(n,hd,'filled');axistight;xlabel('n');ylabel('hd(n)');[Hrw1]=zerophase(h);subplot(222);plot(w1/pi,Hr);axistight;xlabel('\omega/\pi');ylabel('H(\omega)');subplot(223);stem(n,h,'filled');axistight;xlabel('n');ylabel('h(n)');[H,w]=freqz(h,1);subplot(224);plot(w/pi,20*log10(abs(H)/max(abs(H))));xlabel('\omega/\pi');ylabel('dB');gridon實(shí)驗(yàn)結(jié)果：N=19結(jié)果分析：阻帶衰減為20dB，不滿足設(shè)計(jì)指標(biāo)。漢寧窗：程序代碼：wp=0.2*pi;wst=0.3*pi;tr_width=wst-wp;N=ceil(6.2*pi/tr_width)+1;Nn=0:N-1;wc=wp+wst;wc=wc/2;alpha=(N-1)/2;hd=(wc/pi)*sinc((wc/pi)*(n-alpha));w_hanning=hanning(N)';h=hd.*w_hanning;subplot(221);stem(n,hd,'filled');axistight;xlabel('n');ylabel('hd(n)');[Hrw1]=zerophase(h);subplot(222);plot(w1/pi,Hr);axistight;xlabel('\omega/\pi');ylabel('H(\omega)');subplot(223);stem(n,h,'filled');axistight;xlabel('n');ylabel('h(n)');[H,w]=freqz(h,1);subplot(224);plot(w/pi,20*log10(abs(H)/max(abs(H))));xlabel('\omega/\pi');ylabel('dB');gridon實(shí)驗(yàn)結(jié)果N=63結(jié)果分析：阻帶衰減為50dB，符合設(shè)計(jì)指標(biāo)。海明窗：程序代碼：wp=0.2*pi;wst=0.3*pi;tr_width=wst-wp;N=ceil(6.6*pi/tr_width)+1;Nn=0:N-1;wc=wp+wst;wc=wc/2;alpha=(N-1)/2;hd=(wc/pi)*sinc((wc/pi)*(n-alpha));w_hamming=hamming(N)';h=hd.*w_hamming;subplot(221);stem(n,hd,'filled');axistight;xlabel('n');ylabel('hd(n)');[Hrw1]=zerophase(h);subplot(222);plot(w1/pi,Hr);axistight;xlabel('\omega/\pi');ylabel('H(\omega)');subplot(223);stem(n,h,'filled');axistight;xlabel('n');ylabel('h(n)');[H,w]=freqz(h,1);subplot(224);plot(w/pi,20*log10(abs(H)/max(abs(H))));xlabel('\omega/\pi');ylabel('dB');gridon實(shí)驗(yàn)結(jié)果N=67結(jié)果分析：阻帶衰減約為50dB，濾波器滿足給定的設(shè)計(jì)指標(biāo)。布萊克曼窗：程序代碼：wp=0.2*pi;wst=0.3*pi;tr_width=wst-wp;N=ceil(11*pi/tr_width)+1;Nn=0:N-1;wc=wp+wst;wc=wc/2;alpha=(N-1)/2;hd=(wc/pi)*sinc((wc/pi)*(n-alpha));w_blackman=blackman(N)';h=hd.*w_blackman;subplot(221);stem(n,hd,'filled');axistight;xlabel('n');ylabel('hd(n)');[Hrw1]=zerophase(h);subplot(222);plot(w1/pi,Hr);axistight;xlabel('\omega/\pi');ylabel('H(\omega)');subplot(223);stem(n,h,'filled');axistight;xlabel('n');ylabel('h(n)');[H,w]=freqz(h,1);subplot(224);plot(w/pi,20*log10(abs(H)/max(abs(H))));xlabel('\omega/\pi');ylabel('dB');gridon實(shí)驗(yàn)結(jié)果：N=111結(jié)果分析：阻帶衰減約為80dB，濾波器滿足給定的設(shè)計(jì)指標(biāo)。凱瑟窗：程序代碼：wp=0.2*pi;wst=0.3*pi;As=50;N=ceil((As-7.95)/(2.285*(wst-wp)))+1;Nn=0:N-1;wc=wp+wst;wc=wc/2;alpha=(N-1)/2;hd=(wc/pi)*sinc((wc/pi)*(n-alpha));w_kaiser=kaiser(N)';h=hd.*w_kaiser;subplot(221);stem(n,hd,'filled');axistight;xlabel('n');ylabel('hd(n)');[Hrw1]=zerophase(h);subplot(222);plot(w1/pi,Hr);axistight;xlabel('\omega/\pi');ylabel('H(\omega)');subplot(223);stem(n,h,'filled');axistight;xlabel('n');ylabel('h(n)');[H,w]=freqz(h,1);subplot(224);plot(w/pi,20*log10(abs(H)/max(abs(H))));xlabel('\omega/\pi');ylabel('dB');gridon實(shí)驗(yàn)結(jié)果：N=60結(jié)果分析：阻帶衰減約為20dB，濾波器不滿足給定的設(shè)計(jì)指標(biāo)。2.設(shè)計(jì)一個(gè)數(shù)字帶通FIR濾波器，其技術(shù)指標(biāo)如下：下阻帶邊緣：，下通帶邊緣：，上通帶邊緣：，下阻帶邊緣：，程序代碼：wls=0.2*pi;wlp=0.35*pi;whp=0.65*pi;wc=[wlp/pi,whp/pi];B=wlp-wls;N=ceil(8/0.15);n=0:N-1;window=hanning(N);[h1,w]=freqz(window,1);hn=fir1(N-1,wc,hanning(N));[h2,w]=freqz(hn,1,512);subplot(221);stem(window,'filled');axis([06001.2]);gridon;xlabel('n');subplot(222);plot(w/pi,20*log(abs(h1)/abs(h1(1))));axis([01-3500]);gridon;xlabel('w/pi');ylabel('dB');subplot(223);stem(n,hn,'filled');axis([060-0.250.25]);gridon;xlabel('n');ylabel('h(n)');subplot(224);plot(w/pi,20*log(abs(h2)/abs(h2(1))));gridon;xlabel('w/pi');ylabel('dB');實(shí)驗(yàn)結(jié)果：結(jié)果分析：如圖所示，該帶通濾波器的各項(xiàng)指標(biāo)均滿足要求。3.采用頻率抽樣設(shè)計(jì)法設(shè)計(jì)FIR數(shù)字低通濾波器，滿足以下指標(biāo)，，（1）取N=10，過渡帶沒有樣本（2）取N=40，過渡帶有一個(gè)樣本，T=0.39（3）取N=60，過渡帶有兩個(gè)樣本，T1=0.5925，T2=0.1099（4）分別討論采用上述方法設(shè)計(jì)的數(shù)字濾波器是否都能滿足給定指標(biāo)要求。答：（1）取N=20，過渡帶沒有樣本程序代碼：N=20;alpha=(N-1)/2;L=0:N-1;wL=(2*pi/N)*L;Hrs=[1,1,1,zeros(1,15),1,1];Hdr=[1,1,0,0];wdL=[0,0.25,0.25,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=ifft(H,N);w=[0:500]*pi/500;[Hr,wr]=zerophase(h);subplot(221);plot(wdL,Hdr,wL(1:N/2+1)/pi,Hrs(N/2+1),'o');axis([0,1,-0.1,1.1]);xlabel('\omega(\pi)');ylabel('Hr(k)');subplot(222);stem(L,h,'filled');axis([0,N-1,-0.1,0.3]);xlabel('n');ylabel('h(n)');subplot(223);plot(wr/pi,Hr,wL(1:N/2+1)/pi,Hrs(1:N/2+1),'o');axis([0,1,-0.2,1.2]);gridon;xlabel('\omega(\pi)');ylabel('Hr(\omega)');subplot(224);plot(wr/pi,20*log10((abs(Hr)/max(abs(Hr)))));axis([0,1,-50,5]);gridon;xlabel('\omega(\pi)');ylabel('dB');程序結(jié)果：當(dāng)N=20時(shí)通帶波動(dòng)和阻帶衰減都不符合指標(biāo)。（2）當(dāng)N=4時(shí)，過渡帶有一個(gè)樣本，T=0.39程序代碼：N=40;alpha=(N-1)/2;L=0:N-1;wL=(2*pi/N)*L;Hrs=[1,1,1,1,1,0.39,zeros(1,29),0.39,1,1,1,1];Hdr=[1,1,0.39,0,0];wdL=[0,0.2,0.25,0.3,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=ifft(H,N);w=[0:500]*pi/500;[Hr,wr]=zerophase(h);subplot(221);plot(wdL,Hdr,wL(1:N/2+1)/pi,Hrs(N/2+1),'o');axis([0,1,-0.1,1.1]);xlabel('\omega(\pi)');ylabel('Hr(k)');subplot(222);stem(L,h,'filled');axis([0,N-1,-0.1,0.3]);xlabel('n');ylabel('h(n)');subplot(223);plot(wr/pi,Hr,wL(1:N/2+1)/pi,Hrs(1:N/2+1),'o');axis([0,1,-0.2,1.2]);gridon;xlabel('\omega(\pi)');ylabel('Hr(\omega)');subplot(224);plot(wr/pi,20*log10((abs(Hr)/max(abs(Hr)))));axis([0,1,-50,5]);gridon;xlabel('\omega(\pi)');ylabel('dB');程序結(jié)果：N=40時(shí)阻帶衰減指標(biāo)滿足要求，通帶波動(dòng)指標(biāo)不滿足。（3）當(dāng)N=60時(shí)，過渡帶有兩個(gè)樣本，T1=0.5925，T2=0.1099程序代碼：N=60;alpha=(N-1)/2;L=0:N-1;wL=(2*pi/N)*L;Hrs=[1,1,1,1,1,1,1,0.5925,0.1099,zeros(1,43),0.1099,0.5925,1,1,