ftp.nice.ch/pub/next/unix/audio/cmusic.bs.N.tar.gz#/src/pv/lpa.c

This is lpa.c in view mode; [Download] [Up]

#include "pv.h"

/*
 * given N samples of digital waveform x, lpa computes M+1 coefficients
 * by maximum entropy (autocorrelation) method for spectral estimation--
 * these are returned in b[] (b[0] is always set to 1). lpa itself returns
 * the a0 (residual energy) coefficient.
 */
float lpa( x, N, b, M ) float x[], b[] ; int N, M ; {
 int i, j ;
 float s, at ;
 float a0, *rx, *rc ;
 char *malloc() ;
    rx = (float *) malloc( (M+2)*sizeof(float) ) ;
    rc = (float *) malloc( (M+2)*sizeof(float) ) ;
    for ( i = 0 ; i <= M + 1 ; i++ )
	for ( rx[i] = j = 0 ; j < N - i ; j++ )
	    rx[i] += x[j]*x[i + j] ;
    b[0] = 1. ;
    b[1] = rc[0] = rx[0] != 0. ? -rx[1]/rx[0] : 0. ;
    a0 = rx[0] + rx[1]*rc[0] ;
    for ( i = 1 ; i < M ; i++ ) {
	for ( s = j = 0 ; j <= i ; j++ )
	    s += rx[i - j + 1]*b[j] ;
	rc[i] = a0 != 0. ? -s/a0 : 0. ;
	for ( j = 1 ; j <= (i + 1)>>1 ; j++ ) {
	    at = b[j] + rc[i]*b[i-j+1] ;
	    b[i-j+1] += rc[i]*b[j] ;
	    b[j] = at ;
	}
	b[i+1] = rc[i] ;
	a0 += rc[i]*s ;
/*
	if ( a0 <= 0. ) {
	    fprintf( stderr, "lpa: matrix singularity" ) ;
	    return( a0 ) ;
	}
*/
    }
    free( rx ) ;
    free( rc ) ;
    return( a0 ) ;
}

These are the contents of the former NiCE NeXT User Group NeXTSTEP/OpenStep software archive, currently hosted by Netfuture.ch.