數(shù)字信號(hào)處理課程設(shè)計(jì)---正余弦信號(hào)的譜分析

1、<p>　　《數(shù)字信號(hào)處理》課程設(shè)計(jì)</p><p>　　題 目： 正余弦信號(hào)的譜分析 </p><p>　　姓 名： </p><p>　　院 系： 電子信息工程系 </p><p>　　專 業(yè)： 通信工程

2、 </p><p>　　班 級(jí)： 通信091 </p><p>　　學(xué) 號(hào)： </p><p>　　指導(dǎo)教師： </p><p><b>　　2012年 6 月</b></p>&l

3、t;p><b>　　正余弦信號(hào)的譜分析</b></p><p>　　中文摘要：使用MATLAB軟件，通過(guò)編寫程序，對(duì)正余弦信號(hào)進(jìn)行傅里葉變換。用DFT和FFT實(shí)現(xiàn)對(duì)正余弦信號(hào)的譜分析，并且分析DFT長(zhǎng)度對(duì)頻譜的影響。</p><p>　　關(guān)鍵詞：matlab； 傅里葉變換； DFT; FFT; </p><p><b>　　一、

4、概述</b></p><p>　　數(shù)字信號(hào)處理方法的一個(gè)重要用途是在離散時(shí)間域中確定一個(gè)連續(xù)時(shí)間信號(hào)的頻譜，通常稱為頻譜分析，更具體的說(shuō)，它也包括確定能量譜或功率譜。數(shù)字頻譜分析可以應(yīng)用在很廣闊的領(lǐng)域。</p><p><b>　　二、設(shè)計(jì)目的</b></p><p>　　1.用DFT實(shí)現(xiàn)對(duì)正余弦信號(hào)的譜分析；</p>

5、<p>　　2.觀察DFT長(zhǎng)度和窗函數(shù)長(zhǎng)度對(duì)頻譜的影響；</p><p>　　3.對(duì)DFT進(jìn)行譜分析中的誤差現(xiàn)象獲得感性認(rèn)識(shí)。</p><p><b>　　三、設(shè)計(jì)原理</b></p><p><b>　　1、譜分析原理</b></p><p>　　頻譜分析方法是基于以下的觀測(cè)：如果連續(xù)

6、時(shí)間信號(hào)是頻帶有限的，那么對(duì)其離散時(shí)間等效信號(hào)的DFT進(jìn)行譜分析。它的離散時(shí)間等效物g(n)應(yīng)當(dāng)能給出(t)頻譜的一個(gè)很近似的估計(jì)兩者之間只差一個(gè)帶數(shù)因子T。然而，在大多數(shù)情況下，是在范圍內(nèi)定義的，因此也就定義在的無(wú)線范圍內(nèi)，要估計(jì)一個(gè)無(wú)限長(zhǎng)信號(hào)的頻譜是不可能的。實(shí)用的方法是：先讓模擬連續(xù)信號(hào)通過(guò)一個(gè)抗混疊的模擬濾波器，然后把它采樣成一個(gè)離散序列。假定反混疊濾波器的設(shè)計(jì)是正確的，則混疊效應(yīng)可以忽略，又假設(shè)A/D變換器的字長(zhǎng)足夠長(zhǎng)，則A/

7、D變換中的量化噪聲也可忽略。</p><p>　　假定表征正余弦信號(hào)的基本參數(shù)，如振幅、頻率和相位不隨時(shí)間改變，則此信號(hào)的傅立葉變換可以用計(jì)算它的DTFT得到</p><p><b>　?。?）</b></p><p>　　實(shí)際上無(wú)限長(zhǎng)序列首先乘以一個(gè)長(zhǎng)度為M的窗函數(shù)，使它變成一個(gè)長(zhǎng)為M的有限長(zhǎng)序列，，對(duì)求出的DTFT 應(yīng)該可以作為原連續(xù)模擬信

8、號(hào)的頻譜估計(jì)，然后求出在區(qū)間等分為N點(diǎn)的離散傅立葉變換DFT。為保證足夠的分辨率，DFT的長(zhǎng)度N選的比窗長(zhǎng)度M大，其方法是在截?cái)嗔说男蛄泻竺嫜a(bǔ)上N－M個(gè)零。計(jì)算采用FFT算法。</p><p>　　更詳細(xì)地考察一下上面的方法。這樣才能了解它的限制，并正確利用它所得出的結(jié)果。特別要分析加窗的效果，以及和由DFT樣本來(lái)估計(jì)DTFT頻率采樣值的問(wèn)題。</p><p>　　在討論由來(lái)估計(jì)頻譜和時(shí)，

