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

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

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

4、頻譜，得到的結(jié)果只是N個頻譜樣本值，樣本值之間的頻譜是未知的，就像通過一個柵欄觀察頻譜，稱為“柵欄效應(yīng)”。頻譜分辨率與記錄長度成正比，提高頻譜分辨率，就要增加記錄時間。（3）頻譜泄露。對于信號截?cái)鄷汛昂瘮?shù)的頻譜會引入到信號頻譜中，造成頻譜泄露。解決這問題的主要辦法是采用旁瓣小的窗函數(shù)，頻譜泄露和窗函數(shù)均會引起誤差。 因此，要合理選取采樣間隔和截取長度，必要時還需考慮適當(dāng)?shù)拇啊?對于連續(xù)周期信號，我們在采用計(jì)算機(jī)進(jìn)行計(jì)算時，也總是要進(jìn)行截?cái)啵蛄锌偸怯邢揲L的，仍然可以采用上訴方法近似計(jì)算。4、可能用到MATLAB函數(shù)與代碼實(shí)驗(yàn)中的DFT運(yùn)算可以采用MATLAB中提供的FFT來實(shí)現(xiàn)。DTFT可以

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

6、到的頻譜的頻率分辨率符合需要。4.利用DFT近似分析連續(xù)時間信號xt=e-0.1tu(t)的頻譜（幅度譜）。分析采用不同的采樣間隔和截取長度進(jìn)行計(jì)算的結(jié)果，并最終確定合適的參數(shù)。四、實(shí)驗(yàn)代碼、實(shí)驗(yàn)結(jié)果及分析1.（1）>> 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);>> xlabel('Omega/pi');>> title('

7、Magnitude');>> axis tight;>> subplot(212);>> plot(w,angle(X);>> xlabel('Omega/pi');>> title('Phase');>> axis tight;（2）>> Xk=fft(x);>> subplot(211);>> hold on;>> stem(2*pi/4*n,abs(Xk),'filled');>> subplot(21

8、2);>> hold on;>> stem(2*pi/4*n,angle(Xk),'filled');（3）>> x=2 -1 1 1 zeros(1,60);>> n=0:63;>> Xk=fft(x);>> subplot(211);>> stem(n,abs(Xk),'filled');>> xlabel('Omega/pi');>> title('Magnitude');>> axis tight;>

