基于Burg算法的最大熵譜估計(jì)

1、精選優(yōu)質(zhì)文檔-傾情為你奉上基于Burg算法的最大熵譜估計(jì)一、 實(shí)驗(yàn)?zāi)康氖褂肕atlab平臺(tái)實(shí)現(xiàn)基于Burg算法的最大熵譜估計(jì)二、 Burg算法原理現(xiàn)代譜估計(jì)是針對(duì)經(jīng)典譜估計(jì)方差性能較差、分辨率較低的缺點(diǎn)提出并逐漸發(fā)展起來(lái)的，其分為參數(shù)模型譜估計(jì)和非參數(shù)模型譜估計(jì)。而參數(shù)模型譜估計(jì)主要有AR模型、MA模型、ARMA模型等，其中AR模型應(yīng)用最多。ARMA模型功率譜的數(shù)學(xué)表達(dá)式為：其中，P(ej)為功率譜密度；s2是激勵(lì)白噪聲的方差；ai和bi為模型參數(shù)。若ARMA 模型中bi全為0，就變成了AR模型，又稱線性自回歸模型，其是一個(gè)全極點(diǎn)模型：研究表明，ARMA 模型和MA 模型均可用無(wú)限階的AR模型

2、來(lái)表示。且AR模型的參數(shù)估計(jì)計(jì)算相對(duì)簡(jiǎn)單。同時(shí)，實(shí)際的物理系統(tǒng)通常是全極點(diǎn)系統(tǒng)。要利用AR模型進(jìn)行功率譜估計(jì)，必須由Yule - Walker 方程求得AR模型的參數(shù)。而目前求解Yule - Walker方程主要有三種方法: Levinson-Durbin遞推算法、Burg算法和協(xié)方差方法。其中Burg算法計(jì)算結(jié)果較為準(zhǔn)確，且對(duì)于短的時(shí)間序列仍能得到較正確的估計(jì)，因此應(yīng)用廣泛。研究最大熵譜估計(jì)時(shí)，Levinson遞推一直受制于反射系數(shù)Km的求出。而B(niǎo)urg算法秉著使前、后向預(yù)測(cè)誤差平均功率最小的基本思想，不直接估計(jì)AR模型的參數(shù)，而是先估計(jì)反射系數(shù)Km，再利用Levinson關(guān)系式求得AR模型

3、的參數(shù)，繼而得到功率譜估計(jì)。Burg定義m階前、后向預(yù)測(cè)誤差為： (1) (2)由式（1）和（2）又可得到前、后預(yù)測(cè)誤差的階數(shù)遞推公式： (3) (4)定義m階前、后向預(yù)測(cè)誤差平均功率為： (5)將階數(shù)遞推公式（3）和（4）代入（5），并令，可得 (6)三、 Burg算法遞推步驟Burg算法的具體實(shí)現(xiàn)步驟：步驟1 計(jì)算預(yù)測(cè)誤差功率的初始值和前、后向預(yù)測(cè)誤差的初始值，并令m = 1。步驟2 求反射系數(shù)步驟3 計(jì)算前向預(yù)測(cè)濾波器系數(shù) 步驟4 計(jì)算預(yù)測(cè)誤差功率步驟5計(jì)算濾波器輸出步驟6 令m m+1，并重復(fù)步驟2至步驟5，直到預(yù)測(cè)誤差功率Pm不再明顯減小。最后，再利用Levinson遞推關(guān)系式估計(jì)A

4、R參數(shù)，繼而得到功率譜估計(jì)。四、 程序?qū)崿F(xiàn)%基于Burg算法的最大熵譜估計(jì)的Matlab實(shí)現(xiàn)%設(shè)置兩正弦小信號(hào)的歸一化頻率分別為0.175和0.20，信噪比SNR=30dB、N=32%clear,clc; %清空內(nèi)存及變量N=32; %設(shè)置離散傅里葉變換點(diǎn)數(shù)，即最大階數(shù)N為32 SNR=30; %信噪比SNR取為30dBfs=1; %采樣頻率取為1Hzt=1:N; %采樣時(shí)間點(diǎn)從1變化到Nt=t/fs; %得到歸一化頻率采樣點(diǎn)y=sin(2*pi*0.175*t)+sin(2*pi*0.20*t); %信號(hào)歸一化頻率分別取為0.175和0.20x=awgn(y,SNR); %在信號(hào)y中加入高斯

