Making new PLC code work in fixed-point even though it's still using float

arithmetic.

2 files changed, 53 insertions, 56 deletions

diff --git a/libcelt/celt.c b/libcelt/celt.c
index 642099f..41effc4 100644
--- a/libcelt/celt.c
+++ b/libcelt/celt.c
@@ -55,6 +55,9 @@
 #include <stdarg.h>
 
 #define LPC_ORDER 24
+#if !defined(FIXED_POINT) || defined(NEW_PLC)
+#include "plc.c"
+#endif
 
 static const celt_word16 preemph = QCONST16(0.8f,15);
 
@@ -1092,8 +1095,8 @@ struct CELTDecoder {
 
    celt_word16 *oldBandE;
    
-#ifndef FIXED_POINT
-   celt_word16 *lpc;
+#ifdef NEW_PLC
+   float *lpc;
 #endif
 
    int last_pitch_index;
@@ -1170,14 +1173,14 @@ CELTDecoder *celt_decoder_create(const CELTMode *mode, int channels, int *error)
    
    st->preemph_memD = (celt_sig*)celt_alloc(C*sizeof(celt_sig));
 
-#ifndef FIXED_POINT
-   st->lpc = (celt_word16*)celt_alloc(C*LPC_ORDER*sizeof(celt_word16));
+#ifdef NEW_PLC
+   st->lpc = (float*)celt_alloc(C*LPC_ORDER*sizeof(float));
 #endif
 
    st->loss_count = 0;
 
    if ((st->decode_mem!=NULL) && (st->out_mem!=NULL) && (st->oldBandE!=NULL) &&
-#ifndef FIXED_POINT
+#ifdef NEW_PLC
          (st->lpc!=NULL) &&
 #endif
        (st->preemph_memD!=NULL))
@@ -1223,7 +1226,7 @@ void celt_decoder_destroy(CELTDecoder *st)
    celt_free(st->oldBandE);
    celt_free(st->preemph_memD);
 
-#ifndef FIXED_POINT
+#ifdef NEW_PLC
    celt_free(st->lpc);
 #endif
    
@@ -1232,9 +1235,6 @@ void celt_decoder_destroy(CELTDecoder *st)
    celt_free(st);
 }
 
