電子信息工程數(shù)字譜分析實踐及誤差討論

1、汕頭大學工學院三級項目報告課程名稱：數(shù) 字 信 號 處 理課程設計題目：數(shù)字譜分析實踐及誤差討論指導教師：系別：電子工程系專業(yè)：電子信息工程學號：姓名：合作者完成時間：2010年10月6日至10月26日成績：評閱人：一 項目意義與目標數(shù)字譜分析是DSP應用系統(tǒng)最基礎算法。如語音壓縮（MP3、無線通信）用到FFT，頻譜分析儀用到FFT，4G的關鍵技術OFDM用到FFT在實際DSP系統(tǒng)中能夠正確運用數(shù)字譜分析算法。通過親自編程、繪頻譜圖，掌握數(shù)字譜分析有關算法DTFT、DFT和FFT，并能夠恰當舉例說明數(shù)字譜分析算法產生誤差原因及程度，找到減少誤差的方法。理解數(shù)字譜分析中物理分辨力和計算分辨力的概

2、念，討論兩者間的關系。二 項目內容1. 舉例說明時窗寬度、類型與物理分辨率的關系。得出相應結論，給出提高物理分辨率的方法。2. 舉例說明時窗類型與頻譜泄露量的關系，說明頻譜泄露的危害。給出降低頻譜泄露量的方法。3. 舉例說明物理分辨率與計算分辨率的關系，如何避免頻率分量的漏判與誤判。4. 進行課堂演示及討論(PPT)。5. 分組討論，獨立撰寫項目報告，提交Word電子文檔。三 項目報告正文1. 數(shù)字譜分析的意義信號特征提取，語音、圖象處理，通信信號處理等。應用范圍：家用電器，儀器儀表，醫(yī)療儀器，有線通信，無線通信，生物生理信息處理，地球物理信息探測，雷達，電子對抗，航空航天，宇宙探索，動力控制

3、系統(tǒng)幾乎所有電子設備和信息語音壓縮（MP3、無線通信）用到FFT，頻譜分析儀用到FFT，4G的關鍵技術OFDM用到FFT2. DTFT、DFT和FFT算法概述用Equation編輯器編輯公式，下同。3. 時窗效應的仿真分析3.1 矩形窗和海明窗各自與物理分辨率的關系3.1.1Rectangular Window當時窗寬度L=64時，設定f1=1kHz，f2=1.05kHz，f3=1.10kHz；fs=5kHz。代碼為：clear;close all;L=64;n=0:L-1;f1=1000;f2=1050;f3=1100;fs=5000;w1=2*pi*f1/fs;w2=2*pi*f2/fs;

4、w3=2*pi*f3/fs;xL=sin(w1*n)+sin(w2*n)+sin(w3*n);w=0:2*pi/300:2*pi;yL=dtft(xL,w);plot(w/2/pi*fs,abs(yL),'k')axis(500,2500,0,L/1.5)xlabel('f/Hz')gtext('L=64');仿真結果如下：當時窗寬度L=200時，設定f1=1kHz，f2=1.05kHz，f3=1.10kHz；fs=5kHz。代碼：（略）仿真結果如下：仿真結果分析：由物理分辨率計算公式知：（c=1），則當L=64，fs=5kHz時=78.125

