Forums

levinson durbin algo (again)

Started by John McDermick January 20, 2011
Hi,


I keep on getting a divide-by-zero exception when I execute the
following code:

Can somebody please tell me where my fixed-point math is "broken" ???
Thank you!

void levinsonDurbin(SINT32 *a, SINT32 *k, SINT32 *r, SINT32 n)
{
   SINT32 a_temp[10+1], alpha, epsilon;
   int i, j;
   SINT32 tmpval;

   k[0]      = 0;
   a[0]      = 32768;
   a_temp[0] = 32768;
   alpha     = r[0];

   for(i = 1; i <= n; i++)
   {
      epsilon = 0;
      for(j = 0; j < i; j++)
      {
         epsilon += ((a[j]*r[i-j]) >> 0);
      }

      epsilon = epsilon >> 15;

      a[i] = ((-epsilon) << 10);
      a[i] = a[i] / alpha;
      a[i] = (a[i] << 5);
      k[i] = a[i];
      tmpval = 32768 - ((k[i]*k[i]) >> 15);
      alpha = ((alpha*tmpval) >> 15);


      /* update a[ ] array into temporary array */
      for(j=1; j<i; j++)
      {
         a_temp[j] = a[j] + ((k[i]*a[i-j]) >> 15);
      }

      /* update a[ ] array */
      for(j=1; j<i; j++)
      {
         a[j] = a_temp[j];
      }
   }
}
On 01/20/2011 03:55 PM, John McDermick wrote:
> Hi, > > > I keep on getting a divide-by-zero exception when I execute the > following code: > > Can somebody please tell me where my fixed-point math is "broken" ??? > Thank you! > > void levinsonDurbin(SINT32 *a, SINT32 *k, SINT32 *r, SINT32 n) > { > SINT32 a_temp[10+1], alpha, epsilon; > int i, j; > SINT32 tmpval; > > k[0] = 0; > a[0] = 32768; > a_temp[0] = 32768; > alpha = r[0]; > > for(i = 1; i<= n; i++) > { > epsilon = 0; > for(j = 0; j< i; j++) > { > epsilon += ((a[j]*r[i-j])>> 0); > } > > epsilon = epsilon>> 15; > > a[i] = ((-epsilon)<< 10); > a[i] = a[i] / alpha;
Probably here ^^^^^
> a[i] = (a[i]<< 5); > k[i] = a[i]; > tmpval = 32768 - ((k[i]*k[i])>> 15); > alpha = ((alpha*tmpval)>> 15);
Very possibly, alpha is going to zero here ^^^^
> > > /* update a[ ] array into temporary array */ > for(j=1; j<i; j++) > { > a_temp[j] = a[j] + ((k[i]*a[i-j])>> 15); > } > > /* update a[ ] array */ > for(j=1; j<i; j++) > { > a[j] = a_temp[j]; > } > } > }
-- Randy Yates % "My Shangri-la has gone away, fading like Digital Signal Labs % the Beatles on 'Hey Jude'" yates@digitalsignallabs.com % http://www.digitalsignallabs.com % 'Shangri-La', *A New World Record*, ELO
On 01/20/2011 04:09 PM, Randy Yates wrote:
> On 01/20/2011 03:55 PM, John McDermick wrote: >> Hi, >> >> >> I keep on getting a divide-by-zero exception when I execute the >> following code: >> >> Can somebody please tell me where my fixed-point math is "broken" ??? >> Thank you! >> >> void levinsonDurbin(SINT32 *a, SINT32 *k, SINT32 *r, SINT32 n) >> { >> SINT32 a_temp[10+1], alpha, epsilon; >> int i, j; >> SINT32 tmpval; >> >> k[0] = 0; >> a[0] = 32768; >> a_temp[0] = 32768; >> alpha = r[0]; >> >> for(i = 1; i<= n; i++) >> { >> epsilon = 0; >> for(j = 0; j< i; j++) >> { >> epsilon += ((a[j]*r[i-j])>> 0); >> } >> >> epsilon = epsilon>> 15; >> >> a[i] = ((-epsilon)<< 10); >> a[i] = a[i] / alpha; > > Probably here ^^^^^ > >> a[i] = (a[i]<< 5); >> k[i] = a[i]; >> tmpval = 32768 - ((k[i]*k[i])>> 15); >> alpha = ((alpha*tmpval)>> 15); > > Very possibly, alpha is going to zero here ^^^^ > >> >> >> /* update a[ ] array into temporary array */ >> for(j=1; j<i; j++) >> { >> a_temp[j] = a[j] + ((k[i]*a[i-j])>> 15); >> } >> >> /* update a[ ] array */ >> for(j=1; j<i; j++) >> { >> a[j] = a_temp[j]; >> } >> } >> } > >
John, If you were simulating this under cygwin or linux, gdb would be your friend. In fact, the running gdb under emacs would be your bestest buddy! -- Randy Yates % "My Shangri-la has gone away, fading like Digital Signal Labs % the Beatles on 'Hey Jude'" yates@digitalsignallabs.com % http://www.digitalsignallabs.com % 'Shangri-La', *A New World Record*, ELO
On Jan 20, 9:55=A0pm, John McDermick <johnthedsp...@gmail.com> wrote:
> Hi, > > I keep on getting a divide-by-zero exception when I execute the > following code: > > Can somebody please tell me where my fixed-point math is "broken" ???
Don't know about the fixed-point math, but the *algorithm* contains a step where a divide-by-0 type of error might occur. You need to identify that step and test if the result is 0. If it is, you are dealing with a perfect AR(p) model with no noise, and can terminate the function. Rune
>Hi, > > >I keep on getting a divide-by-zero exception when I execute the >following code: > >Can somebody please tell me where my fixed-point math is "broken" ??? >Thank you! > >void levinsonDurbin(SINT32 *a, SINT32 *k, SINT32 *r, SINT32 n) >{ > SINT32 a_temp[10+1], alpha, epsilon; > int i, j; > SINT32 tmpval; > > k[0] = 0; > a[0] = 32768; > a_temp[0] = 32768; > alpha = r[0]; > > for(i = 1; i <= n; i++) > { > epsilon = 0; > for(j = 0; j < i; j++) > { > epsilon += ((a[j]*r[i-j]) >> 0); > } > > epsilon = epsilon >> 15; > > a[i] = ((-epsilon) << 10); > a[i] = a[i] / alpha; > a[i] = (a[i] << 5); > k[i] = a[i]; > tmpval = 32768 - ((k[i]*k[i]) >> 15); > alpha = ((alpha*tmpval) >> 15); > > > /* update a[ ] array into temporary array */ > for(j=1; j<i; j++) > { > a_temp[j] = a[j] + ((k[i]*a[i-j]) >> 15); > } > > /* update a[ ] array */ > for(j=1; j<i; j++) > { > a[j] = a_temp[j]; > } > } >}
what is the application of lavinson and durbin,, i mean where do u imolement
>