matlab 最大熵谱估计,基于Burg算法的最大熵谱估计.doc

基于Burg算法的最大熵谱估计

基于Burg算法的最大熵谱估计

其中，P(ejω)为功率谱密度；s2是激励白噪声的方差；ai和bi为模型参数。

若ARMA 模型中bi全为0，就变成了AR模型，又称线性自回归模型，其是一个全极点模型：

研究表明，ARMA 模型和MA 模型均可用无限阶的AR模型来表示。且AR模型的参数估计计算相对简单。同时，实际的物理系统通常是全极点系统。

要利用AR模型进行功率谱估计，必须由Yule - Walker 方程求得AR模型的参数。而目前求解Yule - Walker方程主要有三种方法: Levinson-Durbin递推算法、Burg算法和协方差方法。其中Burg算法计算结果较为准确，且对于短的时间序列仍能得到较正确的估计，因此应用广泛。

研究最大熵谱估计时，Levinson递推一直受制于反射系数Km的求出。而Burg算法秉着使前、后向预测误差平均功率最小的基本思想，不直接估计AR模型的参数，而是先估计反射系数Km，再利用Levinson关系式求得AR模型的参数，继而得到功率谱估计。

Burg定义m阶前、后向预测误差为：

(1)

(2)

由式(1)和(2)又可得到前、后预测误差的阶数递推公式：

(3)

(4)

定义m阶前、后向预测误差平均功率为：

(5)

将阶数递推公式(3)和(4)代入(5)，并令，可得 (6)

Burg算法递推步骤

Burg算法的具体实现步骤：

步骤1 计算预测误差功率的初始值和前、后向预测误差的初始值，并令m = 1。

步骤2 求反射系数

步骤3 计算前向预测滤波器系数

步骤4 计算预测误差功率

步骤5计算滤波器输出

步骤6 令m ← m+1，并重复步骤2至步骤5，直到预测误差功率Pm不再明显减小。

最后，再利用Levinson递推AR参数，继而得到功率谱估计。

程序实现

%%%%%%%%%%%%基于Burg算法的最大熵谱估计的Matlab实现end

P(M)=P(M)/N; %计算出预测误差功率的初始值

Rx(M)=P(M); %设定自相关函数初始值

M=2; %设置起始计算的阶数M为2

A=0; %微分所得反射系数Km的分子，初始值设为0

D=0; %微分所得反射系数Km的分母，初始值设为0

for n=M:N %AR阶数由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阶前、后向预测误差平均功率)

end

Km=-2*A/D; %计算反射系数Km(此时起始阶数为2)

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); %设置最大预测误差平均功率

TH=FPE(M-1);

for n=M:N %AR阶数由M变化到N

ef(M,n)=ef(M-1,n)+Km*eb(M-1,n-1); %计算滤波器输出的前向预测误差

eb(M,n)=eb(M-1,n-1)+Km*ef(M-1,n); %计算滤波器输出的后向预测误差

end

M=M+1; %阶数叠加，以便递推计算下一阶数据

A=0; %反射系数Km的分子，初始值设为0

D=0; %反射系数Km的分母，初始值设为0

for n=M:N %同前，进行递推运算