cfdlp.c
// levinson-durbin recursion outputting poles. also does gain normalization if
// requested
levinson-durbin recursion outputting poles. also does gain normalization if
requested
levinson - durbin recursion outputting poles . also does gain normalization if requested
void levinson(int p, double *phi, float *poles)
{
    double *alpha = (double *) MALLOC((p+1)*(p+1)*sizeof(double));
    double *E = (double *) MALLOC((p+1)*sizeof(double));
    double *k = (double *) MALLOC((p+1)*sizeof(double));
    float g;
    E[0] = phi[0];
    for ( int i = 1; i <= p; i++ )
    {
	k[i] = -phi[i];
	for ( int j = 1; j <= i-1; j++ )
	{
	    k[i] -= (phi[i-j] * alpha[(i-1)*(p+1)+j]);
	}
	k[i] /= E[i-1];
	alpha[i*(p+1)] = 1;
	alpha[i*(p+1)+i] = k[i];
	for ( int j = 1; j <= i-1; j++ )
	{
	    alpha[i*(p+1)+j] = alpha[(i-1)*(p+1)+j] + k[i]*alpha[(i-1)*(p+1)+i-j];
	}
	E[i] = (1-k[i]*k[i])*E[i-1];
    }
    g = sqrt(E[p]);
    for ( int i = 0; i <= p; i++ )
    {
	poles[i] = alpha[p*(p+1)+i];
	if (!do_gain_norm)
	{
	    if (E[p] <= 0.) {
		poles[i] = 0.;
	    }
	    else
	    {
		poles[i] /= g; 
	    }
	}	
    }
    FREE(k);
    FREE(E);
    FREE(alpha);
}
void levinson ( int p , double * phi , float * poles ) { double * alpha = ( double * ) MALLOC ( ( p + 1 ) * ( p + 1 ) * sizeof ( double ) ) ; double * E = ( double * ) MALLOC ( ( p + 1 ) * sizeof ( double ) ) ; double * k = ( double * ) MALLOC ( ( p + 1 ) * sizeof ( double ) ) ; float g ; E [ 0 ] = phi [ 0 ] ; for ( int i = 1 ; i <= p ; i ++ ) { k [ i ] = - phi [ i ] ; for ( int j = 1 ; j <= i - 1 ; j ++ ) { k [ i ] -= ( phi [ i - j ] * alpha [ ( i - 1 ) * ( p + 1 ) + j ] ) ; } k [ i ] /= E [ i - 1 ] ; alpha [ i * ( p + 1 ) ] = 1 ; alpha [ i * ( p + 1 ) + i ] = k [ i ] ; for ( int j = 1 ; j <= i - 1 ; j ++ ) { alpha [ i * ( p + 1 ) + j ] = alpha [ ( i - 1 ) * ( p + 1 ) + j ] + k [ i ] * alpha [ ( i - 1 ) * ( p + 1 ) + i - j ] ; } E [ i ] = ( 1 - k [ i ] * k [ i ] ) * E [ i - 1 ] ; } g = sqrt ( E [ p ] ) ; for ( int i = 0 ; i <= p ; i ++ ) { poles [ i ] = alpha [ p * ( p + 1 ) + i ] ; if ( ! do_gain_norm ) { if ( E [ p ] <= 0. ) { poles [ i ] = 0. ; } else { poles [ i ] /= g ; } } } FREE ( k ) ; FREE ( E ) ; FREE ( alpha ) ; }
levinson-durbin recursion outputting poles.
levinson - durbin recursion outputting poles .
// Copy final iteration coeffs to output array // nantest // g is nan now (negative) or 0 (if E[p] == 0.) -> complex // squareroot or division by zero. all boils down to all-zero // poles