Main Archive Specials Wiki | FAQ Links Submit Forum

 LPC analysis (autocorrelation + Levinson-Durbin recursion)References : Posted by mail[AT]mutagene[DOT]netNotes : 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; }

 CommentsAdded on : 18/02/04 by abhimj[ AT ]lycos[ DOT ]comComment : 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 ]deComment : 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 ]ukComment : 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