Main Archive Specials Wiki | FAQ Links Submit Forum


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 {
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 {
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


Added on : 18/02/04 by 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?



Added on : 31/03/05 by 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.




Added on : 21/12/07 by 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?




Add your own comment
Comments are displayed in fixed width, no HTML code allowed!
Email:

Comment:

Are you human?



Site created and maintained by Bram
Graphic design by line.out | Server sponsered by fxpansion