### LPC analysis (autocorrelation + Levinson-Durbin recursion)

References : Posted by mail[AT]mutagene[DOT]net
Notes :
The autocorrelation function implements a warped autocorrelation, so that frequency resolution can be specified by the variable 'lambda'. Levinson-Durbin recursion calculates autoregression coefficients a and reflection coefficients (for lattice filter implementation) K. Comments for Levinson-Durbin function implement matlab version of the same function.

No optimizations.
Code :
//find the order-P autocorrelation array, R, for the sequence x of length L and warping of lambda
//wAutocorrelate(&pfSrc[stIndex],siglen,R,P,0);
wAutocorrelate(float * x, unsigned int L, float * R, unsigned int P, float lambda)
{
double    * dl = new double [L];
double    * Rt = new double [L];
double r1,r2,r1t;
R[0]=0;
Rt[0]=0;
r1=0;
r2=0;
r1t=0;
for(unsigned int k=0; k<L;k++)
{
Rt[0]+=double(x[k])*double(x[k]);

dl[k]=r1-double(lambda)*double(x[k]-r2);
r1 = x[k];
r2 = dl[k];
}
for(unsigned int i=1; i<=P; i++)
{
Rt[i]=0;
r1=0;
r2=0;
for(unsigned int k=0; k<L;k++)
{
Rt[i]+=double(dl[k])*double(x[k]);

r1t = dl[k];
dl[k]=r1-double(lambda)*double(r1t-r2);
r1 = r1t;
r2 = dl[k];
}
}
for(i=0; i<=P; i++)
R[i]=float(Rt[i]);
delete[] dl;
delete[] Rt;
}

// Calculate the Levinson-Durbin recursion for the autocorrelation sequence R of length P+1 and return the autocorrelation coefficients a and reflection coefficients K
LevinsonRecursion(unsigned int P, float *R, float *A, float *K)
{
double Am1[62];

if(R[0]==0.0) {
for(unsigned int i=1; i<=P; i++)
{
K[i]=0.0;
A[i]=0.0;
}}
else {
double km,Em1,Em;
unsigned int k,s,m;
for (k=0;k<=P;k++){
A[0]=0;
Am1[0]=0; }
A[0]=1;
Am1[0]=1;
km=0;
Em1=R[0];
for (m=1;m<=P;m++)                    //m=2:N+1
{
double err=0.0f;                    //err = 0;
for (k=1;k<=m-1;k++)            //for k=2:m-1
err += Am1[k]*R[m-k];        // err = err + am1(k)*R(m-k+1);
km = (R[m]-err)/Em1;            //km=(R(m)-err)/Em1;
K[m-1] = -float(km);
A[m]=(float)km;                        //am(m)=km;
for (k=1;k<=m-1;k++)            //for k=2:m-1
A[k]=float(Am1[k]-km*Am1[m-k]);  // am(k)=am1(k)-km*am1(m-k+1);
Em=(1-km*km)*Em1;                //Em=(1-km*km)*Em1;
for(s=0;s<=P;s++)                //for s=1:N+1
Am1[s] = A[s];                // am1(s) = am(s)
Em1 = Em;                        //Em1 = Em;
}
}
return 0;
}

Comments
from : abhimj[AT]lycos[DOT]com
comment : Hi, I the LPC coeffs are to be converted to cepstral coeffs, the formula for cepstral coeffs requires some 'gain term' that is calculated in the LPC analysis phase. Can you please tell how to get this gain term?

from : Christian[AT]savioursofsoul[DOT]de
comment : Blind Object Pascal Translation: -------------------------------- unit Levinson; interface type TDoubleArray = array of Double; TSingleArray = array of Single; implementation //find the P-order autocorrelation array, R, for the sequence x of length L and warping of lambda procedure Autocorrelate(x,R : TSingleArray; P : Integer; lambda : Single; l: Integer = -1); var dl,Rt : TDoubleArray; r1,r2,r1t : Double; k,i : Integer; begin // Initialization if l=-1 then l:=Length(x); SetLength(dl,l); SetLength(Rt,l); R[0]:=0; Rt[0]:=0; r1:=0; r2:=0; r1t:=0; for k:=0 to l-1 do begin Rt[0]:=Rt[0]+x[k]*x[k]; dl[k]:=r1-lambda*(x[k]-r2); r1:= x[k]; r2:= dl[k]; end; for i:=1 to P do begin Rt[i]:=0; r1:=0; r2:=0; for k:=0 to L-1 do begin Rt[i]:=Rt[i]+dl[k]*x[k]; r1t:= dl[k]; dl[k]:=r1-lambda*(r1t-r2); r1:=r1t; r2:=dl[k]; end; end; for i:=1 to P do R[i]:=Rt[i]; setlength(Rt,0); setlength(dl,0); end; // Calculate the Levinson-Durbin recursion for the autocorrelation sequence // R of length P+1 and return the autocorrelation coefficients a and reflection coefficients K procedure LevinsonRecursion(P : Integer; R,A,K : TSingleArray); var Am1 : TDoubleArray; i,j,s,m : Integer; km,Em1,Em : Double; err : Double; begin SetLength(Am1,62); if (R[0]=0.0) then begin for i:=1 to P do begin K[i]:=0.0; A[i]:=0.0; end; end else begin for j:=0 to P do begin A[0]:=0; Am1[0]:=0; end; A[0]:=1; Am1[0]:=1; km:=0; Em1:=R[0]; for m:=1 to P do begin err:=0.0; for j:=1 to m-1 do err:=err+Am1[j]*R[m-j]; km:=(R[m]-err)/Em1; K[m-1]:=-km; A[m]:=km; for j:=1 to m-1 do A[j]:=Am1[j]-km*Am1[m-j]; Em:=(1-km*km)*Em1; for s:=0 to P do Am1[s]:=A[s]; Em1:=Em; end; end; end; end.

from : jamie[AT]postlude[DOT]co[DOT]uk
comment : Hi, This loop: for (k=0;k<=P;k++){ A[0]=0; Am1[0]=0; } Looks like it shouldn't be there?