利用經(jīng)典譜估計(jì)法估計(jì)信號(hào)的功率譜(隨機(jī)信號(hào))

1、精品文檔 1 歡。迎下載 隨機(jī)信號(hào) 利用經(jīng)典譜估計(jì)法估計(jì)信號(hào)的功率譜精品文檔 1 歡。迎下載 作業(yè)綜述： 給出一段信號(hào) “asd.wav asd.wav ”, , 利用經(jīng)典譜估計(jì)法的原理， 通過(guò)不同的譜估計(jì)方法， 求出信 號(hào)的功率譜密度函數(shù)。采用 MATLABMATLAB 語(yǔ)言，利用 MATLAMATLA 萌言強(qiáng)大的數(shù)據(jù)處理和數(shù)據(jù)可視化能 力，通過(guò) GUIGUI 的對(duì)話(huà)框模板，使操作更為簡(jiǎn)便！在一個(gè) GUIGUI 界面中，同時(shí)呈現(xiàn)出不同方法產(chǎn) 生出的功率譜。 這里給出了幾種不同的方法： BTBT 法,周期圖法，平均法以及 WelchWelch 法。把幾種不同方法所 得到的功率譜都呈現(xiàn)在一個(gè)界面

2、中，便于對(duì)幾種不同方法得到的功率譜作對(duì)比。 一 題 目要求 給出一段信號(hào)及采樣率，利用經(jīng)典譜估計(jì)法估計(jì)出信號(hào)的功率譜。 二基本原理及方法 經(jīng)典譜估計(jì)的方法，實(shí)質(zhì)上依賴(lài)于傳統(tǒng)的傅里葉變換法。它是將數(shù)據(jù)工作區(qū)外的未知數(shù) 據(jù)假設(shè)為零，相當(dāng)于數(shù)據(jù)加窗，主要方法有 BTBT 法，周期圖法，平均法以及 WelchWelch 法。 1. BT 法 (Blackman-Tukey) 理論基礎(chǔ)： (1) (1) 隨機(jī)序列的維納 - -辛欽定理 由于隨機(jī)序列X(n)的自相關(guān)函數(shù) Rx(m)=EX(n)X(n+m)定義在離散點(diǎn) m上，設(shè)取 樣間隔為 ，則可將隨機(jī)序列的自相關(guān)函數(shù)用連續(xù)時(shí)間函數(shù)表示為 等式兩邊取傅里葉

3、變換，則隨機(jī)序列的功率譜密度 (2) (2) 譜估計(jì) BTBT 法是先估計(jì)自相關(guān)函數(shù) Rx(m)(|m|=0,1,2,N-1), ,然后再經(jīng)過(guò)離散傅里葉變換求 的功率譜密度的估值 。即 其中 可有式得到。 2. 周期圖法 理論基礎(chǔ)： 周期圖法是根據(jù)各態(tài)歷經(jīng)隨機(jī)過(guò)程功率譜的定義來(lái)進(jìn)行譜估計(jì)的。在前面我們已知， 各態(tài)歷經(jīng)的連續(xù)隨機(jī)過(guò)程的功率譜密度滿(mǎn)足精品文檔 2歡迎下載 葉變換為: 因此有隨機(jī)信號(hào)的觀測(cè)數(shù)據(jù) 的功率譜估計(jì)值（稱(chēng)“周期圖”）如下: 由于上式中的離散傅里葉變換可以用快速傅里葉變換計(jì)算， 因此就可以估計(jì)出功率 譜。 3. 平均法: 理論基礎(chǔ): 平均法可視為周期圖法的改進(jìn)。 目的，有一個(gè)定

4、理：如果 ， 周期圖經(jīng)過(guò)平均后會(huì)使它的方差減少， 達(dá)到一致估計(jì)的 , 是不相關(guān)的隨機(jī)變量，且都有個(gè)均值 及其方差 ， 則可以證明它們的算術(shù)平均的均值為 ，方差為 。 由定理可見(jiàn)：具有 個(gè)獨(dú)立冋分布隨機(jī)變量平均的方差，是單個(gè)隨機(jī)變量方差的 ，時(shí)，方差 ，可以達(dá)到一致估計(jì)的目的。因此，將 個(gè)獨(dú)立的估計(jì)量經(jīng)過(guò)算術(shù) 平均后得到的估計(jì)量的方差也是原估計(jì)量方差的 平均圖法即是將數(shù)據(jù) ， 個(gè)數(shù)據(jù)樣本（平均法適用于數(shù)據(jù)量大的場(chǎng)合） 別計(jì)算每一段的周期圖 然后將這 1010 個(gè)周期圖加以平均得譜估計(jì)值: 1 10 G100 w G100, i w 10 i# 由于這 1010 小段的周期圖取決于同一個(gè)過(guò)程，因而

