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

本文主要是介绍Burg法求解AR(p)模型参数及MATLAB实现，希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！

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

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

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

（1）确定初始条件：

$e_{0}(n)=b_{0}(n)=x(n),0\leq n\leq N-1$ （1-1）

$\sigma _{0}^{2}=\frac{1}{N}\sum_{n=0}^{N-1}x^{2}(n)$ （1-2）

（2） $k=1$ ，根据公式（1-3）计算 $K_{k}$ （反射系数）；

$K_{k}=a_{kk}=-\frac{2\sum_{n=p}^{N-1}[e_{p-1}(n)b_{p-1}(n-1)]}{\sum_{n=p}^{N-1}[e_{p-1}^{2}(n)+b_{p-1}^{2}(n-1)]}$ （1-3）

（3）根据式（1-4）计算 $a_{ki}(i=1,2,...,k-1)$ ；

$a_{ki}=a_{k-1,i}+a_{kk}a_{k-1,k-i}$ （1-4）

（4）根据式（1-5）和（1-6）计算 $e_{k}(n)$ 和 $b_{k}(n)$ ;

$e_{k}(n)=e_{p-1}(n)+K_{p}b_{p-1}(n-1)$ （1-5）

$b_{k}(n)=b_{p-1}(n)+K_{p}e_{p-1}(n)$ （1-6）

（5）根据 $\sigma _{k}^{2}=(1-K_{p}^{2})\sigma _{k-1}^{2}$ ，计算 $\sigma _{k}^{2}$ ；

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

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

function [fi,sita]=burg(data,p,q)
%data:原始数据
%p,q:arma模型的阶数,这里简化为q=1e=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-119if(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);elsekp=-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.计算akifor i=1:k-1 %式5-103fi(k,i)=fi(k-1,i)+kp*fi(k-1,k-i);endfi(k,k)=kp;%4.计算ek和bkfor i=1:Nif(i==1)e(k+1,i)=e(k,i);b(k+1,i)=kp*e(k,i);elsee(k+1,i)=e(k,i)+kp*b(k,i-1);b(k+1,i)=b(k,i-1)+kp*e(k,i);endend%5.计算sitasita(k+1)=(1-kp.^2)*sita(k);%6.K=K+1
end

这篇关于Burg法求解AR(p)模型参数及MATLAB实现的文章就介绍到这儿，希望我们推荐的文章对编程师们有所帮助！