9、需要重新探討一下這些變換和它們所對(duì)應(yīng)的頻率之間的關(guān)系，R點(diǎn)的DFT與它的DTFT的關(guān)系為：</p><p>　　=| (2)</p><p>　　歸一化的數(shù)字角頻率和DFT樣本序號(hào)k的關(guān)系為</p><p>　　= (3)<

10、/p><p>　　同樣的，模擬角頻率和DFT樣本序號(hào)k的關(guān)系為</p><p><b>　　(4)</b></p><p>　　其中T是采樣周期。為了能正確地解釋一DFT為基礎(chǔ)的頻譜分析結(jié)果，首先來(lái)考慮單頻率正余弦序列的頻譜分析。設(shè)一個(gè)具有數(shù)字角頻率的余弦信號(hào)為：</p><p><b>　　(5)</b&g

11、t;</p><p><b>　　把這個(gè)序列表為：</b></p><p><b>　　(6)</b></p><p>　　查表得知它的DTFT為：</p><p><b>　　(7)</b></p><p>　　因此，它是一個(gè)以2為周期的的周期信號(hào)，每

12、個(gè)周期中包含兩個(gè)沖擊信號(hào)。在的頻率范圍內(nèi)，處的沖擊具有復(fù)數(shù)幅特性，而在處具有的復(fù)數(shù)幅特性。</p><p>　　為了用DFT分析個(gè)g(n)的頻域特性，取g(n)的一個(gè)有限長(zhǎng)序列</p><p><b>　　(8)</b></p><p>　　2、MATLAB函數(shù)介紹</p><p>　　1. 輸入函數(shù)input( )&l

13、t;/p><p>　　格式：R=input(string)</p><p>　　功能：在屏幕上顯示input括號(hào)后的’string’內(nèi)容，提示用戶從鍵盤輸入某值，并將輸入的值賦給R。</p><p>　　例如，在命令窗口輸入R=input(‘How many apples’)</p><p>　　會(huì)顯示How many apples</p&

14、gt;<p><b>　　從鍵盤輸入 3</b></p><p><b>　　會(huì)顯示 R=3</b></p><p>　　2. 一維快速傅里葉變換函數(shù)fft( )</p><p>　　格式： y=fft(x)</p><p>　　y=fft(x,n)</p><p&

15、gt;　　說(shuō)明：fft函數(shù)用于計(jì)算矢量或矩陣的傅里葉變換。</p><p>　　３．幅度函數(shù)abs(　)</p><p><b>　　格式：abs（X）</b></p><p>　　功能：對(duì)X取絕對(duì)值，當(dāng)X是復(fù)數(shù)時(shí)，得到X的復(fù)模值。</p><p><b>　　四、設(shè)計(jì)過(guò)程</b></p>

16、;<p>　　1.理解題目，編寫程序；</p><p><b>　　2.運(yùn)行程序；</b></p><p><b>　　3.輸入相關(guān)參數(shù)；</b></p><p><b>　　4.觀察圖像。</b></p><p><b>　　設(shè)計(jì)過(guò)程框圖：</b

17、></p><p><b>　　五、設(shè)計(jì)內(nèi)容</b></p><p>　　1.對(duì)一個(gè)頻率為10Hz，采樣頻率為64Hz的32點(diǎn)余弦序列進(jìn)行譜分析，畫(huà)出其頻譜圖；若將頻率改為11Hz，其他參數(shù)不變，重新畫(huà)出該序列的頻譜圖，觀察頻率泄漏現(xiàn)象，分析原因。</p><p>　?。?）.頻率為10Hz，采樣頻率為64Hz的32點(diǎn)余弦序列</p

18、><p>　　clear all;close all;clc;</p><p>　　N=32;n=0:N-1;Fs=64;T=1/Fs;</p><p>　　x1=cos(10*2*pi*n*T); %定義10Hz序列x1</p><p><b>　　k=0:N-1;</b></p>&

19、lt;p>　　X1=abs(fft(x1,N)); %求余弦序列x1的32點(diǎn)FFT</p><p>　　subplot(1,2,1);stem(n,x1); %繪制10Hz序列x1的波形</p><p>　　xlabel('n');ylabel('x1(n)'); title('余弦序列(f=10

20、Hz)');</p><p>　　subplot(1,2,2);stem(k,X1); %繪制序列x1的幅頻特性曲線</p><p>　　xlabel('k');ylabel('X1(k)'); title('32點(diǎn)FFT幅頻曲線(f=10Hz)');</p><p>　　圖1. 頻率為10H

