在随机信号分析中,可以用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)计算和;
扫描二维码关注公众号,回复:
176060 查看本文章
(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