mdct_backward now does the WOLA, so there isn't much left in compute_inv_mdcts

libcelt/celt.c

@@ -102,7 +102,7 @@ CELTEncoder EXPORT *celt_encoder_create(const CELTMode *mode)
    
    st->in_mem = celt_alloc(st->overlap*C*sizeof(celt_sig_t));
    st->mdct_overlap = celt_alloc(st->overlap*C*sizeof(celt_sig_t));
-   st->out_mem = celt_alloc(MAX_PERIOD*C*sizeof(celt_sig_t));
+   st->out_mem = celt_alloc((MAX_PERIOD+st->overlap)*C*sizeof(celt_sig_t));
 
    st->oldBandE = (celt_word16_t*)celt_alloc(C*mode->nbEBands*sizeof(celt_word16_t));
 
@@ -185,38 +185,40 @@ static void compute_mdcts(const CELTMode *mode, const celt_word16_t * restrict w
 static void compute_inv_mdcts(const CELTMode *mode, const celt_word16_t * restrict window, celt_sig_t *X, celt_sig_t * restrict out_mem, celt_sig_t * restrict mdct_overlap)
 {
    int c, N4;
-   VARDECL(celt_word32_t, x);
-   VARDECL(celt_word32_t, tmp);
    const int C = CHANNELS(mode);
    const mdct_lookup *lookup = MDCT(mode);
    const int N = FRAMESIZE(mode);
    const int overlap = OVERLAP(mode);
-   SAVE_STACK;
-   ALLOC(x, 2*N, celt_word32_t);
-   ALLOC(tmp, N, celt_word32_t);
    N4 = (N-overlap)>>1;
    for (c=0;c<C;c++)
    {
       int j;
       if (C==1) {
-         mdct_backward(lookup, X, x, window, overlap);
+         mdct_backward(lookup, X, out_mem+C*(MAX_PERIOD-N-N4), window, overlap);
       } else {
+         VARDECL(celt_word32_t, x);
+         VARDECL(celt_word32_t, tmp);
+         SAVE_STACK;
+         ALLOC(x, 2*N, celt_word32_t);
+         ALLOC(tmp, N, celt_word32_t);
          /* De-interleaving the sub-frames */
          for (j=0;j<N;j++)
             tmp[j] = X[C*j+c];
+         /* Prevents problems from the imdct doing the overlap-add */
+         CELT_MEMSET(x+N4, 0, overlap);
          mdct_backward(lookup, tmp, x, window, overlap);
-      }
          /* The first and last part would need to be set to zero if we actually
          wanted to use them. */
          for (j=0;j<overlap;j++)
-         out_mem[C*(MAX_PERIOD-N)+C*j+c] = mdct_overlap[C*j+c]+x[j+N4];
+            out_mem[C*(MAX_PERIOD-N)+C*j+c] += x[j+N4];
          for (j=0;j<overlap;j++)
-         mdct_overlap[C*(overlap-j-1)+c] = x[2*N-j-N4-1];
+            out_mem[C*(MAX_PERIOD)+C*(overlap-j-1)+c] = x[2*N-j-N4-1];
          for (j=0;j<2*N4;j++)
             out_mem[C*(MAX_PERIOD-N)+C*(j+overlap)+c] = x[j+N4+overlap];
-   }
          RESTORE_STACK;
       }
+   }
+}
 
 int EXPORT celt_encode(CELTEncoder * restrict st, celt_int16_t * restrict pcm, unsigned char *compressed, int nbCompressedBytes)
 {
@@ -358,7 +360,7 @@ int EXPORT celt_encode(CELTEncoder * restrict st, celt_int16_t * restrict pcm, u
    denormalise_bands(st->mode, X, freq, bandE);
 
 
-   CELT_MOVE(st->out_mem, st->out_mem+C*N, C*(MAX_PERIOD-N));
+   CELT_MOVE(st->out_mem, st->out_mem+C*N, C*(MAX_PERIOD+st->overlap-N));
 
    compute_inv_mdcts(st->mode, st->mode->window, freq, st->out_mem, st->mdct_overlap);
    /* De-emphasis and put everything back at the right place in the synthesis history */
@@ -462,7 +464,7 @@ CELTDecoder EXPORT *celt_decoder_create(const CELTMode *mode)
    st->overlap = mode->overlap;
 
    st->mdct_overlap = celt_alloc(st->overlap*C*sizeof(celt_sig_t));
-   st->out_mem = celt_alloc(MAX_PERIOD*C*sizeof(celt_sig_t));
+   st->out_mem = celt_alloc((MAX_PERIOD+st->overlap)*C*sizeof(celt_sig_t));
    
    st->oldBandE = (celt_word16_t*)celt_alloc(C*mode->nbEBands*sizeof(celt_word16_t));
 
@@ -614,7 +616,7 @@ int EXPORT celt_decode(CELTDecoder * restrict st, unsigned char *data, int len,
    denormalise_bands(st->mode, X, freq, bandE);
 
 
-   CELT_MOVE(st->out_mem, st->out_mem+C*N, C*(MAX_PERIOD-N));
+   CELT_MOVE(st->out_mem, st->out_mem+C*N, C*(MAX_PERIOD+st->overlap-N));
    /* Compute inverse MDCTs */
    compute_inv_mdcts(st->mode, st->mode->window, freq, st->out_mem, st->mdct_overlap);
 

libcelt/mdct.c

@@ -189,18 +189,20 @@ void mdct_backward(const mdct_lookup *l, kiss_fft_scalar *in, kiss_fft_scalar *
    int i;
    int N, N2, N4;
    VARDECL(kiss_fft_scalar, f);
+   VARDECL(kiss_fft_scalar, f2);
    SAVE_STACK;
    N = l->n;
    N2 = N>>1;
    N4 = N>>2;
    ALLOC(f, N2, kiss_fft_scalar);
+   ALLOC(f2, N2, kiss_fft_scalar);
    
    /* Pre-rotate */
    {
       /* Temp pointers to make it really clear to the compiler what we're doing */
       const kiss_fft_scalar * restrict xp1 = in;
       const kiss_fft_scalar * restrict xp2 = in+N2-1;
-      kiss_fft_scalar * restrict yp = out;
+      kiss_fft_scalar * restrict yp = f2;
       kiss_fft_scalar *t = &l->trig[0];
       for(i=0;i<N4;i++) 
       {
@@ -213,7 +215,7 @@ void mdct_backward(const mdct_lookup *l, kiss_fft_scalar *in, kiss_fft_scalar *
    }
 
    /* Inverse N/4 complex FFT. This one should *not* downscale even in fixed-point */
-   cpx32_ifft(l->kfft, out, f, N4);
+   cpx32_ifft(l->kfft, f2, f, N4);
    
    /* Post-rotate */
    {
@@ -235,7 +237,7 @@ void mdct_backward(const mdct_lookup *l, kiss_fft_scalar *in, kiss_fft_scalar *
    {
       const kiss_fft_scalar * restrict fp1 = f;
       const kiss_fft_scalar * restrict fp2 = f+N2-1;
-      kiss_fft_scalar * restrict yp = out+N4;
+      kiss_fft_scalar * restrict yp = f2;
       for(i = 0; i < N4; i++)
       {
          *yp++ =-*fp1*2;
@@ -247,27 +249,44 @@ void mdct_backward(const mdct_lookup *l, kiss_fft_scalar *in, kiss_fft_scalar *
 
    /* Mirror on both sides for TDAC */
    {
+      kiss_fft_scalar * restrict fp1 = f2+N4-1;
       kiss_fft_scalar * restrict xp1 = out+N2-1;
-      kiss_fft_scalar * restrict xp2 = out+N2;
       kiss_fft_scalar * restrict yp1 = out+N4-overlap/2;
+      const celt_word16_t * restrict wp1 = window;
+      const celt_word16_t * restrict wp2 = window+overlap-1;
+      for(i = 0; i< N4-overlap/2; i++)
+      {
+         *xp1 = *fp1;
+         xp1--;
+         fp1--;
+      }
+      for(; i < N4; i++)
+      {
+         kiss_fft_scalar x1;
+         x1 = *fp1--;
+         *yp1++ +=-MULT16_32_Q15(*wp1, x1);
+         *xp1-- += MULT16_32_Q15(*wp2, x1);
+         wp1++;
+         wp2--;
+      }
+   }
+   {
+      kiss_fft_scalar * restrict fp2 = f2+N4;
+      kiss_fft_scalar * restrict xp2 = out+N2;
       kiss_fft_scalar * restrict yp2 = out+N-1-(N4-overlap/2);
       const celt_word16_t * restrict wp1 = window;
       const celt_word16_t * restrict wp2 = window+overlap-1;
       for(i = 0; i< N4-overlap/2; i++)
       {
-         *xp1 = *xp1;
+         *xp2 = *fp2;
-         *xp2 = *xp2;
-         xp1--;
          xp2++;
+         fp2++;
       }
       for(; i < N4; i++)
       {
-         kiss_fft_scalar x1, x2;
+         kiss_fft_scalar x2;
-         x1 = *xp1;
+         x2 = *fp2++;
-         x2 = *xp2;
-         *yp1++ =-MULT16_32_Q15(*wp1, x1);
          *yp2--  = MULT16_32_Q15(*wp1, x2);
-         *xp1-- = MULT16_32_Q15(*wp2, x1);
          *xp2++  = MULT16_32_Q15(*wp2, x2);
          wp1++;
          wp2--;

libcelt/mdct.h

@@ -60,7 +60,7 @@ void mdct_clear(mdct_lookup *l);
 /** Compute a forward MDCT and scale by 2/N */
 void mdct_forward(const mdct_lookup *l, kiss_fft_scalar *in, kiss_fft_scalar *out, const celt_word16_t *window, int overlap);
 
-/** Compute a backward MDCT (no scaling) */
+/** Compute a backward MDCT (no scaling) and performs weighted overlap-add */
 void mdct_backward(const mdct_lookup *l, kiss_fft_scalar *in, kiss_fft_scalar *out, const celt_word16_t * restrict window, int overlap);
 
 #endif