21、z，采樣頻率為64Hz的32點(diǎn)余弦序列頻譜圖</p><p>　?。?）頻率為11Hz，采樣頻率為64Hz的32點(diǎn)余弦序列</p><p>　　clear all;close all;clc;</p><p>　　N=32;n=0:N-1;Fs=64;T=1/Fs;</p><p>　　x2=cos(11*2*pi*n*T);

22、 %定義11Hz序列x2</p><p><b>　　k=0:N-1;</b></p><p>　　X2=abs(fft(x2,N)); %求余弦序列x2的32點(diǎn)FFT</p><p>　　subplot(1,2,1);stem(n,x2); %繪制11Hz序列x2的波形</p>

23、;<p>　　xlabel('n');ylabel('x2(n)'); title('余弦序列(f=11Hz)');</p><p>　　subplot(1,2,2);stem(k,X2); %繪制序列x2的幅頻特性曲線</p><p>　　xlabel('k');ylabel('X2

24、(k)');title('32點(diǎn)FFT幅頻曲線(f=11Hz)');</p><p>　　圖2. 頻率為11Hz，采樣頻率為64Hz的32點(diǎn)余弦序列頻譜圖</p><p><b>　　分析：</b></p><p>　　由圖1可見(jiàn)，它的DFT只有兩個(gè)點(diǎn)不等于零，這樣DFT確實(shí)正確地分辨了余弦信號(hào)的頻率。但是這樣理想的結(jié)果

25、是碰巧得到的，因?yàn)槲覀兦『媒厝×宋鍌€(gè)完整的余弦周期(f*N/Fs=5)。如果截取的不是整數(shù)周期，情況就不同了。把頻率改為11Hz，而采樣頻率仍為64Hz和窗長(zhǎng)度仍為32點(diǎn)，用同樣的程序計(jì)算此余弦信號(hào)的頻譜，可以得到，頻譜圖上有兩個(gè)較大的峰值，而其他的點(diǎn)上幅度也不再為零。，因此可以判斷，信號(hào)的頻譜峰值確實(shí)位于兩者之間。本來(lái)是單一的11Hz頻率的能量會(huì)分布到許多DFT頻率上的現(xiàn)象稱為頻率泄漏，來(lái)源于截?cái)嘈?yīng)。</p><

26、p>　　2. 考察DFT的長(zhǎng)度對(duì)雙頻率信號(hào)頻譜分析的影響。設(shè)待分析的信號(hào)為</p><p>　　令兩個(gè)長(zhǎng)度為16的正余弦序列的數(shù)字頻率為及。取N為四個(gè)不同值16，32，64，128。畫(huà)出四個(gè)DFT幅頻圖，分析DFT長(zhǎng)度對(duì)頻譜分辨率的影響。</p><p>　　n=0:15; N1=16;N2=32;N3=64;N4=128;f1=0.22;f2=0.34;</p>&l

