北京理工大學(xué)數(shù)字信號(hào)處理實(shí)驗(yàn)報(bào)告

1、 本科實(shí)驗(yàn)報(bào)告實(shí)驗(yàn)名稱： 數(shù)字信號(hào)處理實(shí)驗(yàn) 課程名稱：數(shù)字信號(hào)處理實(shí)驗(yàn)實(shí)驗(yàn)時(shí)間：任課教師：實(shí)驗(yàn)地點(diǎn)：4423實(shí)驗(yàn)教師：實(shí)驗(yàn)類型： 原理驗(yàn)證 綜合設(shè)計(jì) 自主創(chuàng)新學(xué)生姓名：學(xué)號(hào)/班級(jí)：組 號(hào)：學(xué) 院：信息與電子學(xué)院同組搭檔：專 業(yè)：信息工程成 績(jī)：實(shí)驗(yàn)1 利用DFT分析信號(hào)頻譜一、實(shí)驗(yàn)?zāi)康?.加深對(duì)DFT原理的理解。2.應(yīng)用DFT分析信號(hào)頻譜。3.深刻理解利用DFT分析信號(hào)頻譜的原理，分析現(xiàn)實(shí)過(guò)程現(xiàn)象及解決辦法。二、實(shí)驗(yàn)原理1、DFT和DTFT的關(guān)系 有限長(zhǎng)序列的離散時(shí)間傅里葉變換在頻率區(qū)間的N個(gè)等分點(diǎn)上的N個(gè)取樣值可以由下式表示： 由上式可知，序列的N點(diǎn)DFT，實(shí)際上就是序列的DTFT在N個(gè)等間

2、隔頻率點(diǎn)上樣本。2、利用DFT求DTFT 方法1：由恢復(fù)出的方法如圖2.1所示： 圖 2.1.由 N點(diǎn)DFT恢復(fù)頻譜DTFT的流程由圖2.1所示流程圖可知：由式2-2可以得到其中為內(nèi)插函數(shù) 方法2：然而在實(shí)際MATLAB計(jì)算中，上訴插值公式不見(jiàn)得是最好的方法。由于DFT是DTFT的取樣值，其相鄰的兩個(gè)頻率樣本點(diǎn)的間距為,所以如果我們?cè)黾訑?shù)據(jù)的長(zhǎng)度N，使得得到的DFT譜線就更加精細(xì)，其包絡(luò)就越接近DTFT的結(jié)果，這樣可以利用DFT來(lái)近似計(jì)算DTFT。如果沒(méi)有更多的數(shù)據(jù)，可以通過(guò)補(bǔ)零來(lái)增加數(shù)據(jù)長(zhǎng)度。3、利用DFT分析連續(xù)時(shí)間信號(hào)的頻譜 采用計(jì)算機(jī)分析連續(xù)時(shí)間信號(hào)的頻譜，第一步就是把連續(xù)時(shí)間信號(hào)離散

3、化，這里需要進(jìn)行連個(gè)操作：一是采樣，二是截?cái)唷?duì)于連續(xù)非周期信號(hào)，按采樣間隔T進(jìn)行采樣，截取長(zhǎng)度為M，那么對(duì)進(jìn)行N點(diǎn)的頻率采樣，得到因此，可以將利用DFT分析連續(xù)非周期信號(hào)頻譜的步驟歸納如下： （1）確定時(shí)域采樣間隔T，得到離散序列； （2）確定截取長(zhǎng)度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ì)算的頻譜，需要注意如下三點(diǎn)問(wèn)題：（1）頻譜混疊。如果不滿足采樣定理的條件，頻譜會(huì)很出現(xiàn)混疊誤差。對(duì)于頻譜無(wú)限寬的信號(hào)，應(yīng)考慮覆蓋大部分主要頻率的范圍。（2）柵欄效

4、應(yīng)和頻譜分辨率。使用DFT計(jì)算頻譜，得到的結(jié)果只是N個(gè)頻譜樣本值，樣本值之間的頻譜是未知的，就像通過(guò)一個(gè)柵欄觀察頻譜，稱為“柵欄效應(yīng)”。頻譜分辨率與記錄長(zhǎng)度成正比，提高頻譜分辨率，就要增加記錄時(shí)間。（3）頻譜泄露。對(duì)于信號(hào)截?cái)鄷?huì)把窗函數(shù)的頻譜會(huì)引入到信號(hào)頻譜中，造成頻譜泄露。解決這問(wèn)題的主要辦法是采用旁瓣小的窗函數(shù)，頻譜泄露和窗函數(shù)均會(huì)引起誤差。 因此，要合理選取采樣間隔和截取長(zhǎng)度，必要時(shí)還需考慮適當(dāng)?shù)拇啊?對(duì)于連續(xù)周期信號(hào)，我們?cè)诓捎糜?jì)算機(jī)進(jìn)行計(jì)算時(shí)，也總是要進(jìn)行截?cái)?，序列總是有限長(zhǎng)的，仍然可以采用上訴方法近似計(jì)算。4、可能用到MATLAB函數(shù)與代碼實(shí)驗(yàn)中的DFT運(yùn)算可以采用MATLAB中提

