#include <stdio.h>
#include <math.h>
#include <alloc.h>
#define N 300

typedef double Real;

Real PI = 3.14159265358979323846;

/*---------- MEM2.C ----------*/
extern void levins(int nmax, Real *g, Real *pm, Real *b1 , Real *b2 ,
					Real *a , int m);
extern void mk_dat(Real *x , int *nmax , int *mmx , int *lmax , Real *dt);
extern void burg(int nmax , int mmx , int lmax , Real *x , Real *g , Real *c ,
				Real *fpe , Real *pm);
extern void mem_1(int nmax , int mmx , Real dt , Real *g , Real *f ,
				Real *spec , Real pm);

Real x[N] , f[N] , fpe[N] , spec[N] , g[N] , c[N];

void main(void)
{
	int nmax , mmx , lmax , i , lag , lg1 , lx1;
	Real dt , sum , av , pm;
	
	mk_dat(x , &nmax , &mmx , &lmax , &dt);
	lag = mmx + 1;
	if((lag+1) > lmax) lmax = lag + 1;
	
	sum = 0.0;
	for(i = 1 ; i <= nmax ; i++) sum += x[i];
	av = sum / (Real)nmax;
	for(i = 1 ; i <= nmax ; i++) x[i] -= av;
	
	burg(nmax , mmx , lmax , x , g , c , fpe , &pm);
	mem_1(nmax , mmx , dt , g , f , spec , pm);
	
	lag = mmx + 1;
	lg1 = lag + 1;
	lx1 = lmax + 1;
	printf("\n frec      spc        fpe    auto c.");
	for(i = 1 ; i <= lag ; i++)
		printf("\n%6.1f    %7.3f   %7.3f   %7.3f" ,
				f[i] , spec[i] , fpe[i] , c[i]);
	if(lag != lmax)
		for(i = lg1 ; i <= lmax ; i++)
			printf("\n%6.1f    %7.3f        %6.3f" ,
					f[i] , spec[i] , c[i]);
	for(i = lx1 ; i <= nmax ; i++)
		printf("\n%6.1f    %7.3f" , f[i] , spec[i]);
}

void levins(int nmax, Real *g, Real *pm, Real *b1 , Real *b2 , Real *a , int m)
{
	int  i , k;
	Real sn , sd;
	
	sn = sd = 0.0;
	for(i = 1 ; i <= (nmax - m) ; i++)
	{
		sn += b1[i] * b2[i];
		sd += b1[i]*b1[i] + b2[i]*b2[i];
	}
	g[m] = -2.0 * sn / sd;
	*pm *= (1.0 - g[m] * g[m]);
	if(m != 1)
		for(k = 1 ; k <= m-1 ; k++) g[k] = a[k] + g[m] * a[m-k];
	for(i = 1 ; i <= m ; i++) a[i] = g[i];
	for(i = 1 ; i < nmax - m ; i++)
	{
		b1[i] += a[m] * b2[i];
		b2[i] = b2[i+1] + a[m] * b1[i+1];
	}
}

void burg(int nmax , int mmx , int lmax , Real *x , Real *g , Real *c ,
				Real *fpe , Real *pm)
{
	Real sum , *b1 , *b2 , *a;
	int  i , m , l , lag , lg1;
	
	b1 = (Real *)malloc(8*(nmax+1));
	b2 = (Real *)malloc(8*(nmax+1));
	a  = (Real *)malloc(8*(nmax+1));

	sum = 0.0;
	for(i = 1 ; i <= nmax ; i++) sum += x[i] * x[i];
	c[1] = sum / (Real)nmax;
	*pm = c[1];
	fpe[1] = (Real)(nmax+1)/(Real)(nmax-1) * (*pm);
	lag = mmx + 1;
	lg1 = lag + 1;
	b1[1] = x[1];
	for(i = 2 ; i <= nmax ; i++)
	{
		b1[i] = x[i];
		b2[i-1] = x[i];
	}
	
	for(m = 1 ; m <= mmx ; m++)
	{
		printf("\n M=%4d " , m);
		levins(nmax , g , pm , b1 , b2 , a , m);
		sum = 0.0;
		for(i = 1 ; i <= m ; i++) sum -= c[m-i+1] * g[i];
		c[m+1] = sum;
		if(m != (nmax-1))
			fpe[m+1] = (Real)(nmax+m+1) / (Real)(nmax-m-1) * (*pm);
	}
	
	if(lag < lmax)
		for(l = lg1 ; l <= lmax ; l++)
		{
			sum = 0.0;
			for(i = 1 ; i <= mmx ; i++) sum -= c[l-1] * g[i];
			c[l] = sum;
		}
	free((char *)b1);
	free((char *)b2);
	free((char *)a);
}

void mem_1(int nmax , int mmx , Real dt , Real *g , Real *f , Real *spec ,
					 Real pm)
{
	int i , j;
	Real f0 , sumr , sumi , cr , ci;
	
	f0 = 1.0 / ((Real)(nmax-1) * dt);
	for(i = 1 ; i <= nmax ; i++)
	{
		f[i] = f0 * (Real)(i-1);
		printf("\n  freq=%6.2f" , f[i]);
		sumr = 1.0;
		sumi = 0.0;
		cr = 0.0;
		ci = 2.0 * PI * f[i] * dt;
		for(j = 1 ; j <= mmx ; j++)
		{
			sumr += g[j] * cos((double)j * ci);
			sumi += g[j] * sin((double)j * ci);
		}
		spec[i] = 2.0 * pm * dt / (sumr*sumr + sumi*sumi);
	}
}

void mk_dat(Real *x , int *nmax , int *mmx , int *lmax , Real *dt)
{
	Real t1 , t2 , t3 , t4 , t , y;
	int  i;
	
	*nmax = 241;
	*mmx = 38;
	*lmax = 50;
	*dt = 1.0 / (Real)(*nmax-1);
	t1 = 1.0/124.5;
	t2 = 1.0/65.0;
	t3 = 1.0 / 64.0;
	t4 = 1.0/4.0;
	for(i = 1 ; i <= *nmax ; i++)
	{
		t = *dt * (Real)(i-1);
		y = sin(2.0 * PI * t / t1) + sin(2.0 * PI * t / t2)
			+ sin(2.0 * PI * t / t3) + sin(2.0 * PI * t / t4);
		x[i] = y;
		/* + 2.0 * drand48()  - 2;*/
	}
}