# 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?

• The first element of the autocorrelation (i.e. its value at lag zero) can only be zero for a signal consisting of zeros only, so you just need to check if the input signal is an all-zero signal. Jun 23 '15 at 15:31

Matt, this is a concise snippet of code from long ago:

/**********************************************************************
*                                                                     *
* Levinson-Durbin Algorithm for LPC Coefficients                      *
*                                                                     *
*   input:                                                            *
*   n  = LPC order                                                    *
*   r -> r[0] to r[n]    autocorrelation values                       *
*                                                                     *
*   output:                                                           *
*   a -> a[0] to a[n]    LPC coefficients, a[0] = 1.0                 *
*   k -> 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 <= 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<=n; i++)
{
epsilon = r[i];                       /* epsilon = a[0]*r[i]; */
for(j=1; j<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<i; j++)
{
a_temp[j] = a[j] + k[i]*a[i-j];   /* update a[] array into temporary array */
}
for(j=1; j<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?

• It shouldn't, because as far as I remember the resulting filter is always stable, and stability means $|k|<1$ (i.e., strictly less than $1$). Of course, if $k$ is close enough to $1$ you might still run into numerical trouble, but the denominator should at least theoretically never become exactly zero. Jun 23 '15 at 18:30