27、t;p>　　x=0.5*sin(2*pi*f1*n)+sin(2*pi*f2*n); %待分析信號(hào)</p><p><b>　　%譜分析</b></p><p>　　X1=abs(fft(x,N1));</p><p>　　X2=abs(fft(x,N2));</p><p>　　X3=abs(fft(x

28、,N3));</p><p>　　X4=abs(fft(x,N4));</p><p>　　k1=0:N1-1;k2=0:N2-1;k3=0:N3-1;k4=0:N4-1;</p><p>　　subplot(4,2,1),stem(n,x); %繪制序列的波形</p><p>　　xlabel('n

29、');ylabel('x1(n)');title('余弦序列');</p><p>　　subplot(4,2,2);stem(k1,X1); %繪制N1的幅頻特性</p><p>　　xlabel('n');ylabel('X(n)'); title('N=16的幅頻特性曲線圖');

30、</p><p>　　subplot(4,2,3),stem(n,x); %繪制序列的波形</p><p>　　xlabel('n');ylabel('x1(n)');title('余弦序列');</p><p>　　subplot(4,2,4);stem(k2,X2);

31、 %繪制N2的幅頻特性</p><p>　　xlabel('n');ylabel('X(n)'); title('N=32的幅頻特性曲線圖');</p><p>　　subplot(4,2,5),stem(n,x); %繪制序列的波形</p><p>　　xlabel('n

32、');ylabel('x1(n)');title('余弦序列');</p><p>　　subplot(4,2,6);stem(k3,X3); %繪制N3的幅頻特性</p><p>　　xlabel('n');ylabel('X(n)'); title('N=64的幅頻特性曲線圖');

33、</p><p>　　subplot(4,2,7),stem(n,x); %繪制序列的波形</p><p>　　xlabel('n');ylabel('x1(n)');title('余弦序列');</p><p>　　subplot(4,2,8);stem(k4,X4);

34、 %繪制N4的幅頻特性</p><p>　　xlabel('n');ylabel('X(n)'); title('N=128的幅頻特性曲線圖')</p><p>　　圖3. N=16的DFT幅頻圖 </p><p>　　圖4. N=32的DFT幅頻圖</p><p>　　圖5. N=64的D

35、FT幅頻圖</p><p>　　圖6. N=128的DFT幅頻圖</p><p><b>　　分析：</b></p><p>　　采樣速率Fs>2fc，譜分辨率F=Fs/N，如果保持采樣點(diǎn)數(shù)N不變，要提高頻譜分辨率（減小F），就必須降低采樣頻率，采樣頻率的降低會(huì)引起譜分析范圍變窄和譜頻率混疊失真。如維持Fs不變，為提高頻譜分辨率可以增加采

36、樣點(diǎn)數(shù)N。更長(zhǎng)的時(shí)域信號(hào)能夠提供更高的頻域分辨率，因?yàn)橐粋€(gè)N點(diǎn)的時(shí)域信號(hào)能被分解為N/2+1個(gè)余弦信號(hào)和N/2+1個(gè)正弦信號(hào)，N增大則(N/2+1)也增大，頻域間隔(1/2的時(shí)域采樣頻率)/(N/2+1)減小，所以頻域分辨率提高了。所以利用DFT計(jì)算頻譜時(shí)增加取樣點(diǎn)的長(zhǎng)度范圍可以提高分辨率。</p><p>　　3. 在上題中若把兩個(gè)正弦波的頻率取得較近，令，，試問(wèn)怎樣選擇FFT參數(shù)才能在頻譜分析中分辨出這兩個(gè)分

37、量？</p><p>　　N=input('輸入譜分析的長(zhǎng)度');</p><p><b>　　n=1:N-1;</b></p><p>　　f1=0.22;f2=0.25;</p><p>　　x=0.5*sin(2*pi*f1*n)+sin(2*pi*f2*n); %待分析信號(hào)</p&

38、gt;<p>　　subplot(1,2,1),stem(n,x); %繪制序列的波形</p><p>　　xlabel('n');ylabel('x1(n)');</p><p>　　title('余弦序列');</p><p>　　X=abs(fft(x,N));

39、 %求余弦序列的N點(diǎn)FFT</p><p>　　subplot(1,2,2);</p><p><b>　　k=0:N-1;</b></p><p>　　stem(k,X); %繪制序列的幅頻特性曲線</p><p&

40、gt;　　xlabel('k');ylabel('X(k)'); </p><p>　　string=[num2str(N),'點(diǎn)FFT幅頻曲線'];</p><p>　　title(string);</p><p>　　圖7. N=16的FFT幅頻曲線</p><p>　　圖8. N=3

41、2的FFT幅頻曲線</p><p>　　圖9. N=64的FFT幅頻曲線</p><p>　　圖10. N=128的FFT幅頻曲線</p><p><b>　　分析：</b></p><p>　　從題可以看出當(dāng)Np=1/|f1-f2|=1/(0.25-0.22)約為33，而N=32最接近33，因?yàn)?2點(diǎn)的分辨率最好。只有

42、當(dāng)Np=2*pi/|f1-f2|時(shí)才能很好的恢復(fù)原信號(hào)的頻譜。因而，對(duì)于雙頻信號(hào)而言，要增大其頻譜的分辨率不僅與其窗函數(shù)的寬度N有關(guān)還與兩頻率差有關(guān)。因此，隨著N的增大，其頻譜的分辨率增大，但是當(dāng)N增加到一定時(shí)，其頻譜分辨率反而下降。因而不能通過(guò)無(wú)線的增大N來(lái)擴(kuò)大頻譜的分辨率。其第2題中Np=1/(0.34-0.22)約等于8，因而四個(gè)中16點(diǎn)的頻譜的分辨率最好。</p><p><b>　　4.擴(kuò)展練