5、供的FFT來(lái)實(shí)現(xiàn)。DTFT可以利用MATLAB矩陣運(yùn)算的方法進(jìn)行計(jì)算。三、實(shí)驗(yàn)內(nèi)容1. ，完成如下要求：（1）計(jì)算其DTFT，并畫(huà)出區(qū)間的波形。（2）計(jì)算4 點(diǎn)DFT，并把結(jié)果顯示在（1）所畫(huà)的圖形中。（3）對(duì)補(bǔ)零，計(jì)算64 點(diǎn)DFT，并顯示結(jié)果。（4）是否可以由DFT 計(jì)算DTFT，如果可以，請(qǐng)編程實(shí)現(xiàn)。2. 考察序列（1）時(shí)，用DFT 估計(jì)的頻譜；將補(bǔ)零加長(zhǎng)到長(zhǎng)度為100點(diǎn)序列用DFT估計(jì)的頻譜。要求畫(huà)出相應(yīng)波形。（2）時(shí)，用DFT 估計(jì)x(n)的頻譜，并畫(huà)出波形。3. 已知信號(hào) ，其中，。從的表達(dá)式可以看出，它包含三個(gè)頻率的正弦波，但是，從其時(shí)域波形來(lái)看，似乎是一個(gè)正弦信號(hào)，利用DFT做

