burg法實(shí)現(xiàn)功率譜估計(jì).doc

《burg法實(shí)現(xiàn)功率譜估計(jì).doc》由會(huì)員上傳分享，免費(fèi)在線閱讀，更多相關(guān)內(nèi)容在教育資源-天天文庫(kù)。

1、用burg法實(shí)現(xiàn)功率譜估計(jì)參數(shù)模型法是現(xiàn)代譜估計(jì)中的主要內(nèi)容，AR模型參數(shù)的求解有三種方法：自相關(guān)法、Burg遞推算法和改進(jìn)協(xié)方差法。Burg算法不是直接估計(jì)AR模型的參數(shù)，而是先估計(jì)反射系數(shù)Km,再利用Levinson關(guān)系式求得AR模型的參數(shù)。Burg算法采用的數(shù)據(jù)加窗方法是協(xié)方差法，不含有對(duì)已知數(shù)據(jù)段之外的數(shù)據(jù)做人為的假設(shè)。1.其原理如下：Burg算法是使前向預(yù)測(cè)誤差和后向預(yù)測(cè)誤差均方誤差之和最小來(lái)求取Km的，它不對(duì)已知數(shù)據(jù)段之外的數(shù)據(jù)做認(rèn)為假設(shè)。計(jì)算m階預(yù)測(cè)誤差的遞推表示公式如下：求取反射系數(shù)的公式如下：對(duì)于平穩(wěn)隨機(jī)過(guò)程

2、，可以用時(shí)間平均代替集合平均，因此上式可寫成：這樣便可求得AR模型的反射系數(shù)。將m階AR模型的反射系數(shù)和m-1階AR模型的系數(shù)代入到Levinson關(guān)系式中，可以求得AR模型其他的p-1個(gè)參數(shù)。Levinson關(guān)系式如下：m階AR模型的第m+1個(gè)參數(shù)G，,其中是預(yù)測(cè)誤差功率，可由遞推公式求得。易知為進(jìn)行該式的遞推，必須知道0階AR模型誤差功率=可知該式由給定序列易于求得。完成上述過(guò)程，即最終求得了表征該隨機(jī)信號(hào)的AR模型的p+1個(gè)參數(shù)。然后根據(jù)=即可求得該隨機(jī)信號(hào)的功率譜密度。一.數(shù)據(jù)為：實(shí)現(xiàn)功率譜估計(jì)的matlan源代碼：c

3、learFs=500;n=0:1/Fs:0.5;w1=200*pi;A1=5;w2=300*pi;A2=12;xn=A1*(cos(w1*n)+j*sin(w1*n))+A2*(j*sin(w2*n)+cos(w2*n))+randn(size(n));%xn=A1*exp(jw1n)+A2*exp(jw2n);subplot(211);plot(n,xn);xlabel('n');ylabel('xn');title('xn=A1*exp(jw1n)+A2*exp(jw2n)+e(n)');ymax_xn=max(xn)+0

4、.2;ymin_xn=min(xn)-0.2;axis([00.5ymin_xnymax_xn]);p=floor(length(xn)/3)+2;nfft=1024;[xpsd,f]=pburg(xn,p,nfft,Fs);pmax=max(xpsd);xpsd=xpsd/pmax;xpsd=10*log10(xpsd);subplot(212);plot(f,xpsd);title('PowerSpectralestimatewithburg');ylabel('PowerSpectralestimate(dB)');xl

5、abel('f(Hz)');gridon;ymin_psd=min(xpsd)-1;ymax_psd=max(xpsd)+1;axis([0Fs/2ymin_psdymax_psd]);實(shí)驗(yàn)結(jié)果：（1）下圖依次是階數(shù)為N/3，N/2二.數(shù)據(jù)為：Matlab源代碼：clearclf'clcN=100;Fs=500;%產(chǎn)生x1信號(hào)即AR信號(hào)vn=rands(1,N);xn=zeros(1,N);xn(1)=vn(1);xn(2)=vn(2);a1=0.78;a2=-0.96;forn=3:Nxn(n)=vn(n)+a1*xn(n-

6、1)+a2*xn(n-2)endsubplot(211);plot(xn);title('xn(n)=noise(n)+a1*xn(n-1)+a2*xn(n-2)');p=floor(length(xn)/5);nfft=1024;[xpsd,f]=pburg(xn,p,nfft,Fs);pmax=max(xpsd);xpsd=xpsd/pmax;xpsd=10*log10(xpsd);subplot(212);plot(f,xpsd);title('PowerSpectralestimatewithburg');ylabel

7、('PowerSpectralestimate(dB)');xlabel('f(Hz)');gridon;ymin_psd=min(xpsd)-1;ymax_psd=max(xpsd)+1;axis([0Fs/2ymin_psdymax_psd]);