5、其均值相同。若這 1010 個(gè)小段的周期圖 是統(tǒng)計(jì)獨(dú)立的，則這 1010 個(gè)小段平均之后的方差卻是單段方差的式中 是連續(xù)隨機(jī)過(guò)程第i個(gè)樣本的截取函數(shù) 的頻譜。對(duì)應(yīng)在隨機(jī)序列中 由于隨機(jī)序列中觀測(cè)數(shù)據(jù) 僅在 的點(diǎn)上存在，貝 U U 的 N N 點(diǎn)離散傅里 分段求周期圖法后再平均。例如，給定 N=1000N=1000 ，則可以將它分成 1010 個(gè)長(zhǎng)度為 100100 的小段，分 Gioo, i 1 100 1001 -4 |n 壬 00(1書(shū) 2 X(ej ) ,1 =1,2,、,10 精品文檔 3歡迎下載 D = DD = D 即：平均法將 ， 的 N N 個(gè)數(shù)據(jù)分成 L L 段(N=MLN=

6、ML,若各數(shù)據(jù)段相互獨(dú)立，則 平方后估計(jì)量的方差是原來(lái)不分段估計(jì)量方差的 。所以當(dāng) 時(shí)，估計(jì)量的方差趨于 0 0， 達(dá)到一致估計(jì)的目的。但是，隨著分段數(shù) L L 的增加，M M 點(diǎn)數(shù)減少，分辨率減少，使估計(jì)變成 有偏估計(jì)。相反，若 L L 減少，M M 增加，雖偏差減少，但方差增大。所以，在應(yīng)用中，必須兼 顧分辨率和方差的要求來(lái)適當(dāng)選擇 M M 和 L L 的值。 4. Welch 法: 理論基礎(chǔ)： WelchWelch 法又稱(chēng)修正周期法，其步驟為： 先將 N N 個(gè)數(shù)據(jù)分成 L L 段，每數(shù)據(jù)段 M M 個(gè)數(shù)據(jù)，N=ML N=ML 選擇適當(dāng)?shù)拇昂瘮?shù) w w(n n),并用該 w(nw(n)依

7、次對(duì)每段數(shù)據(jù)做相應(yīng)的加權(quán)，然后確定 每段的周期圖 U為歸一化因子 1 M丄2 U w2( n), M n 對(duì)每段周期圖進(jìn)行平均得到功率譜估計(jì)： 1 L GM w GM,1 w L i# 當(dāng)數(shù)據(jù)量一定時(shí)，若分段數(shù) L L 增加，M M 點(diǎn)數(shù)減少，則分辨率下降；若 L L 減少，雖 M M 增加，但 方差增大。解決這一矛盾的方法是，讓數(shù)據(jù)段間適當(dāng)“重疊” 三. 算法設(shè)計(jì)與實(shí)現(xiàn) 1 程序流程圖 (1) BT法的流程圖GM ,I w 1 MU Ml J 送 n=M(l d) 2 Xnw(n,I=1,2,、,L 精品文檔 4歡迎下載 (2) 周期圖法的流程圖 開(kāi)始 精品文檔 5歡迎下載 (3) 平均法的

8、流程圖 (4) Welch法的流程圖精品文檔 6歡迎下載 先將 N 個(gè)數(shù)據(jù)分成 L 段 2.主要模塊的設(shè)計(jì)： (1)產(chǎn)生原始信號(hào) Fs=600; n fft=1024; n=0:1/Fs:1; y=cos(2*pi*30* n)+3*cos(2*pi*100* n)+ra ndn (size( n); axes(ha ndles.axes1); Plot (n ,y); title( 原始信號(hào)波形 精品文檔 7歡迎下載 原始信號(hào)波形如下： (2)BT 法 global y Fs nfft n; cxn=xcorr(y,u nbiased); CXk=fft(cx n,n fft); Pxx=a

9、bs(CXk); in dex=0:rou nd( nfft/2-1); k=in dex*Fs/ nfft; plot_Pxx=10*log(Pxx(i ndex+1)/(log(10); axes(ha ndles.axes2); plot(k,plot_Pxx); xlabel( 頻率(hz); ylabel(功率譜密度(Db/Hz); title(BT 法); 周期圖法 global y Fs nfft n; XF=fft(y,nfft); Pxx=abs(XF).A2/le ngth( n); in dex=0:rou nd( nfft/2-1); f=in dex*Fs/nfft;

10、精品文檔 8歡迎下載 axes(ha ndles.axes4); plot(f,10*log(Pxx(i ndex+1); xlabel( 頻率(hz); ylabel(功率譜密度(Db/Hz); title( 周期圖法); 周期圖法所得功率譜波形如下: 周期團(tuán);去 (4)平均法 global y Fs nfft n; wi ndow=hammi ng(nfft); no verlap=0; p=0.9; Pxx,Pxxc=psd(y, nfft,Fs,wi ndow, no verlap,p); in dex=0:rou nd( nfft/2-1); k=in dex*Fs/ nfft; Q

11、xx=10*log10(Pxx(i ndex+1); axes(ha ndles.axes6); plot(k,Qxx ); xlabel( 頻率(hz); ylabel(功率譜密度(Db/Hz); title( 平均法); 平均法所得功率譜波形如下： 平均法 頻率(hz) 100 200 30Q 頻率(hz) 400 500 2020O OO O- -2 2精品文檔 9歡迎下載 (5) Welch 法 global y Fs nfft n; win dow1=ham min g(100); no verlap=20; ran ge=half; Pxx1,f=pwelch(y,wi ndow1

12、, no verlap, nfft,Fs,ra nge); Qxx 1=10*log10(Pxx1); axes(ha ndles.axes5); plot(f,Qxx1); xlabel( 頻率(hz); ylabel(功率譜密度(Db/Hz); title(Welch 法); Welch 法所得功率譜波形如下： 四. 軟件使用說(shuō)明 綜述 從給出一段信號(hào) y=cos(2*pi*30* n)+3*cos(2*pi*100* n),y=cos(2*pi*30* n)+3*cos(2*pi*100* n), 利用經(jīng)典譜估計(jì)法的原理，通 過(guò)不同的譜估計(jì)方法，求出信號(hào)的功率譜密度函數(shù)，并利用 GUIG

13、UI 界面呈現(xiàn)出不同譜估計(jì)方法 所得的結(jié)果。 BT法 BTBT 法是先估計(jì)自相關(guān)函數(shù) Rx(m)(|m|=0,1,2,N-1), ,然后再經(jīng)過(guò)離散傅里葉變換求的功 率譜密度的估值 。 由式- = 估算出 ，再對(duì) 作 FFTFFT 變換，得到 - -2 23030 - 精品文檔 10歡迎下載 平均法精品文檔 11歡迎下載 平均法可視為周期圖法的改進(jìn)。周期圖經(jīng)過(guò)平均后會(huì)使它的方差減少，達(dá)到一致估計(jì)的 目的。 獲取 x(n)x(n)后，將 x(n)x(n)分為 1010 段，對(duì)每段用 計(jì)算出周期圖，對(duì)以上 1010 個(gè)周期圖加以平均得出功率譜。 Welch 法 WelchWelch 法又稱(chēng)修正周期法

14、。 獲取 x(n)x(n)后，先將 N N 個(gè)數(shù)據(jù)分成 L L 段，選擇漢明窗 w(n) )，并用該w(n) )依次對(duì)每段數(shù)據(jù) 做相應(yīng)的加權(quán)，確定每段的周期圖 1 L 、 由GM w GM , i w得出每段譜估計(jì)。 L i 二 GUI界面如下: GUIGUI 界面 五. 結(jié)果分析 由圖可以看出，在頻率 30hz30hz 和 100hz100hz 處功率譜有兩個(gè)峰值， 說(shuō)明信號(hào)中有 30hz30hz 和 100hz 100hz 的周期成分。GM ,i w 1 MU Ml J n=M(l J) 2 xnw(n)eQ ,I=1,2,、,L 經(jīng)典譜估計(jì) -ing 工抱瓦.郭敏 Iff s ft新 1

15、C h 二H.M惻剖日.呼岳 口 0 0 D D 4J?4J? 芒.蘭審也4 精品文檔 12 歡。迎下載 通過(guò) BT BT 法能觀察到兩個(gè)峰值， 但是所呈現(xiàn)的波形不能準(zhǔn)確表達(dá)出信號(hào)的功率譜變化情況。 通過(guò)周期圖法求出的功率譜密度在很大范圍內(nèi)波動(dòng)，而且容易證明，即使增加信號(hào)取樣 點(diǎn)數(shù) N,N,實(shí)驗(yàn)效果依然沒(méi)有明顯改進(jìn)。 用有限長(zhǎng)樣本序列的 DFTDFT 來(lái)表示隨機(jī)序列的功率譜只是一種估計(jì)或近似，不可避免存在 誤差。為了減少誤差，使功率譜估計(jì)更加平滑，可采用平均法。 平均法采用了分段的功率譜估計(jì)，較之于周期圖法，平均法的估計(jì)曲線(xiàn)較為平滑。在程 序中加入窗函數(shù)，使得譜分辨率不會(huì)由于分段而下降。 We

16、lchWelch 法就是用改進(jìn)的平均法來(lái)求取隨機(jī)信號(hào)的功率譜密度估計(jì)。 WelchWelch 法采用信號(hào)重疊 分段，加窗函數(shù)和 FFTFFT 算法等計(jì)算一個(gè)信號(hào)序列的功率譜。 六 . 任務(wù)分工 王龍?jiān)?(學(xué)號(hào)： 3832008009)完成了 BT 法以及周期圖法的程序編寫(xiě)以及 GUI 界面的設(shè)計(jì)和 PPT 的編寫(xiě)； 郭 敏 (學(xué)號(hào)： 3222008054)完成了平均法以及 Welch 法的程序編寫(xiě)以及 實(shí)驗(yàn)報(bào)告的統(tǒng)籌 , 編寫(xiě). 七 . 參考文獻(xiàn) 11 常建平 , , 李林海. . 隨機(jī)信號(hào)分析 . . 科學(xué)出版社 .2010 .2010 年，第 5 5 章第 3 3 節(jié)：1841841891

17、89； 22 黃文梅, , 熊桂林, , 楊勇. .信號(hào)分析與處理 . .國(guó)防科技大學(xué)出版社 .2000 .2000 年, ,第 6 6 章第 3 3 節(jié):222 :222 229229 八附錄 程序源代碼： function varargout = ming(varargin) gui_Singleton = 1； gui_State = struct(gui_Name, mfilename, . gui_Singleton, gui_Singleton, . gui_OpeningFcn, ming_OpeningFcn, . gui_OutputFcn, ming_OutputFcn,

18、. gui_LayoutFcn, , . gui_Callback, )； if nargin & ischar(varargin1) gui_State.gui_Callback = str2func(varargin1)； end if nargout varargout1:nargout = gui_mainfcn(gui_State, varargin:)； else gui_mainfcn(gui_State, varargin:)； end function ming_OpeningFcn(hObject, eventdata, handles, varargin) gui

19、data(hObject, handles); function varargout = ming_OutputFcn(hObject, eventdata, handles) varargout1 = handles.output; function figure1_CreateFcn(hObject, eventdata, handles) function pushbutton1_Callback(hObject, eventdata, handles) global y Fs nfft n; 精品文檔 13 歡。迎下載 Fs=600; nfft=1024; n=0:1/Fs:1; y=

20、cos(2*pi*30*n)+3*cos(2*pi*100*n)+randn(size(n); axes(handles.axes1); plot(n,y); title( 原始信號(hào)波形 ); function pushbutton3_Callback(hObject, eventdata, handles) global y Fs nfft n; window=hamming(nfft); noverlap=0; p=0.9; Pxx,Pxxc=psd(y,nfft,Fs,window,noverlap,p); index=0:round(nfft/2-1); k=index*Fs/nfft

21、; Qxx=10*log10(Pxx(index+1); axes(handles.axes6); plot(k,Qxx ); xlabel( 頻率( hz) ); ylabel( 功率譜密度( Db/Hz); title( 平均法 ); function pushbutton5_Callback(hObject, eventdata, handles) global y Fs nfft n; cxn=xcorr(y,unbiased); CXk=fft(cxn,nfft); Pxx=abs(CXk); index=0:round(nfft/2-1); k=index*Fs/nfft; plot_Pxx=10*log(Pxx(index+1)/(log(10); axes(handles.axes2); plot(k