Signal Processing Library

Program name: lpc

Language: C

In file: parameter_lib.c

Objective: Calculate the Linear Predictive Coefficients of the data in the window

Usage: void lpc(double *rr, double am[], int order, double *gain);

Parameters:

• rr - The autocorrelation coefficients for the window of input data
• am - The output array to put the calculated data
• order - The order
• gain - The calculated gain

• void

Mathematical Description:

Comments: Routine to calculate the LPC of the data in the window

User Comments

Code:

void lpc(double *rr, double am[], int order, double *gain)
{
	int i, m, n, k;
	/* Local arrays need to be at least order+1 size */
	double kk[MAX_ALLOC_SIZE]; /* k(m) on page 343 of book cited below */
	double e[MAX_ALLOC_SIZE];  /* E(m) on page 343 of book cited below */
	double previous_am[MAX_ALLOC_SIZE]; /* am-1(m) on page 343 of book cited below */
	double aux;

	/* Information taken from Douglas O'Shaughnessy's book */
	/*   "Speech Communication, Human and Machine" 1990, pages 341-344 */
	/* A few things to keep in mind while accessing information in the various arrays: */
	/* C is very flexible, but leaves the index limit check up to the programmer */
	/* Limits must be set up and tested so as not to go out of bounds */

	/* the input rr is defined for indexes varying from 0 to order */
	/* output array am[] should have info stored from index 0 thru p, this is what will be saved in outfile */
	/* am[0] = 1.0 */

	/* Test if space allocated for arrays is large enough */
	/* Local arrays need to be at least order+1 size */
	if(order >= MAX_ALLOC_SIZE) {
		fprintf(stderr, "Please increase MAX_ALLOC_SIZE in window_lib.h to at least %d\n", order+1);
		exit(1);
	}

	for(i = 0; i <= order; i++) {
		/* Start with a zero kk vector */
		kk[i] = 0.0;
		/* first "previous_am vector" set to zero */
		previous_am[i] = 0.0;
		/* first coefficient of LPC set to zero */
		am[i] = 0.0;
	}

	//aux = rr[0];

	/* Calculation of the LPC coefficients */
	if(rr[0] <= MINFLOAT) {
		/* no signal. Nothing to do */
	} else {
		e[0]=rr[0];
		am[0] = 1.0;
		/* Interval is 1 to order+1, but indexes in C are from 0 to order */
		/* Care must be taken to not go beyond end of arrays */
		for(m = 1; m <= order; m++) {
			/* equation 8.16.a on page 343 */
			aux=0.0;
			/* formula interval goes from k = 1 to k = m, */
			/* therefore k will always be less than order+1 and less than m+1 */
			for(k = 1; k < m; k++) {
				aux += previous_am[k]*rr[m-k];
			}
			kk[m] = (rr[m]-aux)/e[m-1];

			/* equation 8.16.b on page 343 */
			/* Indexes of am are in the interval between 0 and order */
			am[m] = kk[m];

			/* equation 8.16.c on page 343 */
			/* formula interval goes from k = 1 to k = m - 1, */
			/* therefore k will always be less than order */
			for(k = 1; k < m; k++) {
				/* Indexes of am are in the interval between 1 and order+1 */
				am[k] = previous_am[k] - kk[m]*previous_am[m-k];
			}

			/* equation 8.16.d on page 343 */
			e[m] = (1.0-kk[m]*kk[m])*e[m-1];

			/* save current vector as previous_am vector */
			/* Indexes of am are in the interval between 1 and order+1 */
			for(k = 1; k <= m; k++) {
				previous_am[k] = am[k];
			}
		}
	}
	aux = rr[0];
	for(m = 1; m <= order; m++) {
		aux -= am[m]*rr[m]; /* Equation 8.10 page 341 */
	}
	*gain = sqrt(aux);
} /* end of the Levinson-Durbin LPC routine */