dsp設(shè)計(jì)基于MATLAB的數(shù)字譜分析實(shí)踐及誤差討論

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

2、及程度，找到減少誤差的方法。理解數(shù)字譜分析中物理分辨力和計(jì)算分辨力的概念，討論兩者間的關(guān)系。 二 項(xiàng)目?jī)?nèi)容1. 舉例說明時(shí)窗寬度、類型與物理分辨率的關(guān)系。得出相應(yīng)結(jié)論，給出提高物理分辨率的方法。2. 舉例說明時(shí)窗類型與頻譜泄露量的關(guān)系，說明頻譜泄露的危害。給出降低頻譜泄露量的方法。3. 舉例說明物理分辨率與計(jì)算分辨率的關(guān)系，如何避免頻率分量的漏判與誤判。4. 進(jìn)行課堂演示及討論(PPT)。5. 分組討論，獨(dú)立撰寫項(xiàng)目報(bào)告，提交Word電子文檔。三 項(xiàng)目報(bào)告正文1. 數(shù)字譜分析的意義 信號(hào)特征提取，語音、圖象處理，通信信號(hào)處理等。應(yīng)用范圍：家用電器，儀器儀表，醫(yī)療儀器，有線通信，無線通信，生物生

3、理信息處理，地球物理信息探測(cè)，雷達(dá)，電子對(duì)抗，航空航天，宇宙探索，動(dòng)力控制系統(tǒng)幾乎所有電子設(shè)備和信息語音壓縮（MP3、無線通信）用到FFT，頻譜分析儀用到FFT，4G的關(guān)鍵技術(shù)OFDM用到FFT2. DTFT、DFT和FFT算法概述用Equation編輯器編輯公式，下同。3. 時(shí)窗效應(yīng)的仿真分析3.1 矩形窗和海明窗各自與物理分辨率的關(guān)系3.1.1 Rectangular Window 當(dāng)時(shí)窗寬度L=64時(shí)，設(shè)定f1=1kHz，f2=1.05 kHz，f3=1.10 kHz；fs=5 kHz。代碼為：clear;close all;L=64;n=0:L-1;f1=1000;f2=1050;f3

4、=1100;fs=5000;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/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);仿真結(jié)果如下： 當(dāng)時(shí)窗寬度L=200時(shí)，設(shè)定f1=1kHz，f2=1.05 kHz，f3=1.10 kHz；fs=5 kHz。代碼：（略）仿真結(jié)果如下：仿真結(jié)果分析：由物理分辨率計(jì)算公式知： （c=1），則當(dāng)L=64

5、，fs=5kHz時(shí)=78.125 Hz， 又因?yàn)?個(gè)信號(hào)的譜峰間隔為50Hz，則.物理分辨率大于信號(hào)頻率間隔，故不能分辨出3個(gè)信號(hào)的頻譜；當(dāng)L=200，fs=5kHz時(shí)，=25Hz，則.物理分辨率大于信號(hào)頻率間隔，故不能分辨出3個(gè)信號(hào)的頻譜；當(dāng)L=300，fs=5kHz時(shí)，=33.3Hz，則 物理分辨率小于信號(hào)頻率間隔，故能分辨出3個(gè)信號(hào)的頻譜結(jié)論：取相同的時(shí)窗寬度L，則矩形窗的物理分辨率比海明窗高。由物理分辨率計(jì)算公式可知，在采樣頻率fs一定時(shí)，增加時(shí)窗寬度L，兩種類型的時(shí)窗的物理分辨率均提高。提高物理分辨率的方法是：在采樣率不變的前提下增加時(shí)窗寬度L或者使用矩形窗。3.2矩形窗和海明窗各自

6、與頻譜泄露量的關(guān)系3.2.1 Rectangular leakage 設(shè)時(shí)窗寬度L=100時(shí)，設(shè)定f1=0.5kHz，f2=1 kHz，f3=1.5 kHz；fs=5 kHz。代碼如下：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,200

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

8、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;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);仿真結(jié)果如下：仿真結(jié)果分析：由圖可知，信號(hào)以海明窗截取，頻譜泄露大為減小。結(jié)論：在滿足物理分辨率的前提下，用海明窗截取的信號(hào)的頻譜泄露遠(yuǎn)比用矩形窗截取的少。頻譜泄露不僅會(huì)降低信號(hào)功率

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

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

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

12、理分辨率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_w =0: 2*pi/N

13、: 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)subplot(3,1,2)axis(0 fs/2 0 12)hold onbar (0:N2-1)/N2*fs

14、,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*) plot(fs*vct_w/pi/2,abs(X),b:) legend(DFT, L=10 N=256,DTFT,0)3.計(jì)算分辨率DFT/DTFT仿真clfL2=100; L=800;N1=16

15、;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),0) plot(0:HN1-1)/N1*fs,abs(X1),b*) plot(fs*vct_w/pi/2,abs(X),b:)legend(DFT, L=100 N=16,DTFT,0)subplot(3,1,2)axis(0 fs/2 0 60)hold onNN HN2=size(abs(X2);bar (0:HN2-1)/N2*fs,abs(X2),0) plot(0:HN2-1)/N2*fs,abs(X2),b*) plot(fs*vct_w/pi/2,abs(X),b:)l