lpc

Digital Signal Processing Library

Voice Lab

lpc

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 */