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
Return value:
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 */
|