.module/boot=3/boot=4        durbin_single; 

{ DURBIN.DSP - single precision Levinso-Durbin routine
  
  INPUT : 
    i4 -> buffer with autocorrelated speech (pm)
    l4 =  0
    i0 -> buffer for reflection coeffs
    l0 =  0

  OUTPUT  reflection coeffs calculated
    mr1 = Ep (minimum total squared prediction error

  FUNCTIONS CALLED:
    None

  DESCRIPTON:
    The routine implements Durbins recursion method of solving a
    set of linaer equations forming a Toeplitz matrix. The algorithm in 
    C is as follows:
      Where R[] is the autocorrelation, and k[] the reflection coeffs.
      e[] is the total squared error, and a[][] is the predictor coeff
      matrix (Since only the i'th and the i+1'th column is used at any
      one time, the matrix is implemented as two (a_old and a_new) columns 
      swapping place after each iteration.
     
     e[0] = R[0]     k[1] = R[1] / e[0]     alpha[1][1] = k[1]     e[1] = (1 - k[1]*k[1]) * e[0]   
     for (i=2; i<=N; i++)     begin
       k[i] = 0       for (j=1; j<=i-1; j++)         k[i] = k[i] + R[i-j] * alpha[i-1][j]       k[i] = R[i] - k[i]       k[i] = k[i] / e[i-1]       alpha[i][i] = k[i]       for (j=i-1; j>0; j++)         alpha[i][j] = alpha[i-1][j] - k[i]*alpha[i-1][i-j]       e[i] = (1 - k[i]*k[i]) * e[i-1]     end
}
{Include constant definitions}
#include "lpc.h"

.entry    levinson;

.external overflow;

.global e;

.var/dm/ram i_1;
.var/dm/ram e[N+1];             {error values} 
.var/dm/ram a_new[N],a_old[N];
.var/dm/ram ap_new,ap_old;    {pointers to a_*}
.var/dm/ram p2_k_i;           {pointer to k[i]}
.var/dm     p2_autocor_speech;

  {determines the format that a-values are stored in 
   format: (SBITS+1).(16-SBITS-1)}
.const SBITS  =  3;
.const NSBITS = -SBITS;

levinson:
  i1 = ^a_new; l1 = 0;
  dm(ap_new) = i1;
  i2 = ^a_old; l2 = 0; 
  dm(ap_old) = i2;
  dm(p2_autocor_speech) = i4;
  i5 = ^e;     l5 = 0;

  m2 = -1;
  m6 = -1;

  se = NSBITS;

  /* e[0] = R[0] */
  ax0 = pm(i4,m5);
  dm(i5,m5) = ax0;

  /* k[1] = R[1]/e[0] */
   {ax0 = e[0] = divisor}
  ay1 = pm(i4,m4); {MSW of dividend}
  ay0 = 0000;      {LSW of dividend}
  divide(ax0,ay1);
  ar = -ay0; {reverse sign of k before storing}
  dm(i0,m1) = ar;
  dm(p2_k_i) =  i0;

  /*  a_old[1] = k[1] */
  si = ay0;
  sr = ashift si (hi);  {store in (SBITS+1).(16-SBITS-1) format}
  dm(i2,m0) = sr1;

  /* e[1] = (1 - k[1]*k[1])*e[0] */
   {ay0 = k[1]}
  mx0 = ay0;
  my0 = ay0;
  mr0 = 0xffff;   {mr = 1 in 1.31 format}
  mr1 = 0x7fff;
  mr = mr - mx0*my0 (ss);
   {ax0 = e(0)}
  my0 = ax0;
  mr = mr1 * my0 (ss);
  dm(i5,m4) = mr1; 

  /* for(i = 2; i <= N; i++) */
  cntr = N-1;
  ax0 = 1;    {i-1}
  dm(i_1) = ax0;
  do pass_two until ce;
  
    /* k[i] = 0 */
    mr = 0;
    
    /* for(j = 1; j <= i-1; j++) */
    ay0 = dm(i_1);
    cntr = ay0;
    m3 = ay0;  {i-1}
    m7 = ay0;  {i-1}

    /* prepare: k[i] = k[i] + R[i-j]*a_old[j] */
    i2 = dm(ap_old);            l2 = 0;
    i4 = dm(p2_autocor_speech); l4 = 0; 
    modify(i4,m7);  {->R[i-1]}

    /* loop */
    do calc_ks until ce;
      mx0 = pm(i4,m6);
      my0 = dm(i2,m1);
    calc_ks: mr = mr + mx0*my0 (ss);
    if mv call overflow;

    /* k[i] = R[i] - k[i] */
    i4 = dm(p2_autocor_speech);  l4 = 0;
    modify(i4,m7);
    modify(i4,m5);    {->R[i]}
    si  = pm(i4,m4);  {R[i]}
    sr  = ashift si (hi);  {shift to (SBITS+1).(16-SBITS-1) format}
    ay1 = mr1;        {k[i]}
    ar  = sr1 - ay1;
    if av call overflow;

    /* k[i] = k[i]/e[i-1] */
    i5 = ^e;        l5 = 0;
    modify(i5,m7);
    ax0 = dm(i5,m5); {e[i-1]}
    ay1 = ar;      {MSW of k[i]}
    ay0 = 0000;    {LSW of k[i]}  
     {overflow check}
     si = ax0;
     sr = ashift si (hi);
     ar = sr1 - ay1; {e[i-1] - k[i]}
     if ge jump e_ok;
       {call overflow;}
       si = 0x7fff;   {sat k[i]}
       sr = ashift si (hi);
       {ay1 = sr1;}
     e_ok:
    divide(ax0,ay1);
    si = ay0;
    sr = ashift si by SBITS(hi);  {shift to 1.15 format before storing}
    i0 = dm(p2_k_i);  l0 = 0;
    ay1 = sr1;
    ar  = -ay1; {reverse sign of k before storing}
    dm(i0,m1) = ar;  {k[i] store}
    dm(p2_k_i) = i0;

    /* a_new[i] = k[i] */
    i1 = dm(ap_new);  l1 = 0;
    modify(i1,m3); {->a_new[i]}
    dm(i1,m2) = ay0; 

    /* for(j = i-1; j>0; j--) */
    cntr = dm(i_1);
    
    /*prepare: a_new[j] = a_old[j] - k[i]*a_old[i-j] */
    i2 = dm(ap_old);  l2 = 0;
    modify(i2,m3);   {modify by j (= i-1)}
    modify(i2,m2);   {-> a_old[j]}
    i0 = dm(ap_old); {-> a_old[i-j]} l0 = 0;
    mx0 = sr1;       {k[i]}

    /* loop */
    do calc_as until ce;
      mr0 = 0;
      mr1 = dm(i2,m2);  {a_old[j]}
      my0 = dm(i0,m1);  {a_old[i-j]}
      mr  = mr - mx0*my0 (ss);
      if mv {sat mr} call overflow;
    calc_as: dm(i1,m2) = mr1;

    /* e[i] = (1 - k[i]*k[i]) * e[i-1] */
     {ay0 = k[i]}
    mx0 = sr1;
    my0 = sr1;
    mr0 = 0xffff;  {mr = 1 in 1.31 format}
    mr1 = 0x7fff;
    mr = mr - mx0*my0 (ss);
    if mv call overflow;
     {ax0 = e(i-1)}
    my0 = ax0;
    mr = mr1 * my0 (ss);
    dm(i5,m4) = mr1;

    /* switch the a pointers */
    ax0 = dm(ap_old);
    ay0 = dm(ap_new);
    dm(ap_new) = ax0;
    dm(ap_old) = ay0;
    
    /* i++ */
    ay0 = dm(i_1);
    ar  = ay0 + 1;
pass_two: dm(i_1) = ar;
  
  rts;

.endmod;