decoder excitation now in 16-bit precision (was 32), which saves quite a bit
[speexdsp.git] / libspeex / lpc.c
index 0b8e2c0..f90254f 100644 (file)
    SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */
 
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include "lpc.h"
 
+#ifdef BFIN_ASM
+#include "lpc_bfin.h"
+#endif
 
 /* LPC analysis
  *
 /* Invented by N. Levinson in 1947, modified by J. Durbin in 1959.
  */
 
-#include "lpc.h"
-
-#ifdef FIXED_POINT
-#include <math.h>
-
 /* returns minimum mean square error    */
 spx_word32_t _spx_lpc(
 spx_coef_t       *lpc, /* out: [0...p-1] LPC coefficients      */
@@ -76,7 +79,6 @@ int          p
    int i, j;  
    spx_word16_t r;
    spx_word16_t error = ac[0];
-   spx_word16_t lpcq[10];
 
    if (ac[0] == 0)
    {
@@ -88,39 +90,43 @@ int          p
    for (i = 0; i < p; i++) {
 
       /* Sum up this iteration's reflection coefficient */
-      int rr = -ac[i + 1]<<13;
+      spx_word32_t rr = NEG32(SHL32(EXTEND32(ac[i + 1]),13));
       for (j = 0; j < i; j++) 
-         rr -= lpcq[j] * ac[i - j];
-      r = DIV32_16(rr,error+16);
-
+         rr = SUB32(rr,MULT16_16(lpc[j],ac[i - j]));
+#ifdef FIXED_POINT
+      r = DIV32_16(rr+PSHR32(error,1),ADD16(error,4));
+#else
+      r = rr/(error+.003*ac[0]);
+#endif
       /*  Update LPC coefficients and total error */
-      lpcq[i] = r;
+      lpc[i] = r;
       for (j = 0; j < i>>1; j++) 
       {
-         spx_word16_t tmp  = lpcq[j];
-         lpcq[j]     += MULT16_16_Q13(r,lpcq[i-1-j]);
-         lpcq[i-1-j] += MULT16_16_Q13(r,tmp);
+         spx_word16_t tmp  = lpc[j];
+         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) 
-         lpcq[j] += MULT16_16_Q13(lpcq[j],r);
+         lpc[j] = MAC16_16_P13(lpc[j],lpc[j],r);
 
-      error -= MULT16_16_Q13(r,MULT16_16_Q13(error,r));
+      error = SUB16(error,MULT16_16_Q13(r,MULT16_16_Q13(error,r)));
    }
-   for (i = 0; i < p; i++)
-      lpc[i] = lpcq[i];
    return error;
 }
 
 
+#ifdef FIXED_POINT
+
 /* Compute the autocorrelation
  *                      ,--,
  *              ac(i) = >  x(n) * x(n-i)  for all n
  *                      `--'
  * for lags between 0 and lag-1, and x == 0 outside 0...n-1
  */
-#include <stdio.h>
+
+#ifndef OVERRIDE_SPEEX_AUTOCORR
 void _spx_autocorr(
-const float *x,   /*  in: [0...n-1] samples x   */
+const spx_word16_t *x,   /*  in: [0...n-1] samples x   */
 spx_word16_t       *ac,  /* out: [0...lag-1] ac values */
 int          lag, 
 int          n
@@ -132,8 +138,8 @@ int          n
    int shift, ac_shift;
    
    for (j=0;j<n;j++)
-      ac0 += floor(x[j]*x[j])/256;
-   ac0 += n;
+      ac0 = ADD32(ac0,SHR32(MULT16_16(x[j],x[j]),8));
+   ac0 = ADD32(ac0,n);
    shift = 8;
    while (shift && ac0<0x40000000)
    {
@@ -153,61 +159,19 @@ int          n
       d=0;
       for (j=i;j<n;j++)
       {
-         d += ((int)(floor(x[j]) * floor(x[j-i]))) >> shift;
+         d = ADD32(d,SHR32(MULT16_16(x[j],x[j-i]), shift));
       }
       
-      ac[i] = d >> ac_shift;
+      ac[i] = SHR32(d, ac_shift);
    }
-   /*ac[0] = 8192.;*/
-   /*printf ("%d %d %d, %f\n", ac0, shift, ac_shift, ac[0]);*/
 }
+#endif
 
 
 #else
 
 
 
-/* returns minimum mean square error    */
-spx_word32_t _spx_lpc(
-spx_coef_t       *lpc, /* out: [0...p-1] LPC coefficients      */
-const float *ac,  /* in:  [0...p] autocorrelation values  */
-int          p
-)
-{
-   int i, j;  float r, error = ac[0];
-
-   if (ac[0] == 0)
-   {
-      for (i = 0; i < p; i++)
-         lpc[i] = 0;
-      return 0;
-   }
-
-   for (i = 0; i < p; i++) {
-
-      /* Sum up this iteration's reflection coefficient */
-      r = -ac[i + 1];
-      for (j = 0; j < i; j++) 
-         r -= lpc[j] * ac[i - j];
-      r /= error;
-
-      /*  Update LPC coefficients and total error */
-      lpc[i] = r;
-      for (j = 0; j < i/2; j++) 
-      {
-         float tmp  = lpc[j];
-         lpc[j]     += r * lpc[i-1-j];
-         lpc[i-1-j] += r * tmp;
-      }
-      if (i & 1) 
-         lpc[j] += lpc[j] * r;
-
-      error *= 1.0 - r * r;
-   }
-   return error;
-}
-
-
 /* Compute the autocorrelation
  *                      ,--,
  *              ac(i) = >  x(n) * x(n-i)  for all n
@@ -215,7 +179,7 @@ int          p
  * for lags between 0 and lag-1, and x == 0 outside 0...n-1
  */
 void _spx_autocorr(
-const float *x,   /*  in: [0...n-1] samples x   */
+const spx_word16_t *x,   /*  in: [0...n-1] samples x   */
 float       *ac,  /* out: [0...lag-1] ac values */
 int          lag, 
 int          n