6、頻譜分析，確定適合的參數(shù)，使得到的頻譜的頻率分辨率符合需要。4.利用DFT近似分析連續(xù)時(shí)間信號(hào)xt=e-0.1tu(t)的頻譜（幅度譜）。分析采用不同的采樣間隔和截取長(zhǎng)度進(jìn)行計(jì)算的結(jié)果，并最終確定合適的參數(shù)。四、實(shí)驗(yàn)代碼及實(shí)驗(yàn)結(jié)果實(shí)驗(yàn)1.1.1實(shí)驗(yàn)結(jié)果：實(shí)驗(yàn)代碼：>> n=0:3;>> x=2 -1 1 1;>> w=-pi:0.01*pi:pi;>> X=x*exp(-j*n'*w);>> subplot(211);>> plot(w,abs(X);>> title('幅度');xla

7、bel('w');ylabel('|X|');>> axis tight;>> subplot(212);>> plot(w,angle(X);>> title('相位');xlabel('w');ylabel('Angle（X）');>> axis tight;實(shí)驗(yàn)1.1.2實(shí)驗(yàn)結(jié)果：實(shí)驗(yàn)代碼： n=0:3;x=2 -1 1 1;w=-pi:0.01*pi:pi;X=x*exp(-j*n'*w);subplot(211);plot(w,abs(X

8、);title('幅度');xlabel('w');ylabel('|X|');axis tight;hold on;subplot(212);plot(w,angle(X);title('相位');xlabel('w');ylabel('Angle（X）');axis tight;hold on;H=fft(x);subplot(211);>> stem(n,abs(H),'filled');>> subplot(212);>> stem(n,a

9、ngle(H),'filled');實(shí)驗(yàn)1.1.3實(shí)驗(yàn)結(jié)果：實(shí)驗(yàn)代碼：>> x=2 -1 1 1 zeros(1,60);>> X=fft(x);>> subplot(211);>> n=0:63;>> stem(n,abs(X),'filled');>> title('幅度');xlabel('n');ylabel('|X|');>> subplot(212);>> stem(n,angle(X),'fille

10、d');title('相位');xlabel('n');ylabel('angle(X)');分析：可以由DFT計(jì)算DTFT。通過(guò)補(bǔ)零加長(zhǎng)序列，提高采樣密度，可以由DFT近似計(jì)算DTFT。實(shí)驗(yàn)1.2.1實(shí)驗(yàn)結(jié)果：實(shí)驗(yàn)代碼：n=0:10;x=cos(0.48*n*pi)+cos(0.52*n*pi);X=fft(x);subplot(211);stem(n,abs(X),'filled');title('幅度');xlabel('n');ylabel('|X|');subplo

11、t(212);stem(n,angle(X),'filled');title('相位');xlabel('n');ylabel('angle(X)');補(bǔ)零加長(zhǎng)：實(shí)驗(yàn)結(jié)果：實(shí)驗(yàn)代碼：h=x zeros(1,89);H=fft(h);n=0:99;subplot(211);stem(n,abs(H),'filled');title('幅度');xlabel('n');ylabel('|H|');subplot(212);stem(n,angle(H),'fill

12、ed');title('相位');xlabel('n');ylabel('angle(H)');實(shí)驗(yàn)1.2.2實(shí)驗(yàn)結(jié)果：實(shí)驗(yàn)代碼：n=0:100;>> x=cos(0.48*n*pi)+cos(0.52*n*pi);>> X=fft(x);>> subplot(211);stem(n,abs(X),'filled');title('幅度');xlabel('n');ylabel('|X|');subplot(212);stem(n,angl

13、e(X),'filled');title('相位');xlabel('n');ylabel('angle(X)');分析：可以通過(guò)增大截取長(zhǎng)度和增加補(bǔ)零的個(gè)數(shù)來(lái)提高頻譜分辨率，但是補(bǔ)零不能夠增加分辨力。實(shí)驗(yàn)1.3實(shí)驗(yàn)結(jié)果：實(shí)驗(yàn)代碼：f1=1;f2=2;f3=3;Fs=500;Tp=1;t=0:1/Fs:Tp;x=0.15*sin(2*pi*f1*t)+sin(2*pi*f2*t)-0.1*sin(2*pi*f3*t);N=501;F=Fs/N;f=0:Fs/(N-1):Fs;X=fft(x,N)；stem(f,abs(X)/250

14、,'filled');axis(0,4,0,1);xlabel('f/Hz');title('Magnitude');分析：通過(guò)選取合適的采樣周期，可以完整的恢復(fù)出原信號(hào)的頻譜波形。實(shí)驗(yàn)1.4Tp=1s，0n100實(shí)驗(yàn)結(jié)果：實(shí)驗(yàn)代碼：n=0:1:100;x=exp(-0.1*n);X=fft(x);stem(n,abs(X),'filled');Tp=5s,0n100實(shí)驗(yàn)結(jié)果：實(shí)驗(yàn)代碼：n=0:5:100;x=exp(-0.1*n);X=fft(x);stem(n,abs(X),'filled');Tp=25s，0

15、n100實(shí)驗(yàn)結(jié)果：實(shí)驗(yàn)代碼：n=0:25:100;x=exp(-0.1*n);X=fft(x);stem(n,abs(X),'filled');Tp=1s，0n50實(shí)驗(yàn)結(jié)果：實(shí)驗(yàn)代碼：n=0:1:50;x=exp(-0.1*n);X=fft(x);stem(n,abs(X),'filled');Tp=1,0n25實(shí)驗(yàn)結(jié)果：實(shí)驗(yàn)代碼：n=0:1:25;x=exp(-0.1*n);X=fft(x);stem(n,abs(X),'filled');分析： 最終確定參數(shù)：tp=5s，0n100。五、心得與體會(huì)通過(guò)本次實(shí)驗(yàn)，我們掌握并加深對(duì)DFT原理的理解

16、并且學(xué)會(huì)應(yīng)用DFT分析信號(hào)頻譜，在此基礎(chǔ)上利用DFT分析信號(hào)頻譜的原理，掌握了利用matlab分析現(xiàn)實(shí)問(wèn)題的步驟及辦法。并且，通過(guò)這次的實(shí)驗(yàn)，對(duì)信號(hào)序列有了更加深刻的認(rèn)識(shí)，單純的一個(gè)信號(hào)序列是沒(méi)有意義的，只有配合他本身的時(shí)間序列才是一個(gè)完整的信號(hào)序列，才可以對(duì)其進(jìn)行分析。實(shí)驗(yàn)2 利用FFT計(jì)算線性卷積（選作）一、實(shí)驗(yàn)?zāi)康?.掌握利用FFT計(jì)算線性卷積的原理及具體實(shí)現(xiàn)方法。2.加深理解重疊相加法和重疊保留法。3.考察利用FFT計(jì)算線性卷積各種方法的適用范圍。二、實(shí)驗(yàn)原理1.線性卷積與圓周卷積設(shè)x(n)為L(zhǎng)點(diǎn)序列，h(n)為M點(diǎn)序列，x(n)和h(n)的線性卷積為 （3-1）的長(zhǎng)度為L(zhǎng)+M-1x(

17、n)和h(n)的圓周卷積為 （3-2）圓周卷積與線性卷積相等而不產(chǎn)生交疊的必要條件為NL+M+1 (3-3)圓周卷積定理：根據(jù)DFT性質(zhì)，x(n)和h(n)的N點(diǎn)圓周卷積的DFT等于它們的DFT的乘積： (3-4)2.快速卷積快速卷積發(fā)運(yùn)用圓周卷積實(shí)現(xiàn)線性卷積，根據(jù)圓周卷積定理利用FFT算法實(shí)現(xiàn)圓周卷積。可將快速卷積運(yùn)算的步驟歸納如下：(1)必須選擇；為了能使用基-2算法，要求。采用補(bǔ)零的辦法使得x(n)和h(n)的長(zhǎng)度均為N。(2)計(jì)算x(n)和h(n)的N點(diǎn)FFT。(3)組成乘積(4)利用IFFT計(jì)算Y(k)的IDFT，得到線性卷積y(n)3.分段卷積我們考察單位取樣響應(yīng)為h(n)的線性系

18、統(tǒng)，輸入為x(n)，輸出為y(n)，則yn=xn*h(n)當(dāng)輸入序列x(n)極長(zhǎng)時(shí)，如果要等x(n)全部集齊時(shí)再開(kāi)始進(jìn)行卷積，會(huì)使輸出有較大延時(shí)；如果序列太長(zhǎng)，需要大量存儲(chǔ)單元。為此，我們把x(n)分段，為別求出每段的卷積，合在一起得到最后的總輸出。這稱為分段卷積。分段卷積可以細(xì)分為重疊保留法和重疊相加法。重疊保留法：設(shè)x(n)的長(zhǎng)度為，h(n)的長(zhǎng)度為M。把序列x(n)分成多段N點(diǎn)序列，每段雨前一段重寫(xiě)M-1個(gè)樣本。并在第一個(gè)輸入段前面補(bǔ)M-1個(gè)零。計(jì)算每一段與h(n)的圓周卷積，其結(jié)果中前M-1個(gè)不等與線性卷積，應(yīng)當(dāng)舍去，只保留后面N-M+1個(gè)正確的輸出樣本，把它們合起來(lái)得到總的輸出。利用

19、FFT實(shí)現(xiàn)重疊保留法的步驟如下：(1)在x(n)前面填充M-1個(gè)零，擴(kuò)大以后的序列為(2)將x(n)分為若干段N點(diǎn)子段，設(shè)L=N-M+1為每一段的有效長(zhǎng)度，則第i段的數(shù)據(jù)為：(3)計(jì)算每一段與h(n)的N點(diǎn)圓周卷積，利用FFT計(jì)算圓周卷積(4)舍去每一段卷積結(jié)果的前M-1個(gè)樣本，連接剩下的樣本得到卷積結(jié)果y(n)。重疊相加法：設(shè)h(n)長(zhǎng)度為M，將信號(hào)x(n)分解成長(zhǎng)為L(zhǎng)的子段。以表示沒(méi)斷信號(hào)，則：每一段卷積的長(zhǎng)度為L(zhǎng)+M-1，所以在做求和時(shí)，相鄰兩段序列由M-1個(gè)樣本重疊，即前一段的最后M-1個(gè)樣本和下一段前M-1個(gè)樣本序列重疊，這個(gè)重疊部分相加，再與不重疊的部分共同組成y(n)。利用FFT

20、實(shí)現(xiàn)重疊保留法的步驟如下：(1)將x(n)分為若干L點(diǎn)子段。(2)計(jì)算每一段與h(n)的卷積，根據(jù)快速卷積法利用FFT計(jì)算卷積。(3)將各段相加，得到輸出y(n)。4、可能得到的MATLAB函數(shù) 實(shí)驗(yàn)中FFT運(yùn)算可采用MATLAB中提供的函數(shù)fft來(lái)實(shí)現(xiàn)。三、實(shí)驗(yàn)內(nèi)容假設(shè)要計(jì)算序列x(n)=u(n)-u(n-L),0nL和h(n)=cos(0.2n),0nM的線性卷積完成以下實(shí)驗(yàn)內(nèi)容。1.設(shè)L=M,根據(jù)線性卷積的表達(dá)式和快速卷積的原理分別編程實(shí)現(xiàn)計(jì)算兩個(gè)序列線性卷積的方法,比較當(dāng)序列長(zhǎng)度分別為8,16,32,64,256,512,1024時(shí)兩種方法計(jì)算線性卷積所需時(shí)間。2當(dāng)L=2048且M=2

21、56時(shí)比較直接計(jì)算線性卷積和快速卷積所需的時(shí)間,進(jìn)一步考察當(dāng)L=4096且M=256時(shí)兩種算法所需的時(shí)間。3. 編程實(shí)現(xiàn)利用重疊相加法計(jì)算兩個(gè)序列的線性卷積,考察L=2048且M=256時(shí)計(jì)算線性卷積的時(shí)間,與2題的結(jié)果進(jìn)行比較。4. 編程實(shí)現(xiàn)利用重疊保留法計(jì)算兩個(gè)序列的線性卷積,考察L=2048且M=256時(shí)計(jì)算線性卷積的時(shí)間,與2題的結(jié)果進(jìn)行比較。四、實(shí)驗(yàn)代碼及實(shí)驗(yàn)結(jié)果實(shí)驗(yàn)2.1實(shí)驗(yàn)代碼：for M=8 16 32 64 256 512 1024;L=M;n=0:1:L-1;x=ones(1,L);h=cos(0.2.*pi.*n);disp('L=M=');disp(M)

22、disp('線性卷積用時(shí)：')ticy1=conv(x,h);tocdisp('快速卷積用時(shí)：');ticX1=fft(x);H=fft(h);Y=X1.*H;y=ifft(Y);tocend實(shí)驗(yàn)結(jié)果：L=M= 8線性卷積用時(shí)：時(shí)間已過(guò) 0.000122 秒。快速卷積用時(shí)：時(shí)間已過(guò) 0.000042 秒。L=M= 16線性卷積用時(shí)：時(shí)間已過(guò) 0.000045 秒。快速卷積用時(shí)：時(shí)間已過(guò) 0.000026 秒。L=M= 32線性卷積用時(shí)：時(shí)間已過(guò) 0.000051 秒?？焖倬矸e用時(shí)：時(shí)間已過(guò) 0.000030 秒。L=M= 64線性卷積用時(shí)：時(shí)間已過(guò) 0.0000

23、50 秒?？焖倬矸e用時(shí)：時(shí)間已過(guò) 0.000030 秒。L=M= 256線性卷積用時(shí)：時(shí)間已過(guò) 0.000089 秒。快速卷積用時(shí)：時(shí)間已過(guò) 0.000066 秒。L=M= 512線性卷積用時(shí)：時(shí)間已過(guò) 0.000540 秒?？焖倬矸e用時(shí)：時(shí)間已過(guò) 0.001388 秒。L=M= 1024線性卷積用時(shí)：時(shí)間已過(guò) 0.000617 秒。快速卷積用時(shí)：時(shí)間已過(guò) 0.000115 秒。分析：可見(jiàn)在在相同長(zhǎng)度下，快速卷積比線性卷積快差不多一倍的時(shí)間。實(shí)驗(yàn)2.2實(shí)驗(yàn)代碼：L=2048;M=256;n1=0:1:L;n2=0:1:M;x=ones(1,L);h=cos(0.2.*pi.*n2);disp(

24、'線性卷積：')ticy1=conv(x,h);tocdisp('快速卷積：')N=M+L-1;ticX1=fft(x,N);H=fft(h,N);Y=X1.*H;y=ifft(Y,N);toc實(shí)驗(yàn)結(jié)果： L=2048，M=256線性卷積：時(shí)間已過(guò) 0.000437 秒?？焖倬矸e：時(shí)間已過(guò) 0.000407 秒。L=4096，M=256線性卷積：時(shí)間已過(guò) 0.000445 秒?？焖倬矸e：時(shí)間已過(guò) 0.001395 秒。分析：當(dāng)序列長(zhǎng)度較短時(shí)，快速卷積比線性卷積快；當(dāng)序列長(zhǎng)度過(guò)長(zhǎng)時(shí)，快速卷積比線性卷積慢了許多。實(shí)驗(yàn)2.3實(shí)驗(yàn)代碼：L=2048;M=256;disp

25、('L=');disp(L);disp('M=');disp(M);x=ones(1,L);n=0:1:M-1;h=cos(0.2.*pi.*n);Li=M;Num=L/Li;Ni=Li+M-1;N_yu=L+M-1;y = zeros(1,N_yu);disp('重疊相加法卷積：');ticH = fft(h,Ni);for i=1:Num;i_low=(i-1)*M+1;xi=x(i_low:i_low+M-1);Xi=fft(xi,Ni);Yi=Xi.*H;yi=ifft(Yi,Ni);y(i_low:i_low+Ni-1)=y(i_low

26、:i_low+Ni-1)+yi;endtoc實(shí)驗(yàn)結(jié)果：L= 2048M= 256重疊相加法卷積時(shí)間時(shí)間已過(guò) 0.000626 秒。分析：較第2題結(jié)果快。實(shí)驗(yàn)2.4實(shí)驗(yàn)代碼：L=2048;M=256;disp('L=');disp(L);disp('M=');disp(M);x=ones(1,L);n=0:1:M-1;h=cos(0.2.*pi.*n);Li=M;Num=L/Li;Ni=Li+M-1;N_yu=L+M-1;y=zeros(1,N_yu);x_add=zeros(1,M-1),x;disp('Öصþ&

27、#177;£Áô·¨¾í»ýʱ¼ä');ticH = fft(h,Ni);for i=1:Num;i_low=(i-1)*M+1;xi=x_add(i_low:i_low+M-1);Xi=fft(xi,Ni);Yi=Xi.*H;yi=ifft(Yi,Ni);y(i_low:i_low+M-1)=yi(M:Ni);endxi=x_add(i*M+1):N_yu);Xi=fft(xi,Ni);Y=Xi.*H;yi=ifft(Y,Ni);y(i*M+1):N_y

28、u)=yi(M:Ni-1);toc實(shí)驗(yàn)結(jié)果：L= 2048M= 256重疊保留法卷積時(shí)間時(shí)間已過(guò) 0.000508 秒。分析：較第2題結(jié)果快。五、心得與體會(huì)本次實(shí)驗(yàn)要求我們掌握利用FFT計(jì)算線性卷積的原理及具體實(shí)現(xiàn)方法，以及通過(guò)對(duì)快速卷積和線性卷積運(yùn)算速度的比較，更加直觀的看到利用FFT快速卷積的優(yōu)點(diǎn)。本次實(shí)驗(yàn)讓我切實(shí)看到了FFT算法的高效性及對(duì)于龐大的數(shù)據(jù)的處理實(shí)時(shí)性，同時(shí)對(duì)重疊保留法、重疊相加法也有了更深的認(rèn)識(shí)，對(duì)以后的實(shí)驗(yàn)有了很好的理論基礎(chǔ)，受益頗多。實(shí)驗(yàn)3 IIR數(shù)字濾波器設(shè)計(jì)一、實(shí)驗(yàn)?zāi)康?.掌握利用脈沖響應(yīng)不變法和雙線性變換法設(shè)計(jì)IIR數(shù)字濾波器的原理及具體方法。2.加深理解數(shù)字濾波

29、器和模擬濾波器之間的技術(shù)指標(biāo)轉(zhuǎn)化。3.掌握脈沖響應(yīng)不變法和雙線性變換法設(shè)計(jì)IIR數(shù)字濾波器的優(yōu)缺點(diǎn)及適用范圍。二、實(shí)驗(yàn)內(nèi)容1.設(shè)采樣頻率為fs =4kHz，采用脈沖響應(yīng)不變法設(shè)計(jì)一個(gè)三階巴特沃斯數(shù)字低通濾波器，其3dB截止頻率為fc=1kHz。2.設(shè)采樣頻率為fs=10kHz，設(shè)計(jì)數(shù)字低通濾波器，滿足如下指標(biāo) 通帶截止頻率：fp=1kHz，通帶波動(dòng)：Rp=1dB 阻帶截止頻率：fst=1.5kHz，阻帶衰減：As=15dB要求分別采用巴特沃斯、切比雪夫I型、切比雪夫II型和橢圓模擬原型濾波器及脈沖響應(yīng)不變法進(jìn)行設(shè)計(jì)。結(jié)合實(shí)驗(yàn)結(jié)果，分別討論采用上述設(shè)計(jì)的數(shù)字濾波器是否都能滿足給定指標(biāo)要求，分析脈

30、沖響應(yīng)不變法設(shè)計(jì)IIR數(shù)字濾波器的優(yōu)缺點(diǎn)及適用范圍。四、實(shí)驗(yàn)代碼及實(shí)驗(yàn)結(jié)果巴特沃斯：實(shí)驗(yàn)代碼：fs=10*1000;fp=1*1000;fst=1.5*1000;Rp=1;As=15;Wp = 2*pi*fp/fs;Ws = 2*pi*fst/fs;N = ceil(log10(10(Rp/10)-1)/(10(As/10)-1)/(2*log10(Wp/Ws);Wc = Wp /(10(Rp/10)-1)(1/(2*N);b,a = butter(N,Wc/pi,'s');bz,az = impinvar(b,a,1)w=0:500*pi/500;H,w=freqz(bz,a

31、z); %頻率響應(yīng)subplot(221);plot(w/pi,abs(H);grid on;xlabel('omega(pi)');ylabel('|H(ejomega)|');subplot(222);plot(w/pi,20*log10(abs(H)/max(abs(H);grid on;xlabel('omega(pi)');ylabel('|H(ejomega)|(dB)');subplot(223);plot(w/pi,angle(H)/pi);grid on;xlabel('omega(pi)');y

32、label('Phase of H(ejomega)(pi)');grd=grpdelay(bz,az,w); %群延時(shí)響應(yīng)subplot(224);plot(w/pi,grd);grid on;xlabel('omega(pi)');ylabel('Group delay');實(shí)驗(yàn)結(jié)果：bz = 1.0e-04 * -0.0000 0.0091 0.2033 0.4463 0.1524 0.0051 0az =1.0000 -5.1378 11.0528 -12.7383 8.2921 -2.8898 0.4211切比雪夫型：實(shí)驗(yàn)代碼：fs=1

33、0*1000; fp=1*1000;fst=1.5*1000;Rp=1;As=15;Wp = 2*pi*fp/fs; Ws = 2*pi*fst/fs; e = sqrt(10(Rp/10)-1);A = 10(As/20);N = ceil(acosh( sqrt(A2-1)/e) / ( acosh(Ws/Wp);b,a = cheby1(N,Rp,Wp/pi,'s');bz,az = impinvar(b,a,1)w=0:500*pi/500;H,w=freqz(bz,az); subplot(221);plot(w/pi,abs(H);grid on;xlabel(&#

34、39;omega(pi)');ylabel('|H(ejomega)|');subplot(222);plot(w/pi,20*log10(abs(H)/max(abs(H);grid on;xlabel('omega(pi)');ylabel('|H(ejomega)|(dB)');subplot(223);plot(w/pi,angle(H)/pi);grid on;xlabel('omega(pi)');ylabel('Phase of H(ejomega)(pi)');grd=grpdelay(bz

35、,az,w); subplot(224);plot(w/pi,grd);grid on;xlabel('omega(pi)');ylabel('Group delay');實(shí)驗(yàn)結(jié)果：bz = 1.0e-03 * 0.0000 0.0623 0.2372 0.0567 0az =1.0000 -3.7710 5.3742 -3.4293 0.8265切比雪夫型：實(shí)驗(yàn)代碼：fs=10*1000;fp=1*1000;fst=1.5*1000;Rp=1;As=15;Wp = 2*pi*fp/fs; Ws = 2*pi*fst/fs; e = sqrt(10(Rp/10)

36、-1);A = 10(As/20);N = ceil( acosh( sqrt(A2-1)/e) / ( acosh(Ws/Wp) );b,a = cheby2(N,As,Wp/pi,'s');bz,az = impinvar(b,a,1)w=0:500*pi/500;H,w=freqz(bz,az); subplot(221);plot(w/pi,abs(H);grid on;xlabel('omega(pi)');ylabel('|H(ejomega)|');subplot(222);plot(w/pi,20*log10(abs(H)/max

37、(abs(H);grid on;xlabel('omega(pi)');ylabel('|H(ejomega)|(dB)');subplot(223);plot(w/pi,angle(H)/pi);grid on;xlabel('omega(pi)');ylabel('Phase of H(ejomega)(pi)');grd=grpdelay(bz,az,w); subplot(224);plot(w/pi,grd);grid on;xlabel('omega(pi)');ylabel('Group de

38、lay');實(shí)驗(yàn)結(jié)果：bz = 0.0884 -0.3342 0.5117 -0.3718 0.1076az = 1.0000 -3.4904 4.6022 -2.7148 0.6048橢圓型：實(shí)驗(yàn)代碼：fs=10*1000;fp=1*1000;fst=1.5*1000;Rp=1;As=15;Wp = 2*pi*fp/fs; Ws = 2*pi*fst/fs; e = sqrt(10(Rp/10)-1);A = 10(As/20);N = ceil(acosh( sqrt(A2-1)/e) / ( acosh(Ws/Wp);b,a = ellip(N,Rp,As,Wp/pi,'

39、s');bz,az = impinvar(b,a,1)w=0:500*pi/500;H,w=freqz(bz,az);subplot(221);plot(w/pi,abs(H);grid on;xlabel('omega(pi)');ylabel('|H(ejomega)|');subplot(222);plot(w/pi,20*log10(abs(H)/max(abs(H);grid on;xlabel('omega(pi)');ylabel('|H(ejomega)|(dB)');subplot(223);plot(w

40、/pi,angle(H)/pi);grid on;xlabel('omega(pi)');ylabel('Phase of H(ejomega)(pi)');grd=grpdelay(bz,az,w); subplot(224);plot(w/pi,grd);grid on;xlabel('omega(pi)');ylabel('Group delay');實(shí)驗(yàn)結(jié)果：bz = 0.1470 -0.5636 0.8366 -0.56az =1.0000 -3.7726 5.3933 -3.4604 0.8408分析：脈沖響應(yīng)不變法的

41、優(yōu)點(diǎn)是頻率坐標(biāo)的變換是線性的，因此如果模擬濾波器的頻響是限帶于折疊頻率以內(nèi)的話，通過(guò)變換后的數(shù)字濾波器的頻響可以不失真地反映原響應(yīng)與頻率之間的關(guān)系。但是，其最大的缺點(diǎn)是有頻譜的周期延拓效應(yīng)，會(huì)產(chǎn)生混疊失真。故，脈沖響應(yīng)不變法只能用于限帶的頻響特性，如衰減特性較好的低通，或帶通，而且高頻衰減越大，頻響的混疊效應(yīng)就越小。雙線性變換法的優(yōu)點(diǎn)是消除了脈沖響應(yīng)不變法的固有的混疊失真，但缺點(diǎn)是頻率變換之間的非常嚴(yán)重的非線性，會(huì)使變換后時(shí)域上的圖像產(chǎn)生非常嚴(yán)重的失真。五、心得與體會(huì) 通過(guò)實(shí)驗(yàn)，學(xué)習(xí)了利用脈沖響應(yīng)不變法設(shè)計(jì)IIR數(shù)字濾波器的基本原理為：從時(shí)域響應(yīng)出發(fā)，使數(shù)字濾波器的單位脈沖響應(yīng)h(n)模仿模擬

42、濾波器的單位沖激響應(yīng)ha(t)，h(n)等于ha(t)的取樣值。所以設(shè)計(jì)中，應(yīng)先根據(jù)數(shù)字指標(biāo)轉(zhuǎn)換來(lái)的模擬指標(biāo)，設(shè)計(jì)模擬濾波器，再利用impinvar函數(shù)實(shí)現(xiàn)脈沖響應(yīng)不變法模擬濾波器到數(shù)字濾波器的變換。脈沖響應(yīng)不變法和雙線性法設(shè)計(jì)IIR數(shù)字濾波器設(shè)計(jì)的基礎(chǔ)，掌握了這倆種方法能夠使我們對(duì)數(shù)字濾波器設(shè)計(jì)有更加深刻的認(rèn)識(shí)。實(shí)驗(yàn)4 頻率取樣法設(shè)計(jì)FIR數(shù)字濾波器一、實(shí)驗(yàn)?zāi)康恼莆疹l率取樣法設(shè)計(jì)FIR數(shù)字濾波器的原理及具體方法二、實(shí)驗(yàn)原理1、基本原理 頻率取樣法從頻域出發(fā)，把理想的濾波器等間隔采樣得到,將作為實(shí)際設(shè)計(jì)濾波器的。 （8-1）得到以后可以由來(lái)確定唯一確定濾波器的單位脈沖響應(yīng),可以由求得： （8-

43、2、3）其中為內(nèi)插函數(shù) （8-4）由求得的頻率響應(yīng)將逼近 如果我們?cè)O(shè)計(jì)的是線性相位FIR濾波器，則的幅度和相位滿足線性相位濾波器的約束條件。我們將表示為如下形式： （8-5）當(dāng)為實(shí)數(shù)，則由此得到 （8-6） 即為中心偶對(duì)稱。在利用線性相位條件可知，對(duì)于1型和2型線性相位濾波器 （8-7）對(duì)于3型和4型線性相位濾波器 （8-8）其中，x表示取小于該數(shù)的最大的整數(shù)。設(shè)計(jì)步驟（1）由給定的理想濾波器給出和。（2）由求得（3) 根據(jù)求得或三、實(shí)驗(yàn)內(nèi)容1. 設(shè)計(jì)一個(gè)數(shù)字低通FIR濾波器，其技術(shù)指標(biāo)如下：_p=0.2，R_p=0.25dB_st=0.3，A_s=50dB分別采用矩形窗、漢寧窗、海明窗、布萊

44、克曼窗、凱瑟窗設(shè)計(jì)該濾波器。結(jié)合實(shí)驗(yàn)結(jié)果分別討論上述方法設(shè)計(jì)的數(shù)字濾波器是否符合指標(biāo)。2. 設(shè)計(jì)一個(gè)數(shù)字帶通FIR濾波器，其技術(shù)指標(biāo)如下：3. 采用頻率取樣設(shè)計(jì)法設(shè)計(jì)FIR數(shù)字低通濾波器，滿足以下指標(biāo)_p=0.2，R_p=0.25dB_st=0.3，A_s=50dB（1）取N=20，過(guò)渡帶沒(méi)有樣本。（2）取N=40，過(guò)渡帶有一個(gè)樣本，T=0.39。（3）取N=60，過(guò)渡帶有兩個(gè)樣本，T1=0.5925, T2=0.1099。（4）分別采用上述方法設(shè)計(jì)的數(shù)字濾波器是否都能滿足給定的指標(biāo)要求。4. 采用頻率取樣技術(shù)設(shè)計(jì)下面的高通濾波器_st=0.6，A_s=50dB_p=0.8，R_p=1dB對(duì)于

45、高通濾波器，N必須為奇數(shù)（或1型濾波器）。選擇N=33，過(guò)渡帶有兩個(gè)樣本，過(guò)渡帶樣本最優(yōu)值為T(mén)1=0.1095，T2=0.598。四、實(shí)驗(yàn)代碼及實(shí)驗(yàn)結(jié)果實(shí)驗(yàn)4.1矩形窗：實(shí)驗(yàn)代碼：wp = 0.2 * pi; wst = 0.3 * pi;tr_width = wst - wp;N = ceil(1.8 * pi / tr_width) + 1;n = 0 : N - 1;wc = (wp + wst) / 2;alpha = (N - 1) / 2;hd = (wc / pi) * sinc(wc / pi) * (n - alpha);w_boxcar = boxcar(N)'h

46、= hd .* w_boxcar;subplot(221);stem(n, hd, 'filled');axis tight; xlabel('n'); ylabel('hd(n)');Hr, w1 = zerophase(h);subplot(222); plot(w1/pi, Hr);axis tight; xlabel('omega/pi'); ylabel('H(omega)');subplot(223); stem(n, h, 'filled');axis tight; xlabel(&#

47、39;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');grid on;實(shí)驗(yàn)結(jié)果：分析：阻帶的衰減小于50db，所以，不符合指標(biāo)。漢寧窗：實(shí)驗(yàn)代碼：wp = 0.2 * pi; wst = 0.3 * pi;tr_width = wst - wp;N = ceil(6.2 * pi / tr_width) + 1;n = 0 : N

48、 - 1;wc = (wp + wst) / 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');axis tight; xlabel('n'); ylabel('hd(n)');Hr, w1 = zerophase(h);subplot(222); plot(w1/pi, Hr);axis tight;

49、 xlabel('omega/pi'); ylabel('H(omega)');subplot(223); stem(n, h, 'filled');axis tight; 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');grid on;實(shí)驗(yàn)結(jié)果：分析：阻帶衰

50、減小于50db，故，不符合指標(biāo)。海明窗：實(shí)驗(yàn)代碼：wp = 0.2 * pi; wst = 0.3 * pi;tr_width = wst - wp;N = ceil(6.6 * pi / tr_width) + 1;n = 0 : N - 1;wc = (wp + wst) / 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');ax

51、is tight; xlabel('n'); ylabel('hd(n)');Hr, w1 = zerophase(h);subplot(222); plot(w1/pi, Hr);axis tight; xlabel('omega/pi'); ylabel('H(omega)');subplot(223); stem(n, h, 'filled');axis tight; xlabel('n'); ylabel('h(n)');H, w = freqz(h, 1);subplot(

52、224); plot(w / pi, 20 * log10(abs(H) / max(abs(H);xlabel('omega/pi'); ylabel('dB');grid on;實(shí)驗(yàn)結(jié)果：分析：阻帶衰減大于50db，故，該數(shù)字濾波器符合指標(biāo)。布萊克曼窗：實(shí)驗(yàn)代碼：wp = 0.2 * pi; wst = 0.3 * pi;tr_width = wst - wp;N = ceil(11 * pi / tr_width) + 1;n = 0 : N - 1;wc = (wp + wst) / 2;alpha = (N - 1) / 2;hd = (wc / pi

53、) * sinc(wc / pi) * (n - alpha);w_blackman = blackman(N)'h = hd .* w_blackman;subplot(221);stem(n, hd, 'filled');axis tight; xlabel('n'); ylabel('hd(n)');Hr, w1 = zerophase(h);subplot(222); plot(w1/pi, Hr);axis tight; xlabel('omega/pi'); ylabel('H(omega)')

54、;subplot(223); stem(n, h, 'filled');axis tight; 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');grid on;實(shí)驗(yàn)結(jié)果：分析：阻帶衰減大于50db，故，該數(shù)字濾波器符合要求。凱瑟窗：實(shí)驗(yàn)代碼：wp = 0.2 * pi; wst = 0.3

55、* pi;tr_width = wst - wp;As = 50;N = ceil(As - 7.95) / (2.285 * tr_width) + 1;beta = 0.1102 * (As - 8.7);n = 0 : N - 1;wc = (wp + wst) / 2;alpha = (N - 1) / 2;hd = (wc / pi) * sinc(wc / pi) * (n - alpha);w_kaiser = kaiser(N, beta)'h = hd .* w_kaiser;subplot(221);stem(n, hd, 'filled');axi

56、s tight; xlabel('n'); ylabel('hd(n)');Hr, w1 = zerophase(h);subplot(222); plot(w1/pi, Hr);axis tight; xlabel('omega/pi'); ylabel('H(omega)');subplot(223); stem(n, h, 'filled');axis tight; xlabel('n'); ylabel('h(n)');H, w = freqz(h, 1);subplot(2

57、24); plot(w / pi, 20 * log10(abs(H) / max(abs(H);xlabel('omega/pi'); ylabel('dB');grid on;實(shí)驗(yàn)結(jié)果：分析：阻帶衰減大于50db，故，該濾波器符合指標(biāo)。實(shí)驗(yàn)4.2實(shí)驗(yàn)代碼：ws1 = 0.2 * pi; ws2 = 0.8 * pi;wp1 = 0.35 * pi; wp2 = 0.65 * pi;tr_width = wp1 - ws1;wc1 = (ws1 + wp1) / 2;wc2 = (ws2 + wp2) / 2;N = ceil(11 * pi / tr_width) + 1;n = 0 : N - 1;alpha = (N - 1) / 2