Division by zero in Levinson–Durbin recursion? I am doing a LPC analysis of a speech signal using the autocorrelation method. To calculate the LPCs I am using the Levinson–Durbin recursion. In my literature the error is initialized with the first element of the autocorrelation sequence. The algorithm later uses the error (from the previous pass) to calculate the reflection coefficient. However, the error is used as divisor in a division. Thus, in case the first element of the autocorrelation sequence is zero this leads to a problem in my C++ implementation of LPC analysis.

I have seen implementations setting all LPCs to zero in case the first element of the autocorrelation sequence is zero. Other implementations ensure that the first element of the autocorrelation is always at least 1.
Matlab gives me LPCs of NaN if the first element of the autocorrelation is zero in a call to its levinson function.

What is the "correct" way to cope with this divison-by-zero-problem? Or what is the most suitable for LPC analysis of a speech signal?
Matt, this is a concise snippet of code from long ago:

/**********************************************************************
*                                                                     *
* Levinson-Durbin Algorithm for LPC Coefficients                      *
*                                                                     *
*   input:                                                            *
*   n  = LPC order                                                    *
*   r -&gt; r[0] to r[n]    autocorrelation values                       *
*                                                                     *
*   output:                                                           *
*   a -&gt; a[0] to a[n]    LPC coefficients, a[0] = 1.0                 *
*   k -&gt; k[0] to k[n]    reflection coefficients, k[0] = unused       *
*                                                                     *
*                                                                     *
**********************************************************************/

#define N   17

levinson_durbin(const int n, const double* r, double* a, double* k)
  {
  double a_temp[N], alpha, epsilon;         /* n &lt;= N = constant  */
  int i, j;

  k[0]      = 0.0;                          /* unused */

  a[0]      = 1.0;
  a_temp[0] = 1.0;                          /* unnecessary but consistent */

  alpha     = r[0];

  for(i=1; i&lt;=n; i++)
      {
      epsilon = r[i];                       /* epsilon = a[0]*r[i]; */
      for(j=1; j&lt;i; j++)
          {
          epsilon += a[j]*r[i-j];
          }

      a[i] = k[i] = -epsilon/alpha;
      alpha = alpha*(1.0 - k[i]*k[i]);

      for(j=1; j&lt;i; j++)
          {
          a_temp[j] = a[j] + k[i]*a[i-j];   /* update a[] array into temporary array */
          }
      for(j=1; j&lt;i; j++)
          {
          a[j] = a_temp[j];                 /* update a[] array */
          }
      }
  }


now i see many places (well, one place, but many instances) where division by zero is possible.  any time k[i] = 1 or -1.  can't that conceivably happen with a non-zero input?