burg_ar.cpp

求AR模型系数的burg算法

#include   <stdio.h>
#include   <stdlib.h>


void Burg_AR(double *x,int N,int p,double *a,double *sigma2)
{
	int i,k;
	double R0,sumd,sumn,*b,*ef,*eb,*garma;
    b=(double*)malloc((p+1)*sizeof(double));
    ef=(double*)malloc(N*sizeof(double));
    eb=(double*)malloc(N*sizeof(double));
	garma=(double*)malloc((p+1)*sizeof(double));

	//初始化,0阶AR模型
    a[0]=1.0;
	b[0]=1.0;
	for (i=1;i<=p;i++)
	{
		a[i]=b[i]=0;
	}

    R0=0.0;   
    for   (i=0;i<N;i++)
	{
		R0+=x[i]*x[i]/N;
	}
	sigma2[0]=R0;
	
	for (i=0;i<N;i++)
	{
		ef[i]=x[i];
		eb[i]=x[i];
	}

	//逐阶计算

	for (k=1;k<=p;k++)
	{
		sumn=0.0;
		sumd=0.0;
		for (i=k;i<N;i++)
		{
			sumn+=ef[i]*eb[i-1];
			sumd+=ef[i]*ef[i]+eb[i-1]*eb[i-1];
		}
		garma[k]=-2*sumn/sumd;

		for (i=1;i<=k;i++)
		{
			b[i]=a[i]+garma[k]*a[k-i];//计算新的k阶AR系数
		}
		for (i=1;i<=k;i++)
		{
			a[i]=b[i];//新的k阶AR系数赋值
		}

		sigma2[k]=(1.0-garma[k]*garma[k])*sigma2[k-1];//白噪声功率

		//计算本阶AR模型的前向和后向预测误差
		//从后向前计算可以节约存储空间
		for (i=(N-1);i>=k;i--)
		{
			eb[i]=eb[i-1]+garma[k]*ef[i];
			ef[i]=ef[i]+garma[k]*eb[i-1];	
		}
	}

	free(b);
	free(ef);
	free(eb);
	free(garma);

}


void main()
{
	double sig[10]={1,2,3,78,45,23,-45,-12,2,1};
	int N=10;
	int i;
	double *a,*sigma2;
	int p=4;
	a=(double*)malloc((p+1)*sizeof(double));
	sigma2=(double*)malloc((p+1)*sizeof(double));

	Burg_AR(sig,N,p,a,sigma2);
for (i=0;i<=p;i++)
{
	printf("%f\n", a[i]);
	printf("%f\n\n", sigma2[i]);
}
	free(a);
	free(sigma2);

}