/*		Levnsn.c		*/
/* Levinson-Durbin LPC
 *
 * | R0 R1 R2 ... RN-1 |   | A1 |       | -R1 |
 * | R1 R0 R1 ... RN-2 |   | A2 |       | -R2 |
 * | R2 R1 R0 ... RN-3 |   | A3 |   =   | -R3 |
 * |          ...      |   | ...|       | ... |
 * | RN-1 RN-2... R0   |   | AN |       | -RN |
 *
 * Ref: John Makhoul, "Linear Prediction, A Tutorial Review"
 * Proc. IEEE Vol. 63, PP 561-580 April, 1975.
 *
 * R is the input autocorrelation function.  R0 is the zero lag
 * term.  A is the output array of predictor coefficients.  Note
 * that a filter impulse response has a coefficient of 1.0 preceding
 * A1.  E is an array of mean square error for each prediction order
 * 1 to N.  REFL is an output array of the reflection coefficients.
 */

#define abs(x) ( (x) < 0 ? -(x) : (x) )

int levnsn( n, r, a, e, refl )
int n;
double r[], a[], e[], refl[];
{
int k, km1, i, kmi, j;
double ai, akk, err, err1, r0, t, akmi;
double *pa, *pr;

for( i=0; i<n; i++ )
	{
	a[i] = 0.0;
	e[i] = 0.0;
	refl[i] = 0.0;
	}
r0 = r[0];
e[0] = r0;
err = r0;

akk = -r[1]/err;
err = (1.0 - akk*akk) * err;
e[1] = err;
a[1] = akk;
refl[1] = akk;

if( err < 1.0e-2 )
	return 0;

for( k=2; k<n; k++ )
	{
	t = 0.0;
	pa = &a[1];
	pr = &r[k-1];
	for( j=1; j<k; j++ )
		t += *pa++ * *pr--;
	akk = -( r[k] + t )/err;
	refl[k] = akk;
	km1 = k/2;
	for( j=1; j<=km1; j++ )
		{
		kmi = k-j;
		ai = a[j];
		akmi = a[kmi];
		a[j] = ai + akk*akmi;
		if( i == kmi )
			goto nxtk;
		a[kmi] = akmi + akk*ai;
		}
nxtk:
	a[k] = akk;
	err1 = (1.0 - akk*akk)*err;
	e[k] = err1;
	if( err1 < 0 )
		err1 = -err1;
/*	err1 = abs(err1);*/
/*	if( (err1 < 1.0e-2) || (err1 >= err) )*/
	if( err1 < 1.0e-2 )
		return 0;
	err = err1;
	}
  return 0;
}