9、;> subplot(212);>> stem(n,angle(Xk),'filled');>> xlabel('Omega/pi');>> title('Phase');>> axis tight;（4） 可以由DFT計(jì)算DTFT，只需采樣點(diǎn)足夠即可逼近DTFT，可采用增加補(bǔ)零個數(shù)的方法實(shí)現(xiàn)。2.（1）>> n=0:10;>> x=cos(0.48*pi*n)+cos(0.52*pi*n);>> Xk=fft(x);>> subplot(211

10、);>> stem(n,abs(Xk),'filled');>> xlabel('Omega');>> title('Magnitude');>> n1=0:99;>> x1=x zeros(1,89);>> Xk1=fft(x1);>> subplot(212);>> stem(n,abs(Xk1),'filled');>> xlabel('Omega');>> title('Magni

11、tude');>> n=0:10;>> x=cos(0.48*pi*n)+cos(0.52*pi*n);>> Xk1=fft(x);>> subplot(211);>> stem(n,angle(Xk),'filled');>> xlabel('Omega');>> title('Phase');>> n1=0:99;>> x1=x zeros(1,89);>> Xk1=fft(x1);>> subplot(

12、212);>> stem(n1,angle(Xk1),'filled');>> xlabel('Omega');>> title('Phase');（2）（3） 可通過增大采樣點(diǎn)數(shù)，時域補(bǔ)零提高頻譜分辨率。3.>> n=0:20;>> x=0.15*sin(2*pi*n)+sin(2*pi*2*n)-0.1*sin(2*pi*3*n);>> Xk=fft(x);>> subplot(211);>> stem(n,abs(Xk),'filled&

13、#39;);>> xlabel('Omega');>> title('Magnitude');>> subplot(212);>> stem(n,angle(Xk),'filled');>> xlabel('Omega');>> title('Phase');4.【修改采樣間隔】>> n=0:10;>> x=exp(-0.1*n).*heaviside(n);>> Xk=fft(x);>> ste

14、m(n,abs(Xk),'filled');>> n=0:50;>> x=exp(-0.1*n).*heaviside(n);>> Xk=fft(x);>> stem(n,abs(Xk),'filled');>> n=0:100;>> x=exp(-0.1*n).*heaviside(n);>> Xk=fft(x);>> stem(n,abs(Xk),'filled');【修改截取長度】>> n=0:100;>> x=exp(-

15、0.1*n).*(heaviside(n)-heaviside(n-50);>> Xk=fft(x);>> stem(n,abs(Xk),'filled');>> n=0:100;>> x=exp(-0.1*n).*(heaviside(n)-heaviside(n-10);>> Xk=fft(x);>> stem(n,abs(Xk),'filled');五、實(shí)驗(yàn)總結(jié)第一次數(shù)字信號處理實(shí)驗(yàn)，我了解了如何在MATLAB中使用DFT分析信號的頻譜，對DFT采樣、截取的概念和操作結(jié)果有了更直觀的認(rèn)

16、識。此外，使用MATLAB的FFT大大簡化了DFT的計(jì)算，讓我認(rèn)識到了FFT算法的快捷。實(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點(diǎn)序列，h(n)為M點(diǎn)序列，x(n)和h(n)的線性卷積為 （3-1）的長度為L+M-1x(n)和h(n)的圓周卷積為 （3-2）圓周卷積與線性卷積相等而不產(chǎn)生交疊的必要條件為NL+M+1 (3-3)圓周卷積定理：根據(jù)DFT性質(zhì)，x(n)和h(n)的N點(diǎn)圓周卷積的DFT等于它們的DF

17、T的乘積： (3-4)2.快速卷積快速卷積發(fā)運(yùn)用圓周卷積實(shí)現(xiàn)線性卷積，根據(jù)圓周卷積定理利用FFT算法實(shí)現(xiàn)圓周卷積。可將快速卷積運(yùn)算的步驟歸納如下：(1)必須選擇；為了能使用基-2算法，要求。采用補(bǔ)零的辦法使得x(n)和h(n)的長度均為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)的線性系統(tǒng)，輸入為x(n)，輸出為y(n)，則yn=xn*h(n)當(dāng)輸入序列x(n)極長時，如果要等x(n)全部集齊時再開始進(jìn)行卷積，會使輸出有較大延時；如果序列太長，需要大量存儲單元。為此，我們把x

18、(n)分段，為別求出每段的卷積，合在一起得到最后的總輸出。這稱為分段卷積。分段卷積可以細(xì)分為重疊保留法和重疊相加法。重疊保留法：設(shè)x(n)的長度為，h(n)的長度為M。把序列x(n)分成多段N點(diǎn)序列，每段雨前一段重寫M-1個樣本。并在第一個輸入段前面補(bǔ)M-1個零。計(jì)算每一段與h(n)的圓周卷積，其結(jié)果中前M-1個不等與線性卷積，應(yīng)當(dāng)舍去，只保留后面N-M+1個正確的輸出樣本，把它們合起來得到總的輸出。利用FFT實(shí)現(xiàn)重疊保留法的步驟如下：(1)在x(n)前面填充M-1個零，擴(kuò)大以后的序列為(2)將x(n)分為若干段N點(diǎn)子段，設(shè)L=N-M+1為每一段的有效長度，則第i段的數(shù)據(jù)為：(3)計(jì)算每一段與

19、h(n)的N點(diǎn)圓周卷積，利用FFT計(jì)算圓周卷積(4)舍去每一段卷積結(jié)果的前M-1個樣本，連接剩下的樣本得到卷積結(jié)果y(n)。重疊相加法：設(shè)h(n)長度為M，將信號x(n)分解成長為L的子段。以表示沒斷信號，則：每一段卷積的長度為L+M-1，所以在做求和時，相鄰兩段序列由M-1個樣本重疊，即前一段的最后M-1個樣本和下一段前M-1個樣本序列重疊，這個重疊部分相加，再與不重疊的部分共同組成y(n)。利用FFT實(shí)現(xiàn)重疊保留法的步驟如下：(1)將x(n)分為若干L點(diǎn)子段。(2)計(jì)算每一段與h(n)的卷積，根據(jù)快速卷積法利用FFT計(jì)算卷積。(3)將各段相加，得到輸出y(n)。4、可能得到的MATLAB函

20、數(shù) 實(shí)驗(yàn)中FFT運(yùn)算可采用MATLAB中提供的函數(shù)fft來實(shí)現(xiàn)。三、實(shí)驗(yàn)題目假設(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ì)算兩個序列線性卷積的方法,比較當(dāng)序列長度分別為8,16,32,64,256,512,1024時兩種方法計(jì)算線性卷積所需時間。2當(dāng)L=2048且M=256時比較直接計(jì)算線性卷積和快速卷積所需的時間,進(jìn)一步考察當(dāng)L=4096且M=256時兩種算法所需的時間。3. 編程實(shí)現(xiàn)利用重疊相加法計(jì)算兩個序列的線性卷積,考察L=2048且M=256時計(jì)算

21、線性卷積的時間,與2題的結(jié)果進(jìn)行比較。4. 編程實(shí)現(xiàn)利用重疊保留法計(jì)算兩個序列的線性卷積,考察L=2048且M=256時計(jì)算線性卷積的時間,與2題的結(jié)果進(jìn)行比較。四、實(shí)驗(yàn)代碼、實(shí)驗(yàn)結(jié)果及分析1.【線性卷積】>> for i=1:7L=input('L:');n=0:L;x=heaviside(n)-heaviside(n-L);h=cos(0.2*pi*n);ticy=conv(x,h);tocendL:8時間已過 0.077287 秒。L:16時間已過 0.032760 秒。L:32時間已過 0.000095 秒。L:64時間已過 0.000060 秒。L:256

22、時間已過 0.000131 秒。L:512時間已過 0.000230 秒。L:1024時間已過 0.000846 秒。【快速卷積】>> for i=1:7L=input('L:');n=0:L;x=heaviside(n)-heaviside(n-L);h=cos(0.2*pi*n);ticXk=fft(x,L+1);Hk=fft(h,L+1);Yk=Xk.*Hk;y=ifft(Yk);tocendL:8時間已過 0.002417 秒。L:16時間已過 0.001108 秒。L:32時間已過 0.000384 秒。L:64時間已過 0.000616 秒。L:256時

23、間已過 0.002367 秒。L:512時間已過 0.000702 秒。L:1024時間已過 0.000967 秒。在序列長度較短時，快速卷積優(yōu)勢明顯，序列加長后出現(xiàn)了比線性卷積慢的情況2.線性卷積：>> for i=1:2L=input('L:');n1=0:L;M=input('M:');n2=0:M;x=heaviside(n1)-heaviside(n1-L);h=cos(0.2*pi*n2);ticy=conv(x,h);tocendL:2048M:256時間已過 0.000533 秒。L:4096M:256時間已過 0.000712 秒。

24、快速卷積：for i=1:2L=input('L:');n1=0:L;M=input('M:');n2=0:M;x=heaviside(n1)-heaviside(n1-L);h=cos(0.2*pi*n2);ticXk=fft(x,M);Hk=fft(h,M);Yk=Xk.*Hk;y=ifft(Yk);tocendL:2048M:256時間已過 0.000592 秒。L:4096M:256時間已過 0.001827 秒。可見，當(dāng)序列加長，快速卷積速度開始優(yōu)于線性卷積3.>> tic>> N=512;>> m=0:256;&g

25、t;> h=cos(0.2*pi*m);>> n=0:2048;>> x=heaviside(n)-heaviside(n-2048);>> Lx=length(x);>> Lm=length(h);>> M=Lm-1;>> L=N-M;>> h=fft(h,N);>> K=ceil(Lx/L);>> for i=Lx:K*L-1x(i+1)=0;end>> Y=zeros(K,N);>> Y2=zeros(1,(K-1)*L+N);>> for

26、 k=0:K-1Xk=x(k*L+1:k*L+L),zeros(1,M);Y(k+1, :)=(ifft(fft(Xk).*h);Y2(k*L+1:k*L+N)=Y2(k*L+1:k*L+N)+Y(k+1,:);end>> toc時間已過 0.000575 秒。可見，重疊相加法的速度相比第二題更快4.>> tic>> N=512;>> m=0:256;>> h=cos(0.2*pi*m);>> n=0:2048;>> x=heaviside(n)-heaviside(n-2048);>> Lx=le

27、ngth(x);>> Lm=length(h);>> M=Lm-1;>> L=N-M;>> h=fft(h,N);>> K=floor(Lx+M-1)/L)+1;>> z=(K)*L-Lx;>> x1=zeros(1,M),x,zeros(1,z);>> Y=zeros(K,N);>> for k=0:K-1Xk=fft(x1(k*L+1:k*L+N);Y(k+1, :)=(ifft(Xk).*h);end>> Z=reshape(Y(:,M:N)',1,);>

28、> toc時間已過 0.000452 秒。可見，重疊保留法的速度相比第二題也更快五、實(shí)驗(yàn)總結(jié)這次FFT的實(shí)驗(yàn)讓我對FFT的算法效率有了新的認(rèn)識，并不是任何情況下FFT卷積的速度都優(yōu)于直接計(jì)算線性卷積，當(dāng)序列長度很長時反而會比線性卷積慢。此外，我還對兩種新的FFT方法：重疊相加法、重疊保留法及其MATLAB實(shí)現(xiàn)有了基本的了解。實(shí)驗(yàn)3 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)題目1.設(shè)采

29、樣頻率為fs =4kHz，采用脈沖響應(yīng)不變法設(shè)計(jì)一個三階巴特沃斯數(shù)字低通濾波器，其3dB截止頻率為fc=1kHz。2.設(shè)采樣頻率為fs=10kHz，設(shè)計(jì)數(shù)字低通濾波器，滿足如下指標(biāo) 通帶截止頻率：fp=1kHz，通帶波動：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)要求，分析脈沖響應(yīng)不變法設(shè)計(jì)IIR數(shù)字濾波器的優(yōu)缺點(diǎn)及適用范圍。三、實(shí)驗(yàn)代碼、實(shí)驗(yàn)結(jié)果及分析1. 巴特沃斯模擬原型濾波器 + 脈沖響應(yīng)不變法

30、/ 雙線性變換法>> fp=1000;>> fs=1500;>> Rp=1;>> As=15;>> f=10000;>> Wp=2*pi*fp;>> Ws=2*pi*fs;>> wp=2*f*tan(2*pi*fp/(2*f);>> ws=2*f*tan(2*pi*fs/(2*f);>> N1=ceil(log10(10(Rp/10)-1)/(10(As/10)-1)/(2*log10(Wp/Ws);>> N2=ceil(log10(10(Rp/10)-1)/(1

31、0(As/10)-1)/(2*log10(wp/ws);>> OmegaC1=Wp/(10(Rp/10)-1)(1/(2*N1);>> OmegaC2=wp/(10(Rp/10)-1)(1/(2*N1);>> z,p,k=buttap(N1);>> b1=k*OmegaC1N1;>> a1=poly(p*OmegaC1);>> H,w=freqs(b1,a1);>> subplot(221);>> plot(w/pi,abs(H);>> grid on;>> axis tig

32、ht;>> xlabel('Omega(pi)');>> ylabel('|H(Omega)|');>> subplot(223);>> plot(w/pi,angle(H);>> grid on;>> axis tight;>> xlabel('Omega(pi)');>> ylabel('Phase of |H(Omega)|(pi)');>> b2,a2=butter(N2,OmegaC2,'S');&

33、gt;> H,w=freqs(b2,a2);>> subplot(222);>> plot(w/pi,abs(H);>> grid on;>> axis tight;>> xlabel('Omega(pi)');>> ylabel('|H(Omega)|');>> subplot(224);>> plot(w/pi,angle(H);>> grid on;>> axis tight;>> xlabel('Omega(p

34、i)');>> ylabel('Phase of |H(Omega)|(pi)');2. 切比雪夫I型 + 脈沖響應(yīng)不變法 / 雙線性變換法>> wp=0.2*pi;>> ws=0.3*pi;>> Rp=1;>> As=15;>> f=10000;>> T=1/f;>> e=(10(Rp/10)-1)(1/2);>> A=10(As/20);>> wp1=wp/T;>> ws1=ws/T;>> N=ceil(acosh(A2-1

35、)(1/2)/e)/(acosh(ws1/wp1);>> wc=wp/pi;>> b,a=cheby1(N,Rp,wc);>> w=0:500*pi/500;>> H,w=freqz(b,a);>> subplot(221);>> plot(w/pi,abs(H);>> grid on;>> axis tight;>> xlabel('omega(pi)');>> ylabel('|H(ejomega)|');>> subplot(

36、222);>> plot(w/pi,20*log10(abs(H)/max(abs(H);>> grid on;>> axis tight;>> xlabel('omega(pi)');>> ylabel('|H(ejomega)|(dB)');>> subplot(223);>> plot(w/pi,angle(H)/pi);>> grid on;>> axis tight;>> xlabel('omega(pi)');>

37、;> ylabel('Phase of H(ejomega)');>> grd=grpdelay(b,a,w);>> subplot(224);>> plot(w/pi,grd);>> grid on;>> axis tight;>> xlabel('omega(pi)');>> ylabel('Group Delay');3. 切比雪夫II型 + 脈沖響應(yīng)不變法 / 雙線性變換法>> wp=0.2*pi;>> ws=0.3*pi;&

38、gt;> Rp=1;>> As=15;>> f=10000;>> T=1/f;>> e=(10(Rp/10)-1)(1/2);>> A=10(As/20);>> wp1=wp/T;>> ws1=ws/T;>> N=ceil(acosh(A2-1)(1/2)/e)/(acosh(ws1/wp1);>> wc=wp/pi;>> b,a=cheby2(N,Rp,wc);>> w=0:500*pi/500;>> H,w=freqz(b,a);>&g

39、t; subplot(221);>> plot(w/pi,abs(H);>> grid on;>> axis tight;>> xlabel('omega(pi)');>> ylabel('|H(ejomega)|');>> subplot(222);>> plot(w/pi,20*log10(abs(H)/max(abs(H);>> grid on;>> axis tight;>> xlabel('omega(pi)');&g

40、t;> ylabel('|H(ejomega)|(dB)');>> subplot(223);>> plot(w/pi,angle(H)/pi);>> grid on;>> axis tight;>> xlabel('omega(pi)');>> ylabel('Phase of H(ejomega)');>> grd=grpdelay(b,a,w);>> subplot(224);>> plot(w/pi,grd);>>

41、 grid on;>> axis tight;>> xlabel('omega(pi)');>> ylabel('Group Delay');4. 橢圓型 + 脈沖響應(yīng)不變法 / 雙線性變換法>> wp=0.2*pi;>> ws=0.3*pi;>> Rp=1;>> As=15;>> f=10000;>> T=1/f;>> e=(10(Rp/10)-1)(1/2);>> A=10(As/20);>> wp1=wp/T;&g

42、t;> ws1=ws/T;>> N=ceil(acosh(A2-1)(1/2)/e)/(acosh(ws1/wp1);>> wc=wp/pi;>> b,a=ellip(N,Rp,wc);>> w=0:500*pi/500;>> H,w=freqz(b,a);>> subplot(221);>> plot(w/pi,abs(H);>> grid on;>> axis tight;>> xlabel('omega(pi)');>> ylabel(

43、'|H(ejomega)|');>> subplot(222);>> plot(w/pi,20*log10(abs(H)/max(abs(H);>> grid on;>> axis tight;>> xlabel('omega(pi)');>> ylabel('|H(ejomega)|(dB)');>> subplot(223);>> plot(w/pi,angle(H)/pi);>> grid on;>> axis tight

44、;>> xlabel('omega(pi)');>> ylabel('Phase of H(ejomega)');>> grd=grpdelay(b,a,w);>> subplot(224);>> plot(w/pi,grd);>> grid on;>> axis tight;>> xlabel('omega(pi)');>> ylabel('Group Delay');五、實(shí)驗(yàn)總結(jié)通過這次實(shí)驗(yàn)，我了解了巴特沃斯、切比雪夫I

45、/II型、橢圓型濾波器的MATLAB設(shè)計(jì)方法，對課本上介紹的濾波器設(shè)計(jì)步驟有了實(shí)踐，對其思路有了更深的認(rèn)識，譬如從數(shù)字域指標(biāo)出發(fā)設(shè)計(jì)模擬域的指標(biāo)，求出模擬濾波器的系統(tǒng)函數(shù)再轉(zhuǎn)換回?cái)?shù)字濾波器。此外，MATLAB提供的濾波器設(shè)計(jì)函數(shù)大大簡化了設(shè)計(jì)，十分便捷。實(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）得到以后可以由來確定唯一確定濾波器的單位脈沖響應(yīng),可以由求得： （8-2、3）其中為內(nèi)插函數(shù) （8-4）由求得的頻率響應(yīng)將逼近 如果

46、我們設(shè)計(jì)的是線性相位FIR濾波器，則的幅度和相位滿足線性相位濾波器的約束條件。我們將表示為如下形式： （8-5）當(dāng)為實(shí)數(shù)，則由此得到 （8-6） 即為中心偶對稱。在利用線性相位條件可知，對于1型和2型線性相位濾波器 （8-7）對于3型和4型線性相位濾波器 （8-8）其中，x表示取小于該數(shù)的最大的整數(shù)。設(shè)計(jì)步驟（1）由給定的理想濾波器給出和。（2）由求得（3) 根據(jù)求得或三、實(shí)驗(yàn)題目1. 設(shè)計(jì)一個數(shù)字低通FIR濾波器，其技術(shù)指標(biāo)如下：_p=0.2，R_p=0.25dB_st=0.3，A_s=50dB分別采用矩形窗、漢寧窗、海明窗、布萊克曼窗、凱瑟窗設(shè)計(jì)該濾波器。結(jié)合實(shí)驗(yàn)結(jié)果分別討論上述方法設(shè)計(jì)的

47、數(shù)字濾波器是否符合指標(biāo)。2. 設(shè)計(jì)一個數(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，過渡帶沒有樣本。（2）取N=40，過渡帶有一個樣本，T=0.39。（3）取N=60，過渡帶有兩個樣本，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對于高通濾波器，N必須為奇數(shù)（或1型濾波器）。選擇N=33，過渡帶

48、有兩個樣本，過渡帶樣本最優(yōu)值為T1=0.1095，T2=0.598。四、實(shí)驗(yàn)代碼、實(shí)驗(yàn)結(jié)果及分析1.【矩形窗】>> 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=hd .*

49、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

50、, 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;顯然，矩形窗組帶衰減不滿足指標(biāo)【漢寧窗】>> wp=0.2*pi;>> wst=0.

51、3*pi;>> tr_width=wst-wp;>> N=ceil(6.2*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_hanning= hanning (N)'>> h=hd .*w_hanning;>> subplot(221);stem(n, hd,'filled');>> axis tight; xla

52、bel('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)'

53、;);>> H,w=freqz(h,1);>> subplot(224); plot(w/pi,20*log10(abs(H)/max(abs(H);>> xlabel('omega/pi'); ylabel('dB');>> grid on;顯然，漢寧窗組帶衰減不滿足指標(biāo)【海明窗】>> 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

54、=(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');>> axis tight; xlabel('n'); ylabel('hd(n)');>> Hr,w1 = zerophase(h);>> sub

55、plot(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(224); plot(w/pi,20*log10(abs(H)/max(abs(

56、H);>> xlabel('omega/pi'); ylabel('dB');>> grid on;可見，海明窗滿足指標(biāo)要求【布萊克曼窗】>> 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)*sinc(wc/pi)*(n-alpha);>>

57、; 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(

58、'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;可見，布萊

59、克曼窗滿足指標(biāo)要求【凱瑟窗】>> wp=0.2*pi;>> wst=0.3*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)'

60、>> h=hd .*w_kaiser;>> 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)');>> subp

61、lot(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;可見，凱瑟窗滿足指標(biāo)要求2.【不妨使用布萊克曼窗設(shè)計(jì)】>> wp

62、1=0.35*pi;>> wp2=0.65*pi;>> wst1=0.2*pi;>> wst2=0.8*pi;>> tr_width=wp1-wst1;>> N=ceil(11*pi/tr_width)+1;>> n=0:N-1; wc1=(wp1+wst1)/2; wc2=(wp2+wst2)/2;>> alpha=(N-1)/2;>> hd=(wc2/pi)*sinc(wc2/pi)*(n-alpha)-(wc1/pi)*sinc(wc1/pi)*(n-alpha);>> w_bla

63、ckman=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(ome

64、ga)');>> 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è)計(jì)指標(biāo)要求3. （1）>> 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*