-#ifndef FIXED_POINT
-#include "plc.c"
-#endif
 static void celt_decode_lost(CELTDecoder * restrict st, celt_word16 * restrict pcm)
 {
    int c, N;
@@ -1270,7 +1270,7 @@ static void celt_decode_lost(CELTDecoder * restrict st, celt_word16 * restrict p
          fade = 0;
    }
 
-#ifdef FIXED_POINT
+#ifndef NEW_PLC
    offset = MAX_PERIOD-pitch_index;
    ALLOC(freq,C*N, celt_sig); /**< Interleaved signal MDCTs */
    while (offset+len >= MAX_PERIOD)
@@ -1285,12 +1285,12 @@ static void celt_decode_lost(CELTDecoder * restrict st, celt_word16 * restrict p
 #else
    for (c=0;c<C;c++)
    {
-      float e[MAX_PERIOD];
+      celt_word32 e[MAX_PERIOD];
       float exc[MAX_PERIOD];
       float ac[LPC_ORDER+1];
       float decay = 1;
       float S1=0;
-      celt_word32 mem[LPC_ORDER]={0};
+      float mem[LPC_ORDER]={0};
 
       offset = MAX_PERIOD-pitch_index;
       for (i=0;i<MAX_PERIOD;i++)
@@ -1341,7 +1341,7 @@ static void celt_decode_lost(CELTDecoder * restrict st, celt_word16 * restrict p
             decay *= decay;
          }
          e[i] = decay*exc[offset+i];
-         S1 += st->out_mem[offset+i]*st->out_mem[offset+i];
+         S1 += st->out_mem[offset+i]*1.*st->out_mem[offset+i];
       }
 
       iir(e, st->lpc, e, len+st->mode->overlap, LPC_ORDER, mem);
@@ -1349,7 +1349,7 @@ static void celt_decode_lost(CELTDecoder * restrict st, celt_word16 * restrict p
       {
          float ratio, S2=0;
          for (i=0;i<len+overlap;i++)
-            S2 += e[i]*e[i];
+            S2 += e[i]*1.*e[i];
          ratio = sqrt((S1+1)/(S2+1));
          if (ratio < 1)
             for (i=0;i<len+overlap;i++)
@@ -1363,20 +1363,20 @@ static void celt_decode_lost(CELTDecoder * restrict st, celt_word16 * restrict p
          previous and next frames */
       for (i=0;i<overlap/2;i++)
       {
-         float tmp1, tmp2;
-         tmp1 = e[i          ]*st->mode->window[i          ] -
-                e[overlap-i-1]*st->mode->window[overlap-i-1];
-         tmp2 = e[N+overlap-1-i]*st->mode->window[i] +
-                e[N+i          ]*st->mode->window[overlap-i-1];
-         tmp1 *= fade;
-         tmp2 *= fade;
-         st->out_mem[C*(MAX_PERIOD+i)+c] = tmp2*st->mode->window[overlap-i-1];
-         st->out_mem[C*(MAX_PERIOD+overlap-i-1)+c] = tmp2*st->mode->window[i];
-         st->out_mem[C*(MAX_PERIOD-N+i)+c] += tmp1*st->mode->window[i];
-         st->out_mem[C*(MAX_PERIOD-N+overlap-i-1)+c] -= tmp1*st->mode->window[overlap-i-1];
+         celt_word32 tmp1, tmp2;
+         tmp1 = MULT16_32_Q15(st->mode->window[i          ], e[i          ]) -
+                MULT16_32_Q15(st->mode->window[overlap-i-1], e[overlap-i-1]);
+         tmp2 = MULT16_32_Q15(st->mode->window[i],           e[N+overlap-1-i]) +
+                MULT16_32_Q15(st->mode->window[overlap-i-1], e[N+i          ]);
+         tmp1 = MULT16_32_Q15(fade, tmp1);
+         tmp2 = MULT16_32_Q15(fade, tmp2);
+         st->out_mem[C*(MAX_PERIOD+i)+c] = MULT16_32_Q15(st->mode->window[overlap-i-1], tmp2);
+         st->out_mem[C*(MAX_PERIOD+overlap-i-1)+c] = MULT16_32_Q15(st->mode->window[i], tmp2);
+         st->out_mem[C*(MAX_PERIOD-N+i)+c] += MULT16_32_Q15(st->mode->window[i], tmp1);
+         st->out_mem[C*(MAX_PERIOD-N+overlap-i-1)+c] -= MULT16_32_Q15(st->mode->window[overlap-i-1], tmp1);
       }
       for (i=0;i<N-overlap;i++)
-         st->out_mem[C*(MAX_PERIOD-N+overlap+i)+c] = fade*e[overlap+i];
+         st->out_mem[C*(MAX_PERIOD-N+overlap+i)+c] = MULT16_32_Q15(fade, e[overlap+i]);
    }
 #endif
 
diff --git a/libcelt/plc.c b/libcelt/plc.c
index 8f5100c..8eb12d5 100644
--- a/libcelt/plc.c
+++ b/libcelt/plc.c
@@ -1,16 +1,17 @@
 
+#ifndef NEW_PLC
+#define NEW_PLC
+#endif
 
-
-
-celt_word32 _celt_lpc(
-celt_word16       *lpc, /* out: [0...p-1] LPC coefficients      */
-const celt_word16 *ac,  /* in:  [0...p] autocorrelation values  */
+float _celt_lpc(
+      float       *lpc, /* out: [0...p-1] LPC coefficients      */
+const float *ac,  /* in:  [0...p] autocorrelation values  */
 int          p
 )
 {
    int i, j;  
-   celt_word16 r;
-   celt_word16 error = ac[0];
+   float r;
+   float error = ac[0];
 
    if (ac[0] == 0)
    {
@@ -22,38 +23,34 @@ int          p
    for (i = 0; i < p; i++) {
       
       /* Sum up this iteration's reflection coefficient */
-      celt_word32 rr = NEG32(SHL32(EXTEND32(ac[i + 1]),13));
+      float rr = -ac[i + 1];
       for (j = 0; j < i; j++) 
-         rr = SUB32(rr,MULT16_16(lpc[j],ac[i - j]));
-#ifdef FIXED_POINT
-      r = DIV32_16(rr+PSHR32(error,1),ADD16(error,1));
-#else
+         rr = rr - lpc[j]*ac[i - j];
       r = rr/(error+1e-15);
-#endif
       /*  Update LPC coefficients and total error */
       lpc[i] = r;
       for (j = 0; j < i>>1; j++) 
       {
-         celt_word16 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);
+         float tmp  = lpc[j];
+         lpc[j]     = lpc[j    ] + r*lpc[i-1-j];
+         lpc[i-1-j] = lpc[i-1-j] + r*tmp;
       }
       if (i & 1) 
-         lpc[j] = MAC16_16_P13(lpc[j],lpc[j],r);
+         lpc[j] = lpc[j] + lpc[j]*r;
       
-      error = SUB16(error,MULT16_16_Q13(r,MULT16_16_Q13(error,r)));
+      error = error - r*r*error;
       if (error<.00001*ac[0])
          break;
    }
    return error;
 }
 
-void fir(const celt_word16 *x,
-         const celt_word16 *num,
-         celt_word16 *y,
+void fir(const float *x,
+         const float *num,
+         float *y,
          int N,
          int ord,
-         celt_word32 *mem)
+         float *mem)
 {
    int i,j;
 
@@ -73,12 +70,12 @@ void fir(const celt_word16 *x,
    }
 }
 
-void iir(const celt_word16 *x,
-         const celt_word16 *den,
-         celt_word16 *y,
+void iir(const celt_word32 *x,
+         const float *den,
+         celt_word32 *y,
          int N,
          int ord,
-         celt_word32 *mem)
+         float *mem)
 {
    int i,j;
    for (i=0;i<N;i++)
@@ -98,9 +95,9 @@ void iir(const celt_word16 *x,
 }
 
 void _celt_autocorr(
-                   const celt_word16 *x,   /*  in: [0...n-1] samples x   */
+                   const float *x,   /*  in: [0...n-1] samples x   */
                    float       *ac,  /* out: [0...lag-1] ac values */
-                   const float       *window,
+                   const celt_word16       *window,
                    int          overlap,
                    int          lag, 
                    int          n
@@ -115,8 +112,8 @@ void _celt_autocorr(
       xx[i] = x[i];
    for (i=0;i<overlap;i++)
    {
-      xx[i] *= window[i];
-      xx[n-i-1] *= window[i];
+      xx[i] *= (1./Q15ONE)*window[i];
+      xx[n-i-1] *= (1./Q15ONE)*window[i];
    }
    while (lag>=0)
    {