5、Hz， 又因為3個信號的譜峰間隔為50Hz，則>.物理分辨率大于信號頻率間隔，故不能分辨出3個信號的頻譜；當L=200，fs=5kHz時，=25Hz，則<物理分辨率小于信號頻率間隔，故能分辨出3個信號的頻譜3.1.2 Hamming window當時窗寬度L=100時，設定f1=1kHz，f2=1.05kHz，f3=1.10kHz；fs=5kHz。代碼為：clear;close all;L=100;n=0:L-1;fs=5000;f1=1000;f2=1050;f3=1100;w1=2*pi*f1/fs;w2=2*pi*f2/fs;w3=2*pi*f3/fs;xL=sin(w1*n

6、)+sin(w2*n)+sin(w3*n);w=0:2*pi/500:2*pi;wh=0.54-0.46*cos(2*pi*n/(L-1);xL_H=xL.*wh;yL=dtft(xL_H,w);plot(w/2/pi*fs,abs(yL),'k')axis(500,1500,0,L/3)xlabel('f/Hz')gtext('L=100');仿真結果如下：當時窗寬度L=300時，設定f1=1kHz，f2=1.05kHz，f3=1.10kHz；fs=5kHz。代碼：clear;close all;L=300;n=0:L-1;fs=5000;f1

7、=1000;f2=1050;f3=1100;w1=2*pi*f1/fs;w2=2*pi*f2/fs;w3=2*pi*f3/fs;xL=sin(w1*n)+sin(w2*n)+sin(w3*n);w=0:2*pi/1000:2*pi;wh=0.54-0.46*cos(2*pi*n/(L-1);xL_H=xL.*wh;yL=dtft(xL_H,w);plot(w/2/pi*fs,abs(yL),'k')axis(500,1500,0,L/3)xlabel('f/Hz')gtext('L=300');仿真結果如下：仿真結果分析：由物理分辨率計算公式知：

8、 ，（c=2）則當L=100，fs=5kHz時=100 Hz， 又因為3個信號的譜峰間隔為50Hz，則>.物理分辨率大于信號頻率間隔，故不能分辨出3個信號的頻譜；當L=300，fs=5kHz時，=33.3Hz，則<物理分辨率小于信號頻率間隔，故能分辨出3個信號的頻譜結論：取相同的時窗寬度L，則矩形窗的物理分辨率比海明窗高。由物理分辨率計算公式可知，在采樣頻率fs一定時，增加時窗寬度L，兩種類型的時窗的物理分辨率均提高。提高物理分辨率的方法是：在采樣率不變的前提下增加時窗寬度L或者使用矩形窗。3.2矩形窗和海明窗各自與頻譜泄露量的關系 Rectangular leakage設時窗寬度

9、L=100時，設定f1=0.5kHz，f2=1kHz，f3=1.5kHz；fs=5kHz。代碼如下：clear;close all;L=100;n=0:L-1;f1=500;f2=1000;f3=1500;fs=5000;w1=2*pi*f1/fs;w2=2*pi*f2/fs;w3=2*pi*f3/fs;xL=sin(w1*n)+0.1*sin(w2*n)+sin(w3*n);w=0:2*pi/500:2*pi;yL=dtft(xL,w);plot(w/2/pi*fs,abs(yL),'k')axis(0,2000,0,L/1.5)xlabel('f/Hz')g

10、text('L=100');仿真結果如下：仿真結果分析：由圖可知，信號以矩形窗截取，有較大的頻譜泄露。頻率為0.5kHz及1.5kHz的信號頻譜，最靠近譜峰的副瓣幅值比頻率為1kHz的信號的譜峰值還要高，導致難以根據譜峰值區(qū)分信號，從而引致信號的缺失。Hamming leakage設時窗寬度L=200時，設定f1=0.5kHz，f2=1kHz，f3=1.5kHz；fs=5kHz。代碼為：clear;close all;L=200;n=0:L-1;fs=5000;f1=500;f2=1000;f3=1500;w1=2*pi*f1/fs;w2=2*pi*f2/fs;w3=2*pi*

11、f3/fs;xL=sin(w1*n)+0.1*sin(w2*n)+sin(w3*n);w=0:2*pi/500:2*pi;wh=0.54-0.46*cos(2*pi*n/(L-1);xL_H=xL.*wh;yL=dtft(xL_H,w);plot(w/2/pi*fs,abs(yL),'k')axis(0,2000,0,L/3)xlabel('f/Hz')gtext('L=200');仿真結果如下：仿真結果分析：由圖可知，信號以海明窗截取，頻譜泄露大為減小。結論：在滿足物理分辨率的前提下，用海明窗截取的信號的頻譜泄露遠比用矩形窗截取的少。頻譜泄露不

12、僅會降低信號功率而且還會掩蓋一些微弱信號，導致信號缺失。我們可以通過增加采樣的長度或者加海明窗來降低頻譜泄露量。4. DFT/DTFT計算分辨率與物理分辨率關系的仿真分析4.1 物理分辨率仿真 代碼見附錄。仿真結果如下： 圖14.2 計算分辨率仿真：（代碼見附錄）仿真結果如下： 圖2仿真結果分析：由上圖1可知，當物理分辨率低于信號頻率間隔時，無論計算分辨率再大，也無法從頻譜圖中把信號的頻率分量分辨出來；相反，由圖2知，當物理分辨率高于信號頻率間隔時，隨著計算分辨率的加大頻譜圖像取點更豐富，N=16時，3個頻率分量中只有頻率f=2.5KHz有對應的頻譜點，很可能造成頻率分量的漏判；N=64時，基

13、本可以判斷有3個頻率分量，除了頻率f=2.5KHz能正確對應頻譜，其余兩個分量的頻譜點與其對應頻率都有偏差，這時候就很可能造成誤判；N=256時，清楚顯示3個頻率分量及其對應頻率。綜上所述，在滿足物理分辨率的前提下，增加計算分辨率，可以有效避免頻率分量的漏判與誤判。5. 數(shù)字譜分析中存在的誤差及減小誤差的方法Aliasing, Frequency resolution, Frequency leakage四 總結 通過撰寫這次三級項目的報告，我深深體會到我們CDIO理念中的Conceive和Implement的意義。隨著這次項目的開展，我不斷的自我增值，獨立撰寫報告使我的自學研究能力有所提高，

14、和同學們的討論在增強我的溝通及語言表達能力的同時也加強了我的團隊合作精神。開始時，為了更深入了解數(shù)字譜分析，我翻閱英文教材以及研讀老師的課件。為了實現(xiàn)數(shù)字譜分析的MATLAB仿真，我從圖書館借來MATLAB教程自學MATLAB軟件，并不斷向其他精通MATLAB的同學請教。從開始構思如何分析數(shù)字譜到最終通過MATLAB軟件仿真來實現(xiàn)構想，期間我不但加深了對數(shù)字譜分析的理解，還學會了運用MATLAB來對某些理論進行仿真驗證，真是可謂獲益良多。五 附錄部分代碼1.DTFT函數(shù)代碼：function X=dtft(x,w)%x為輸入離散序列X = x*exp(-j*1:length(x)'*w

15、); %w為頻率數(shù)組2.物理分辨率DFT/DTFT仿真clear allclfj=sqrt(-1);f1=2; f2=2.5; f3=3;fs=10;w1=2*pi*f1/fs; w2=2*pi*f2/fs; w3=2*pi*f3/fs;L1=10; L=800;N1=16; N2=64; N3=256; N=800;n=0:L1-1;x=sin(w1*n)+sin(w2*n)+sin(w3*n);vct_w1=0: 2*pi/N1: 2*pi-2*pi/N1;vct_w2=0: 2*pi/N2: 2*pi-2*pi/N2;vct_w3=0: 2*pi/N3: 2*pi-2*pi/N3;vct

16、_w =0: 2*pi/N: 2*pi-2*pi/N;X1=dtft(x,vct_w1); X2=dtft(x,vct_w2);X3=dtft(x,vct_w3);X =dtft(x,vct_w);subplot(3,1,1)axis(0 fs/2 0 12)hold onbar (0:N1-1)/N1*fs,abs(X1),0) plot(0:N1-1)/N1*fs,abs(X1),'b*') plot(fs*vct_w/pi/2,abs(X),'b:')legend('DFT, L=10 N=16','DTFT',0)subp

17、lot(3,1,2)axis(0 fs/2 0 12)hold onbar (0:N2-1)/N2*fs,abs(X2),0) plot(0:N2-1)/N2*fs,abs(X2),'b*') plot(fs*vct_w/pi/2,abs(X),'b:')legend('DFT, L=10 N=64','DTFT',0)subplot(3,1,3)axis(0 fs/2 0 12)hold onbar (0:N3-1)/N3*fs,abs(X3),0) plot(0:N3-1)/N3*fs,abs(X3),'b*'

18、) plot(fs*vct_w/pi/2,abs(X),'b:')legend('DFT, L=10 N=256','DTFT',0)3.計算分辨率DFT/DTFT仿真clfL2=100; L=800;N1=16;N2=64;N3=256;n=0:L2-1;x=cos(w1*n)+cos(w2*n)+cos(w3*n);vct_w1=0: 2*pi/N1: pi;vct_w2=0: 2*pi/N2: pi;vct_w3=0: 2*pi/N3: pi;vct_w =0: 2*pi/N: pi;X1=dtft(x,vct_w1);X2=dtft(x,vct_w2);X3=dtft(x,vct_w3);X =dtft(x,vct_w);subplot(3,1,1)axis(0 fs/2 0 60)hold onNN HN1=size(abs(X1);bar (0:HN1-1)/N1*fs,abs(X1),