# levinson durbin algo (again)

Started by 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;
a      = 32768;
a_temp = 32768;
alpha     = r;

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;
>     a      = 32768;
>     a_temp = 32768;
>     alpha     = r;
>
>     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;
>> a = 32768;
>> a_temp = 32768;
>> alpha = r;
>>
>> 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;
>   a      = 32768;
>   a_temp = 32768;
>   alpha     = r;
>
>   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
>
```