43、習(xí)：</b></p><p>　　DFT長(zhǎng)度對(duì)單頻率信號(hào)頻譜分析的影響。</p><p>　　對(duì)連續(xù)的單一頻率周期信號(hào)按采樣頻率采樣，截取長(zhǎng)度N，分別選N=20和N=16，觀察其DFT結(jié)果的幅度譜。此時(shí)離散序列，即k=8。</p><p>　　N=input('輸入譜分析的長(zhǎng)度');</p><p><b&g

44、t;　　n=1:N-1;</b></p><p><b>　　f1=1/8;</b></p><p>　　x=sin(2*pi*f1*n); %待分析信號(hào)</p><p>　　subplot(1,2,1),stem(n,x); %繪制序列的波形<

45、;/p><p>　　xlabel('n');ylabel('x1(n)');</p><p>　　title('余弦序列');</p><p>　　X=abs(fft(x,N)); %求余弦序列的N點(diǎn)DFT</p><p>　　subplot(1,2

46、,2);</p><p><b>　　k=0:N-1;</b></p><p>　　stem(k,X); %繪制序列的幅頻特性曲線</p><p>　　xlabel('k');ylabel('X(k)'); </p><p&

47、gt;　　string=[num2str(N),'點(diǎn)DFT幅頻曲線'];</p><p>　　title(string);</p><p>　　圖11. N=20的單頻DFT頻譜曲線</p><p>　　圖12. N=16的單頻DFT頻譜曲線</p><p><b>　　分析：</b></p>

48、<p>　　計(jì)算結(jié)果示于圖11圖12，圖11分別是N=20時(shí)的截取信號(hào)和DFT結(jié)果，由于截取了兩個(gè)半周期，頻譜出現(xiàn)泄漏；圖12分別是N=16時(shí)的截取信號(hào)和DFT結(jié)果，由于截取了兩個(gè)整周期得到單一譜線的頻譜。上述頻譜的誤差主要是由于時(shí)域中對(duì)信號(hào)的非整周期截?cái)喈a(chǎn)生的頻譜泄漏 。</p><p><b>　　結(jié)束語(yǔ)：</b></p><p>　　通過(guò)此次課程設(shè)

49、計(jì)，使我更加扎實(shí)的掌握了有關(guān)數(shù)字信號(hào)處理以及MATLAB方面的知識(shí)，在設(shè)計(jì)過(guò)程中雖然遇到了一些問(wèn)題，但經(jīng)過(guò)一次又一次的思考，一遍又一遍的檢查以及與同學(xué)老師的交流終于找出了原因所在，也暴露出了前期我在這方面的知識(shí)欠缺和經(jīng)驗(yàn)不足。實(shí)踐出真知，通過(guò)親自動(dòng)手制作，使我們掌握的知識(shí)不再是紙上談兵。</p><p>　　在課程設(shè)計(jì)過(guò)程中，我們不斷發(fā)現(xiàn)錯(cuò)誤，不斷改正，不斷領(lǐng)悟，不斷獲取。最終的檢測(cè)調(diào)試環(huán)節(jié)，本身就是在踐行“過(guò)而

50、能改，善莫大焉”的知行觀。這次課程設(shè)計(jì)終于順利完成了，在設(shè)計(jì)中遇到了很多問(wèn)題，最后在老師的指導(dǎo)下，終于游逆而解。在今后社會(huì)的發(fā)展和學(xué)習(xí)實(shí)踐過(guò)程中，一定要不懈努力，不能遇到問(wèn)題就想到要退縮，一定要不厭其煩的發(fā)現(xiàn)問(wèn)題所在，然后一一進(jìn)行解決，只有這樣，才能成功的做成想做的事，才能在今后的道路上劈荊斬棘，而不是知難而退，那樣永遠(yuǎn)不可能收獲成功，收獲喜悅，也永遠(yuǎn)不可能得到社會(huì)及他人對(duì)你的認(rèn)可。</p><p><b&

51、gt;　　參考文獻(xiàn)：</b></p><p>　　[1]  張威 . MATLAB基礎(chǔ)與編程入門 . 西安電子科技大學(xué)出版社，20008.1[2]  程佩青  數(shù)字信號(hào)處理教程  清華大學(xué)出版社，2010.5[3]  曹志剛、錢亞生  現(xiàn)代

52、通信原理  清華大學(xué)出版社  2010.6[4]  王福昌  通信原理學(xué)習(xí)輔導(dǎo)  華中科技大學(xué)出版社  2008.8[5]  高成  matlab圖像處理與應(yīng)用  國(guó)防工業(yè)出版社  2007.4 </p><p>