Burg法求解AR(p)模型参数及MATLAB实现

    在随机信号分析中，可以用AR模型进行功率谱估计。在求解Yule-Walker方程中的AR系数可用Levinson递推算法简化计算，但它需要知道自相关序列。自相关序列实际上只能从随机序列x(n)的有限个观测数据估计得到。当时间序列较短时，的估计误差很大，这将对AR参数的计算引入较大的误差，导致谱估计性能下降，甚至出现谱线分裂与谱峰偏移等现象。如果利用观测到的数据直接计算AR模型的参数，则能克服上述方法的缺点，得到性能较好的谱估计结果。这种方法是由Burg提出的，称为Burg法，也叫做最大熵谱算法。

    Burg递推算法的优点是不需要估计自相关函数，它直接从已知序列x(n)求得反射系数，然后利用Levinson递推算法由反射系数来求得AR参数。

    Burg法估计AR(p)模型参数的具体计算步骤可归纳如下：

    （1）确定初始条件：

        （1-1）

        （1-2）

    （2），根据公式（1-3）计算（反射系数）；

        （1-3）

    （3）根据式（1-4）计算；

        （1-4）

    （4）根据式（1-5）和（1-6）计算和;

        （1-5）

        （1-6）

    （5）根据，计算；

    （6），重复（2）（3）（4）直至所需要的阶数为止。

    用MATLAB实现Burg算法如下所示：

function [fi,sita]=burg(data,p,q)
%data:原始数据
%p,q:arma模型的阶数,这里简化为q=1

e=zeros(p+1,length(data));
b=zeros(p+1,length(data));
sita=zeros(1,p+1);

N=length(data);
%1.确定初始条件
fi=zeros(p,p);
e(1,:)=data;
b(1,:)=data;
sita(1,1)=mean(data.^2);

for k=1:p
    %2.计算Kp 式5-119
    if(k==1)
        kp=-2*(sum(e(k,k+1:N).*(b(k,(k:N-1)))))/sum(e(k,k+1:N).^2+b(k,k:N-1).^2);
    else
        kp=-2*(sum(e(k,k:N).*(b(k,k-1:N-1))))/sum(e(k,k:N).^2+b(k,k-1:N-1).^2);
    end
    %3.计算aki
    for i=1:k-1 %式5-103
        fi(k,i)=fi(k-1,i)+kp*fi(k-1,k-i);
    end
    fi(k,k)=kp;
    %4.计算ek和bk
    for i=1:N
        if(i==1)
            e(k+1,i)=e(k,i);
            b(k+1,i)=kp*e(k,i);
        else
            e(k+1,i)=e(k,i)+kp*b(k,i-1);
            b(k+1,i)=b(k,i-1)+kp*e(k,i);
        end
    end
    %5.计算sita
    sita(k+1)=(1-kp.^2)*sita(k);
    %6.K=K+1
end