5、白噪聲,信噪比為SNR設(shè)定的數(shù)值 M=1; %設(shè)置起始計(jì)算的階數(shù)M為1P(M)=0; %預(yù)測(cè)誤差功率初值設(shè)為0Rx(M)=0; %自相關(guān)函數(shù)初值設(shè)為0for n=1:N %樣本數(shù)從1變化到N P(M)=P(M)+(abs(x(n)2; %計(jì)算預(yù)測(cè)誤差功率和的初始值 ef(1,n)=x(n); %計(jì)算前向預(yù)測(cè)誤差初值，令其等于此時(shí)的信號(hào)序列 eb(1,n)=x(n); %計(jì)算后向預(yù)測(cè)誤差初值，令其等于此時(shí)的信號(hào)序列endP(M)=P(M)/N; %計(jì)算出預(yù)測(cè)誤差功率的初始值 Rx(M)=P(M); %設(shè)定自相關(guān)函數(shù)初始值 M=2; %設(shè)置起始計(jì)算的階數(shù)M為2A=0; %微分所得反射系數(shù)Km的分子

6、，初始值設(shè)為0D=0; %微分所得反射系數(shù)Km的分母，初始值設(shè)為0for n=M:N %AR階數(shù)由M變化到N A=A+ef(M-1,n)*eb(M-1,n-1); %計(jì)算分子的和 D=D+(abs(ef(M-1,n)2+(abs(eb(M-1,n-1)2; %計(jì)算分母的和（即M階前、后向預(yù)測(cè)誤差平均功率)endKm=-2*A/D; %計(jì)算反射系數(shù)Km（此時(shí)起始階數(shù)為2）a(M-1,M-1)=-2*A/D; %計(jì)算前向預(yù)測(cè)濾波器系數(shù)P(M)=P(M-1)*(1-(abs(Km)2); %計(jì)算預(yù)測(cè)誤差功率FPE(M-1)=P(M)*(N+M)/(N-M); %設(shè)置最大預(yù)測(cè)誤差平均功率 TH=FPE

7、(M-1); for n=M:N %AR階數(shù)由M變化到N ef(M,n)=ef(M-1,n)+Km*eb(M-1,n-1); %計(jì)算濾波器輸出的前向預(yù)測(cè)誤差 eb(M,n)=eb(M-1,n-1)+Km*ef(M-1,n); %計(jì)算濾波器輸出的后向預(yù)測(cè)誤差 end M=M+1; %階數(shù)疊加，以便遞推計(jì)算下一階數(shù)據(jù) A=0; %反射系數(shù)Km的分子，初始值設(shè)為0 D=0; %反射系數(shù)Km的分母，初始值設(shè)為0 for n=M:N %同前，進(jìn)行遞推運(yùn)算 A=A+ef(M-1,n)*eb(M-1,n-1); D=D+(abs(ef(M-1,n)2+(abs(eb(M-1,n-1)2; end Km=-2

8、*A/D; a(M-1,M-1)=-2*A/D; P(M)=P(M-1)*(1-(abs(Km)2); FPE(M-1)=P(M)*(N+M)/(N-M); for m=1:M-2 %AR階數(shù)m由1變化到（M-2） a(M-1,m)=a(M-2,m)+Km*a(M-2,M-1-m); %遞推計(jì)算各階前向預(yù)測(cè)濾波器的系數(shù) endwhile FPE(M-1)<TH %比較此刻階數(shù)的誤差平均功率與之前設(shè)置的平均功率的大小 TH=FPE(M-1); %小于之前數(shù)值時(shí)，覆蓋得到新的最小平均功率，并進(jìn)行遞推運(yùn)算 for n=M:N ef(M,n)=ef(M-1,n)+Km*eb(M-1,n-1);

9、%遞推計(jì)算濾波器輸出的各階前向預(yù)測(cè)誤差 eb(M,n)=eb(M-1,n-1)+Km*ef(M-1,n); %遞推計(jì)算濾波器輸出的各階后向預(yù)測(cè)誤差 end Km=-2*A/D; %反射系數(shù) a(M-1,M-1)=-2*A/D; %前向預(yù)測(cè)濾波器系數(shù) P(M)=P(M-1)*(1-(abs(Km)2); %預(yù)測(cè)誤差功率 FPE(M-1)=P(M)*(N+M)/(N-M); for m=1:M-2 %AR階數(shù)m由1變化到（M-2） a(M-1,m)=a(M-2,m)+Km*a(M-2,M-1-m); %遞推得到各階前向預(yù)測(cè)濾波器的系數(shù) endendT=1/fs;sum1=0; %采樣周期T賦值；功

10、率譜初值設(shè)為0f=0.01:0.01:0.5; %選取數(shù)據(jù)采樣點(diǎn)，歸一化頻率0.010.5，間隔為0.01 for m=1:M-1; %AR階數(shù)為1(M-1) sum1=sum1+a(M-1,m)*exp(-j*2*pi*m*f*T); %傅里葉變換，得到AR參數(shù)的估計(jì) ends1=(abs(1+sum1).2; %由Fejer-Riesz定理，得到最大熵譜估計(jì)，即ARMA功率譜s=P(M)*T./s1; %求得各階功率譜的矩陣 plot(f,10*log10(s),'k'); %畫(huà)出功率譜隨頻率變化的曲線圖xlabel('f/fs'); %X軸坐標(biāo)名稱ylabel('功率譜/dB'); %Y軸坐標(biāo)名稱五、 仿真結(jié)果及分析結(jié)果分析：如上，是階數(shù)分別為16和32時(shí)所得的功率譜曲線。如圖可知，Burg算法得到的譜線分辨率很高，曲線平滑而且干擾