DSPRelated.com
Forums

Fixed-point implementation of levinson durbin algorithm

Started by John McDermick January 3, 2011
Hi,

Any links to a fixed point implementation of the Levinson Durbin
algorithm ...something like this:

void levdur(pLpcCoeff, pAutoCorr, pReflec, nOrder)

where pLpcCoeff is a pointer to the lpc coefficients
where pAutoCorr is a pointer to the autocorrelation coefficients
(range -40 to 40)
where pReflec is a pointer to the reflection coefficients
where nOrder is 10

The autocorrelation sequence is a 32bit array where values are in Q15
format.

Thanks

>Hi, > >Any links to a fixed point implementation of the Levinson Durbin >algorithm ...something like this: > >void levdur(pLpcCoeff, pAutoCorr, pReflec, nOrder) > >where pLpcCoeff is a pointer to the lpc coefficients >where pAutoCorr is a pointer to the autocorrelation coefficients >(range -40 to 40) >where pReflec is a pointer to the reflection coefficients >where nOrder is 10 > >The autocorrelation sequence is a 32bit array where values are in Q15 >format. > >Thanks > >
You may want to take a look at the fixed point code for the ITU-T G.729A or G.723.1 speech codecs (go to http://www.itu.int/rec/T-REC-G/en). Both have levinson algorithms of 10th order if I recall correctly. -Shawn
On Jan 3, 6:35=A0pm, "Shawn Stevenson" <s_stevenson@n_o_s_p_a_m.shaw.ca>
wrote:
> >Hi, > > >Any links to a fixed point implementation of the Levinson Durbin > >algorithm ...something like this: > > >void levdur(pLpcCoeff, pAutoCorr, pReflec, nOrder) > > >where pLpcCoeff is a pointer to the lpc coefficients > >where pAutoCorr is a pointer to the autocorrelation coefficients > >(range -40 to 40) > >where pReflec is a pointer to the reflection coefficients > >where nOrder is 10 > > >The autocorrelation sequence is a 32bit array where values are in Q15 > >format. > > >Thanks > > You may want to take a look at the fixed point code for the ITU-T G.729A =
or
> G.723.1 speech codecs (go tohttp://www.itu.int/rec/T-REC-G/en). Both have > levinson algorithms of 10th order if I recall correctly. > > -Shawn- Hide quoted text - > > - Show quoted text -
Thank you, I had a look at it already ...but I don't like that implementation because the autocorrelation values are split up in MSBs and LSBs (as far as I can see)....
On Jan 3, 4:53=A0pm, John McDermick <johnthedsp...@gmail.com> wrote:
> Hi, > > Any links to a fixed point implementation of the Levinson Durbin > algorithm ...something like this: > > void levdur(pLpcCoeff, pAutoCorr, pReflec, nOrder) > > where pLpcCoeff is a pointer to the lpc coefficients > where pAutoCorr is a pointer to the autocorrelation coefficients > (range -40 to 40) > where pReflec is a pointer to the reflection coefficients > where nOrder is 10 > > The autocorrelation sequence is a 32bit array where values are in Q15 > format.
now this is floating point. but it's pretty simple. do you know what to do to make it fixed? (addition/subtraction stays the same. multiplication/division will need some arithmetic shifting after/ before. pretty straight-forward.) /********************************************************************* * * * Levinson-Durbin algorithm for LPC coefficients * * * * input: * * n =3D LPC order * * r -> r[0] to r[n] autocorrelation values * * * * output: * * a -> a[0] to a[n] LPC coefficients, a[0] =3D 1.0 * * k -> k[0] to k[n] reflection coefficients, k[0] =3D unused * * * * * *********************************************************************/ #define N 32 /* greater than any conceivable LPC order */ levison_durbin(int n, double* r, double* a, double* k) { double a_temp[N+1], alpha, epsilon; /* n <=3D N =3D constant */ int i, j; k[0] =3D 0.0; /* unused */ a[0] =3D 1.0; a_temp[0] =3D 1.0; /* unnecessary but consistent */ alpha =3D r[0]; for(i=3D1; i<=3Dn; i++) { epsilon =3D 0.0; for(j=3D0; j<i; j++) { epsilon +=3D a[j]*r[i-j]; } a[i] =3D k[i] =3D -epsilon/alpha; alpha =3D alpha*(1.0 - k[i]*k[i]); /* update a[ ] array into temporary array */ for(j=3D1; j<i; j++) { a_temp[j] =3D a[j] + k[i]*a[i-j]; } /* update a[ ] array */ for(j=3D1; j<i; j++) { a[j] =3D a_temp[j]; } } }
> Thanks
FWIW. r b-j
On Jan 3, 9:04=A0pm, robert bristow-johnson <r...@audioimagination.com>
wrote:
> On Jan 3, 4:53=A0pm, John McDermick <johnthedsp...@gmail.com> wrote: > > > Hi, > > > Any links to a fixed point implementation of the Levinson Durbin > > algorithm ...something like this: > > > void levdur(pLpcCoeff, pAutoCorr, pReflec, nOrder) > > > where pLpcCoeff is a pointer to the lpc coefficients > > where pAutoCorr is a pointer to the autocorrelation coefficients > > (range -40 to 40) > > where pReflec is a pointer to the reflection coefficients > > where nOrder is 10 > > > The autocorrelation sequence is a 32bit array where values are in Q15 > > format. > > now this is floating point. =A0but it's pretty simple. =A0do you know wha=
t
> to do to make it fixed? =A0(addition/subtraction stays the same. > multiplication/division will need some arithmetic shifting after/ > before. =A0pretty straight-forward.) > > /********************************************************************* > * =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0=
=A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0*
> * =A0 Levinson-Durbin algorithm for LPC coefficients =A0 =A0 =A0 =A0 =A0 =
=A0 =A0 =A0 =A0 *
> * =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0=
=A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0*
> * =A0 input: =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =
=A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 *
> * =A0 n =A0=3D LPC order =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =
=A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 *
> * =A0 r -> r[0] to r[n] =A0 =A0autocorrelation values =A0 =A0 =A0 =A0 =A0=
=A0 =A0 =A0 =A0 =A0 =A0*
> * =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0=
=A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0*
> * =A0 output: =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0=
=A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0*
> * =A0 a -> a[0] to a[n] =A0 =A0LPC coefficients, a[0] =3D 1.0 =A0 =A0 =A0=
=A0 =A0 =A0 =A0 =A0*
> * =A0 k -> k[0] to k[n] =A0 =A0reflection coefficients, k[0] =3D unused =
=A0 =A0 =A0*
> * =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0=
=A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0*
> * =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0=
=A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0*
> *********************************************************************/ > > #define N =A0 =A0 =A0 32 =A0 =A0 =A0 =A0/* greater than any conceivable L=
PC order */
> > levison_durbin(int n, double* r, double* a, double* k) > =A0 =A0 { > =A0 =A0 double a_temp[N+1], alpha, epsilon; =A0 =A0 =A0 /* n <=3D N =3D c=
onstant =A0*/
> =A0 =A0 int i, j; > > =A0 =A0 k[0] =A0 =A0 =A0=3D 0.0; =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =
=A0 =A0 =A0/* unused */
> > =A0 =A0 a[0] =A0 =A0 =A0=3D 1.0; > =A0 =A0 a_temp[0] =3D 1.0; =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0/* unnecess=
ary but consistent */
> > =A0 =A0 alpha =A0 =A0 =3D r[0]; > > =A0 =A0 for(i=3D1; i<=3Dn; i++) > =A0 =A0 =A0 =A0 { > =A0 =A0 =A0 =A0 epsilon =3D 0.0; > =A0 =A0 =A0 =A0 for(j=3D0; j<i; j++) > =A0 =A0 =A0 =A0 =A0 =A0 { > =A0 =A0 =A0 =A0 =A0 =A0 epsilon +=3D a[j]*r[i-j]; > =A0 =A0 =A0 =A0 =A0 =A0 } > > =A0 =A0 =A0 =A0 a[i] =3D k[i] =3D -epsilon/alpha; > =A0 =A0 =A0 =A0 alpha =3D alpha*(1.0 - k[i]*k[i]); > > =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0/* update a[ ] array i=
nto temporary array */
> =A0 =A0 =A0 =A0 for(j=3D1; j<i; j++) > =A0 =A0 =A0 =A0 =A0 =A0 { > =A0 =A0 =A0 =A0 =A0 =A0 a_temp[j] =3D a[j] + k[i]*a[i-j]; > =A0 =A0 =A0 =A0 =A0 =A0 } > =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0/* update a[ ] array *=
/
> =A0 =A0 =A0 =A0 for(j=3D1; j<i; j++) > =A0 =A0 =A0 =A0 =A0 =A0 { > =A0 =A0 =A0 =A0 =A0 =A0 a[j] =3D a_temp[j]; > =A0 =A0 =A0 =A0 =A0 =A0 } > =A0 =A0 =A0 =A0 } > =A0 =A0 } > > > Thanks > > FWIW. > > r b-j
Thank you. I made a similar implementation and it worked :o) The only problem is that it crashes for autocorrelation sequences containing the same values (alpha is zero, division by zero!). However, my guess is that there is no signal x of length greater than 1 which has a corresponding autocorrelation sequence containing identical values.
 Can you tell me the whole code of vowel recognition ?
its urgent
>On Jan 3, 6:35=A0pm, "Shawn Stevenson" <s_stevenson@n_o_s_p_a_m.shaw.ca> >wrote: >> >Hi, >> >> >Any links to a fixed point implementation of the Levinson Durbin >> >algorithm ...something like this: >> >> >void levdur(pLpcCoeff, pAutoCorr, pReflec, nOrder) >> >> >where pLpcCoeff is a pointer to the lpc coefficients >> >where pAutoCorr is a pointer to the autocorrelation coefficients >> >(range -40 to 40) >> >where pReflec is a pointer to the reflection coefficients >> >where nOrder is 10 >> >> >The autocorrelation sequence is a 32bit array where values are in Q15 >> >format. >> >> >Thanks >> >> You may want to take a look at the fixed point code for the ITU-T G.729A
=
>or >> G.723.1 speech codecs (go tohttp://www.itu.int/rec/T-REC-G/en). Both
have
>> levinson algorithms of 10th order if I recall correctly. >> >> -Shawn- Hide quoted text - >> >> - Show quoted text - > >Thank you, I had a look at it already ...but I don't like that >implementation because the autocorrelation values are split up in MSBs >and LSBs (as far as I can see)....
There is floating point reference code for those codecs which doesn't do this kind of thing. Most of the codec fixed point reference code around is oriented towards 16 bit DSPs, so it usually has the kind of messy LSB/MSB stuff you've seen, as soon as the dynamic range of a section of code starts to grow. Steve _____________________________ Posted through www.DSPRelated.com
On 9/7/13 10:29 PM, Shacz wrote:
> Can you tell me the whole code of vowel recognition ? > its urgent
rots o' ruk, dude. if it were me, i would not be using Levinson-Durbin (which is a method for doing LPC) for the vowel/fricative decision, but would base that decision on how periodic the signal is. if you're trying to recognize *which* vowel it is, then LPC is useful but you better check out the research lit on it. LPC will tell you where the resonances (formants) are and someone, somewhere, has established a sorta dictionary that relates the formant space to the various shapes of the mouth and lips (and nasal coupling) associated with different vowel, and i am sure that is not independent on cultural/regional/language/accent issues (like the difference between "aye" and "ahh" for the same intended vowel). -- r b-j rbj@audioimagination.com "Imagination is more important than knowledge."
On Monday, January 3, 2011 4:53:28 PM UTC-5, John McDermick wrote:
> Hi, > > Any links to a fixed point implementation of the Levinson Durbin > algorithm ...something like this: > > void levdur(pLpcCoeff, pAutoCorr, pReflec, nOrder) > > where pLpcCoeff is a pointer to the lpc coefficients > where pAutoCorr is a pointer to the autocorrelation coefficients > (range -40 to 40) > where pReflec is a pointer to the reflection coefficients > where nOrder is 10 > > The autocorrelation sequence is a 32bit array where values are in Q15 > format. > > Thanks
Hello John, Finding the LPCs directly in fixed point is somewhat tricky since they often have a large dynamic range; however, an end run around the problem can be had by finding the PARCORS via the Leroux-Gueguen algorithm. This works very well in fixed point. This is detailed in Panos E. Papamichalis' book Practical Approaches to Speech Coding. IHTH, Clay
clay@claysturner.com wrote:

(snip)
> Finding the LPCs directly in fixed point is somewhat tricky > since they often have a large dynamic range; however, an end > run around the problem can be had by finding the PARCORS via > the Leroux-Gueguen algorithm. This works very well in fixed point.
How large is large? Floating point is designed for problems having a dynamic range of hundred of orders of magnitude, but it a reasonable relative error. The mass of a proton and the mass of the earth differ by over 50 orders of magnitude, the mass of the earth is known to about six digits and of the proton to about nine digits. Note that the relative uncertainties are about the same, where the absolute uncertainties are very different. Seems to me that in DSP 24 bits (about 7 decimal digits) is a large dynamic range. -- glen