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

{ DURBIN2.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 = Vp minimum total squared prediction error (normalized)

  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

  In this version the alpha's (a's) are stored as 32 bit numbers.  
}

#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[2*N],a_old[2*N]; {msw0, lsw0, msw1, lsw1,.....}
.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).(32-SBITS-1)}
.const SBITS  =  4;
.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).(32-SBITS-1) format}
  dm(i2,m1) = sr1;
  dm(i2,m1) = sr0;

  /* 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 */
    ay0 = 0;  {LSW}
    ay1 = 0;  {MSW}

    /* for(j = 1; j <= i-1; j++) */
    ax0 = dm(i_1);
    cntr = ax0;
    m3 = ax0;  {i-1}
    m7 = ax0;  {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;
      my1 = pm(i4,m6);         {R[i-j]}
      mx1 = dm(i2,m1);         {msw of a_old[j]}
      mx0 = dm(i2,m1);         {lsw of a_old[j]}
      mr  = mx0 * my1 (us);    {lsw * msw}
      mr0 = mr1;               {shift down 16 bits}
      mr1 = mr2;
      mr  = mr + mx1*my1 (ss); {msw * msw}
      if mv call overflow;
      ar  = mr0 + ay0;          {acum. lsw's}  
      ay0 = ar;
      ar  = mr1 + ay1 + c;      {acum. msw's}  
      if av call overflow;
    calc_ks: ay1 = ar;
      
    /* 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).(32-SBITS-1) format}
    {ay0 = LSW of k[i]}
    ar   = sr0 - ay0;
    si   = ar;         {store for double precision upshift}
    {ay1 = MSW of k[i]}
    ar   = sr1 - ay1 + c - 1;
    if av call overflow;
    
/*    sr  = lshift si by SBITS (lo);
      si  = ar;
      se  = exp si (hi);
      ay1 = se;
        se  = NSBITS;
      ax1 = SBITS;
      ar  = ax1 + ay1;
      if gt call overflow;
    sr  = sr or ashift si by SBITS (hi);
*/    
    /* k[i] = k[i]/e[i-1] */
    i5 = ^e;     l5 = 0;
    modify(i5,m7);   {->e[i-1]} 
    ax0 = dm(i5,m5); {e[i-1]}
    ay1 = ar {sr1};    {MSW of k[i]}
    {overflow check}
      ar  = abs ax0;
      ay0 = ar;
      ar  = pass ay1;
      ar  = abs ar;
      ar  = ar - ay0;  {abs(k[i]) - abs(e[i-1])}
      if gt call overflow;
    ay0 = si {sr0};    {LSW of k[i]}  
    divide(ax0,ay1);
    
    si = ay0;
    sr = ashift si by SBITS (hi);
    ay0 = sr1;

    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] */
    si = ay0;
    sr = ashift si (hi);  {store in (SBITS+1).(32-SBITS-1) format}
    i1 = dm(ap_new);   l1 = 0;
    modify(i1,m3); 
    modify(i1,m3);   {->a_new[i].msw}
    modify(i1,m1);   {->a_new[i].lsw}
    dm(i1,m2) = sr0; {store lsw}
    dm(i1,m2) = sr1; {store msw} 

    /* 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(i2,m3);    {-> a_old[j+1].msw}
    modify(i2,m2);    {-> a_old[j].lsw} 
    i0  = dm(ap_old); {-> a_old[i-j].msw}   l0 = 0;
    my1 = ay0;        {k[i]}

    /* loop */
    do calc_as until ce;
      ay0 = dm(i2,m2);    {a_old[j].lsw}
      ay1 = dm(i2,m2);    {a_old[j].msw}
      mx1 = dm(i0,m1);    {a_old[i-j].msw}
      mx0 = dm(i0,m1);    {a_old[i-j].lsw}
      mr  = mx0*my1 (us); {lsw * msw}
      mr0 = mr1;          {shift down by 16 bits}  
      mr1 = mr2;  
      mr  = mr + mx1*my1 (ss); {msw * msw}
      if mv call overflow;
      ar  = ay0 - mr0;         {acum. lsw's}  
      dm(i1,m2) = ar;
      ar  = ay1 - mr1 + c - 1; {acum. msw's}  
      if av call overflow;
    calc_as: dm(i1,m2) = ar;

    /* e[i] = (1 - k[i]*k[i]) * e[i-1] */
     {my1 = k[i]}
    mx0 = my1;
    mr0 = 0xffff;  {mr = 1 in 1.31 format}
    mr1 = 0x7fff;
    mr = mr - mx0*my1 (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;