基于Burg算法的最大熵譜估計

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

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

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

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

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

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

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

8、-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); %遞推計算各階前向預(yù)測濾波器的系數(shù) endwhile FPE(M-1)TH %比較此刻階數(shù)的誤差平均功率與之前設(shè)置的平均功率的大小 TH=FPE(M-1); %小于之前數(shù)值時，覆蓋得到新的最小平均功率，并進行遞推運算 for n=M:N ef(M,n)=ef(M-1,n)+Km*eb(M-1,n-1); %遞推計算濾波器輸出的各階前向預(yù)測

9、誤差 eb(M,n)=eb(M-1,n-1)+Km*ef(M-1,n); %遞推計算濾波器輸出的各階后向預(yù)測誤差 end Km=-2*A/D; %反射系數(shù) a(M-1,M-1)=-2*A/D; %前向預(yù)測濾波器系數(shù) P(M)=P(M-1)*(1-(abs(Km)2); %預(yù)測誤差功率 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ù)測濾波器的系數(shù) endendT=1/fs;sum1=0; %采樣周期T賦值；功率譜初值設(shè)為0f=0.01:0.0

10、1:0.5; %選取數(shù)據(jù)采樣點，歸一化頻率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ù)的估計 ends1=(abs(1+sum1).2; %由Fejer-Riesz定理，得到最大熵譜估計，即ARMA功率譜s=P(M)*T./s1; %求得各階功率譜的矩陣 plot(f,10*log10(s),k); %畫出功率譜隨頻率變化的曲線圖xlabel(f/fs); %X軸坐標名稱ylabel(功率譜/dB); %Y軸坐標名稱五、 仿真結(jié)果及分析結(jié)果分析：如上，是階數(shù)分別為16和32時所得的功率譜曲線。如圖可知，Burg算法得到的譜線分辨率很高，曲線平