decoder excitation now in 16-bit precision (was 32), which saves quite a bit
[speexdsp.git] / libspeex / lpc.c
index d53c438..f90254f 100644 (file)
 #include "config.h"
 #endif
 
+#include "lpc.h"
 
+#ifdef BFIN_ASM
+#include "lpc_bfin.h"
+#endif
 
 /* LPC analysis
  *
@@ -65,8 +69,6 @@
 /* Invented by N. Levinson in 1947, modified by J. Durbin in 1959.
  */
 
-#include "lpc.h"
-
 /* returns minimum mean square error    */
 spx_word32_t _spx_lpc(
 spx_coef_t       *lpc, /* out: [0...p-1] LPC coefficients      */
@@ -88,11 +90,11 @@ int          p
    for (i = 0; i < p; i++) {
 
       /* Sum up this iteration's reflection coefficient */
-      spx_word32_t rr = -SHL(ac[i + 1],13);
+      spx_word32_t rr = NEG32(SHL32(EXTEND32(ac[i + 1]),13));
       for (j = 0; j < i; j++) 
          rr = SUB32(rr,MULT16_16(lpc[j],ac[i - j]));
 #ifdef FIXED_POINT
-      r = DIV32_16(rr,ADD16(error,16));
+      r = DIV32_16(rr+PSHR32(error,1),ADD16(error,4));
 #else
       r = rr/(error+.003*ac[0]);
 #endif
@@ -101,11 +103,11 @@ int          p
       for (j = 0; j < i>>1; j++) 
       {
          spx_word16_t tmp  = lpc[j];
-         lpc[j]     = MAC16_16_Q13(lpc[j],r,lpc[i-1-j]);
-         lpc[i-1-j] = MAC16_16_Q13(lpc[i-1-j],r,tmp);
+         lpc[j]     = MAC16_16_P13(lpc[j],r,lpc[i-1-j]);
+         lpc[i-1-j] = MAC16_16_P13(lpc[i-1-j],r,tmp);
       }
       if (i & 1) 
-         lpc[j] = MAC16_16_Q13(lpc[j],lpc[j],r);
+         lpc[j] = MAC16_16_P13(lpc[j],lpc[j],r);
 
       error = SUB16(error,MULT16_16_Q13(r,MULT16_16_Q13(error,r)));
    }
@@ -122,6 +124,7 @@ int          p
  * for lags between 0 and lag-1, and x == 0 outside 0...n-1
  */
 
+#ifndef OVERRIDE_SPEEX_AUTOCORR
 void _spx_autocorr(
 const spx_word16_t *x,   /*  in: [0...n-1] samples x   */
 spx_word16_t       *ac,  /* out: [0...lag-1] ac values */
@@ -135,7 +138,7 @@ int          n
    int shift, ac_shift;
    
    for (j=0;j<n;j++)
-      ac0 = ADD32(ac0,SHR(MULT16_16(x[j],x[j]),8));
+      ac0 = ADD32(ac0,SHR32(MULT16_16(x[j],x[j]),8));
    ac0 = ADD32(ac0,n);
    shift = 8;
    while (shift && ac0<0x40000000)
@@ -156,12 +159,13 @@ int          n
       d=0;
       for (j=i;j<n;j++)
       {
-         d = ADD32(d,SHR(MULT16_16(x[j],x[j-i]), shift));
+         d = ADD32(d,SHR32(MULT16_16(x[j],x[j-i]), shift));
       }
       
-      ac[i] = SHR(d, ac_shift);
+      ac[i] = SHR32(d, ac_shift);
    }
 }
+#endif
 
 
 #else