float vorbis_lpc_from_data()

in nlsCppSdk/encoder/lpc.cpp [58:157]

• 74 lines of code
• 16 McCabe index (conditional complexity)

float vorbis_lpc_from_data(float *data, float *lpci, int n, int m, int stride) {
  double *aut = (double *)malloc(sizeof(*aut) * (m + 1));
  if (NULL == aut) {
    return 0;
  } else {
    memset(aut, 0, sizeof(*aut) * (m + 1));
  }
  double *lpc = (double *)malloc(sizeof(*lpc) * (m));
  if (!lpc) {
    free(aut);
    aut = NULL;
    return 0;
  } else {
    memset(lpc, 0, sizeof(*lpc) * m);
  }
  double error;
  double epsilon;
  int i, j;

  if (!aut || !lpc) {
    if (aut) {
      free(aut);
      aut = NULL;
    }
    if (lpc) {
      free(lpc);
      lpc = NULL;
    }
    return 0;
  }

  /* autocorrelation, p+1 lag coefficients */
  j = m + 1;
  while (j--) {
    double d = 0; /* double needed for accumulator depth */
    for (i = j; i < n; i++) {
      d += (double)data[i * stride] * data[(i - j) * stride];
    }
    aut[j] = d;
  }

  /* Generate lpc coefficients from autocorr values */

  /* set our noise floor to about -100dB */
  error = aut[0] * (1. + 1e-10);
  epsilon = 1e-9 * aut[0] + 1e-10;

  for (i = 0; i < m; i++) {
    double r = -aut[i + 1];

    if (error < epsilon) {
      memset(lpc + i, 0, (m - i) * sizeof(*lpc));
      goto done;
    }

    /* Sum up this iteration's reflection coefficient; note that in
       Vorbis we don't save it.  If anyone wants to recycle this code
       and needs reflection coefficients, save the results of 'r' from
       each iteration. */

    for (j = 0; j < i; j++) r -= lpc[j] * aut[i - j];
    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 & 1) {
      lpc[j] += lpc[j] * r;
    }

    error *= 1. - r * r;
  }

done:

  /* slightly damp the filter */
  {
    double g = .99;
    double damp = g;
    for (j = 0; j < m; j++) {
      lpc[j] *= damp;
      damp *= g;
    }
  }

  for (j = 0; j < m; j++) {
    lpci[j] = (float)lpc[j];
  }

  /* we need the error value to know how big an impulse to hit the
     filter with later */
  free(aut);
  free(lpc);
  return error;
}