#define	P_MAX	8	/* order p of LPC analysis, typically 8..14 */

double levinson_durbin(double const * ac, double * ref, double * lpc);
double schur(double const * ac, double * ref);
void autocorrelation(int n, double const * x, int lag, double * ac);

/* LPC- and Reflection Coefficients
 * 
 * The next two functions calculate linear prediction coefficients
 * and/or the related reflection coefficients from the first P_MAX+1
 * values of the autocorrelation function.
 */

/* Invented by N. Levinson in 1947, modified by J. Durbin in 1959.
 */
double            	    /* returns minimum mean square error    */
levinson_durbin(
	double const * ac,  /*  in: [0...p] autocorrelation values  */
	double 	     * ref, /* out: [0...p-1] reflection coef's	    */
	double	     * lpc) /*      [0...p-1] LPC coefficients      */
{
	int i, j;  double r, error = ac[0];

	if (ac[0] == 0) {
		for (i = 0; i < P_MAX; i++) ref[i] = 0; return 0; }

	for (i = 0; i < P_MAX; i++) {

		/* Sum up this iteration's reflection coefficient.
		 */
		r = -ac[i + 1];
		for (j = 0; j < i; j++) r -= lpc[j] * ac[i - j];
		ref[i] = r /= error;    

		/*  Update LPC coefficients and total error.
		 */
		lpc[i] = r;
		for (j = 0; j < i/2; j++) {
			double tmp  = lpc[j];
			lpc[j]     += r * lpc[i-1-j];
			lpc[i-1-j] += r * tmp;
		}
		if (i % 2) lpc[j] += lpc[j] * r;

		error *= 1.0 - r * r;
	}
	return error;
}

/* I. Sch"ur's recursion from 1917 is related to the Levinson-Durbin
 * method, but faster on parallel architectures; where Levinson-Durbin
 * would take time proportional to p * log(p), Sch"ur only requires
 * time proportional to p.  The GSM coder uses an integer version of
 * the Sch"ur recursion.
 */
double			   /* returns the minimum mean square error */
schur(  double const * ac,	/*  in: [0...p] autocorrelation c's */
	double       * ref)	/* out: [0...p-1] reflection c's    */
{
	int i, m;  double r, error = ac[0], G[2][P_MAX];

	if (ac[0] == 0.0) {
		for (i = 0; i < P_MAX; i++) ref[i] = 0;
		return 0;
	}

	/*  Initialize rows of the generator matrix G to ac[1...P]
	 */
	for (i = 0; i < P_MAX; i++) G[0][i] = G[1][i] = ac[i+1];

	for (i = 0;;) {

		/*  Calculate this iteration's reflection
		 *  coefficient and error.
		 */
		ref[i] = r = -G[1][0] / error;
		error += G[1][0] * r;

		if (++i >= P_MAX) return error;

		/*  Update the generator matrix.
		 *
		 *  Unlike the Levinson-Durbin summing of reflection
		 *  coefficients, this loop could be distributed to
		 *  p processors who each take only constant time.
		 */
		for (m = 0; m < P_MAX - i; m++) {
			G[1][m] = G[1][m+1] + r * G[0][m];
			G[0][m] = G[1][m+1] * r + G[0][m];
		}
	}
}

/* Compute the autocorrelation 
 *			,--,
 *              ac(i) = >  x(n) * x(n-i)  for all n
 *	       		`--'
 * for lags between 0 and lag-1, and x == 0 outside 0...n-1
 */
void autocorrelation(
	int   n, double const * x,   /*  in: [0...n-1] samples x   */
	int lag, double       * ac)  /* out: [0...lag-1] ac values */
{
	double d; int i;
	while (lag--) {
		for (i = lag, d = 0; i < n; i++) d += x[i] * x[i-lag];
		ac[lag] = d;
	}
}