burg.c

常用算法与数据结构原代码

#include "stdlib.h"
#include "stdio.h"

double uniform(double a, double b, long *s)
{
	double t;
	
	*s = 2045 * (*s) +1;
	*s = *s - (*s / 1048576) * 1048576;
	t  = (*s) / 1048576.0;
	t  = a + (b - a) * t;
	
	return t;
}

double gauss(double mean, double sigma, long *s)
{
	int	   i;
	double x, y;
	
	for (x=0, i=0; i<12; i++)
	{
		x += uniform(0.0, 1.0, s);
	}
	x -= 6.0;
	y =  mean + x * sigma;

	return y;
}

void arma(double a[], double b[], int p, int q, double mean,
		  double sigma, long *seed, double x[], int n)
{
	int i, k, m;
	double s, *w;

	w = malloc(n * sizeof(double));
	for (k=0; k<n; k++)
	{
		w[k] = gauss(mean, sigma, seed);
	}

	x[0] = b[0] * w[0];
	for (k=1; k<=p; k++)
	{
		s = 0.0;
		for (i=1; i<=k; i++)
		{
			s += a[i] * x[k-i];
		}
		s = b[0] * w[k] - s;
		if (q == 0)
		{
			x[k] =s;
			continue;
		}
		m = (k>q) ? q : k;
		for (i=1; i<=m; i++)
		{
			s += b[i] * w[k-i];
		}
		x[k] = s;
	}

	for (k=p+1; k<n; k++)
	{
		s = 0.0;
		for (i=1; i<=p; i++)
		{
			s += a[i] * x[k-i];
		}
		s = b[0] * w[k] -s;
		if (q == 0)
		{
			x[k] =s;
			continue;
		}
		for (i=1; i<=q; i++)
		{
			s += b[i] * w[k-i];
		}
		x[k] = s;
	}
	free(w);
}

void burg(double x[], int n, int p, double a[], double *v)
{
	int i, k;
	double r0, sumd, sumn, *b, *ef, *eb;

	b = malloc((p+1)*sizeof(double));
	ef = malloc(n*sizeof(double));
	eb = malloc(n*sizeof(double));
	
	a[0] =1.0;
	r0 = 0.0;
	for (i=0; i<n; i++)
	{
		r0 += x[i]*x[i]/n;
	}
	v[0]=r0;

	for(i=1; i<n; i++)
	{
		ef[i] = x[i];
		eb[i-1] = x[i-1];
	}
	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];
		}
		a[k] = -2 * sumn/sumd;
		for (i=1; i<k; i++)
		{
			b[i] = a[i] + a[k]*a[k-i];
		}
		for (i=1;i<k; i++)
		{
			a[i] = b[i];
		}
		v[0] = (1.0-a[k]*a[k])*v[0];
		for (i=n-1; i>=k+1; i--)
		{
			ef[i] = ef[i]+ a[k]*eb[i-1];
			eb[i-1] = eb[i-2] + a[k]*ef[i-1];
		}
	}
	free(b);
	free(ef);
	free(eb);
}

void main()
{
	double a[6] = {1.0, 0.1, -0.47, -0.15, -0.007, 0.051};
	double b[1] = {1.0};
	int p= 5;
	int q = 0;
	long seed = 13579;
	double mean = 0.0;
	double var = 1.0;
	double x[1000];
	double c[20];
	int i, n = 1000;
	double v;

	arma(a,b,p,q,mean,var, &seed, x,n);
	burg(x,n,p,c,&v);
	for (i=0;i<=p;i++)
	{
		printf("a(%d) = %10.7lf\n",i,c[i]);
	}
}