No longer storing the zero parts of the in[] array in the encoder
[opus.git] / libcelt / celt.c
1 /* (C) 2007 Jean-Marc Valin, CSIRO
2 */
3 /*
4    Redistribution and use in source and binary forms, with or without
5    modification, are permitted provided that the following conditions
6    are met:
7    
8    - Redistributions of source code must retain the above copyright
9    notice, this list of conditions and the following disclaimer.
10    
11    - Redistributions in binary form must reproduce the above copyright
12    notice, this list of conditions and the following disclaimer in the
13    documentation and/or other materials provided with the distribution.
14    
15    - Neither the name of the Xiph.org Foundation nor the names of its
16    contributors may be used to endorse or promote products derived from
17    this software without specific prior written permission.
18    
19    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
20    ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
21    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
22    A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE FOUNDATION OR
23    CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
24    EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
25    PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
26    PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
27    LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
28    NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
29    SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
30 */
31
32 #ifdef HAVE_CONFIG_H
33 #include "config.h"
34 #endif
35
36 #define CELT_C
37
38 #include "os_support.h"
39 #include "mdct.h"
40 #include <math.h>
41 #include "celt.h"
42 #include "pitch.h"
43 #include "kiss_fftr.h"
44 #include "bands.h"
45 #include "modes.h"
46 #include "entcode.h"
47 #include "quant_pitch.h"
48 #include "quant_bands.h"
49 #include "psy.h"
50 #include "rate.h"
51
52 #define MAX_PERIOD 1024
53
54 /** Encoder state 
55  @brief Encoder state
56  */
57 struct CELTEncoder {
58    const CELTMode *mode;     /**< Mode used by the encoder */
59    int frame_size;
60    int block_size;
61    int nb_blocks;
62    int overlap;
63    int channels;
64    
65    ec_byte_buffer buf;
66    ec_enc         enc;
67
68    celt_word16_t preemph;
69    celt_sig_t *preemph_memE;
70    celt_sig_t *preemph_memD;
71
72    kiss_fftr_cfg fft;
73    struct PsyDecay psy;
74
75    celt_sig_t *in_mem;
76    celt_sig_t *mdct_overlap;
77    celt_sig_t *out_mem;
78
79    celt_word16_t *oldBandE;
80 };
81
82
83
84 CELTEncoder *celt_encoder_create(const CELTMode *mode)
85 {
86    int N, B, C;
87    CELTEncoder *st;
88
89    if (check_mode(mode) != CELT_OK)
90       return NULL;
91
92    N = mode->mdctSize;
93    B = mode->nbMdctBlocks;
94    C = mode->nbChannels;
95    st = celt_alloc(sizeof(CELTEncoder));
96    
97    st->mode = mode;
98    st->frame_size = B*N;
99    st->block_size = N;
100    st->nb_blocks  = B;
101    st->overlap = mode->overlap;
102
103    ec_byte_writeinit(&st->buf);
104    ec_enc_init(&st->enc,&st->buf);
105
106    st->fft = kiss_fftr_alloc(MAX_PERIOD, 0, 0);
107    psydecay_init(&st->psy, MAX_PERIOD/2, st->mode->Fs);
108    
109    st->in_mem = celt_alloc(N*C*sizeof(celt_sig_t));
110    st->mdct_overlap = celt_alloc(N*C*sizeof(celt_sig_t));
111    st->out_mem = celt_alloc(MAX_PERIOD*C*sizeof(celt_sig_t));
112
113    st->oldBandE = (celt_word16_t*)celt_alloc(C*mode->nbEBands*sizeof(celt_word16_t));
114
115    st->preemph = QCONST16(0.8f,15);
116    st->preemph_memE = (celt_sig_t*)celt_alloc(C*sizeof(celt_sig_t));;
117    st->preemph_memD = (celt_sig_t*)celt_alloc(C*sizeof(celt_sig_t));;
118
119    return st;
120 }
121
122 void celt_encoder_destroy(CELTEncoder *st)
123 {
124    if (st == NULL)
125    {
126       celt_warning("NULL passed to celt_encoder_destroy");
127       return;
128    }
129    if (check_mode(st->mode) != CELT_OK)
130       return;
131
132    ec_byte_writeclear(&st->buf);
133
134    kiss_fft_free(st->fft);
135    psydecay_clear(&st->psy);
136
137    celt_free(st->in_mem);
138    celt_free(st->mdct_overlap);
139    celt_free(st->out_mem);
140    
141    celt_free(st->oldBandE);
142    
143    celt_free(st->preemph_memE);
144    celt_free(st->preemph_memD);
145    
146    celt_free(st);
147 }
148
149 inline celt_int16_t SIG2INT16(celt_sig_t x)
150 {
151    x = PSHR32(x, SIG_SHIFT);
152    if (x>32767)
153       x = 32767;
154    else if (x<-32767)
155       x = -32767;
156 #ifdef FIXED_POINT
157    return EXTRACT16(x);
158 #else
159    return (celt_int16_t)floor(.5+x);
160 #endif
161 }
162
163 /** Apply window and compute the MDCT for all sub-frames and all channels in a frame */
164 static celt_word32_t compute_mdcts(const mdct_lookup *mdct_lookup, celt_word16_t *window, celt_sig_t *in, celt_sig_t *out, int N, int overlap, int B, int C)
165 {
166    int i, c, N4;
167    celt_word32_t E = 0;
168    VARDECL(celt_word32_t *x);
169    VARDECL(celt_word32_t *tmp);
170    SAVE_STACK;
171    N4 = (N-overlap)/2;
172    ALLOC(x, 2*N, celt_word32_t);
173    ALLOC(tmp, N, celt_word32_t);
174    for (c=0;c<C;c++)
175    {
176       for (i=0;i<B;i++)
177       {
178          int j;
179          for (j=0;j<2*N-2*N4;j++)
180             x[j+N4] = in[C*i*N+C*j+c];
181          for (j=0;j<overlap;j++)
182          {
183             x[j+N4] = MULT16_32_Q15(window[j],x[j+N4]);
184             x[2*N-j-N4-1] = MULT16_32_Q15(window[j],x[2*N-j-N4-1]);
185          }
186          for (j=0;j<N4;j++)
187          {
188             x[j] = 0;
189             x[2*N-j-1] = 0;
190          }
191          for (j=0;j<2*N;j++)
192             E += MULT16_16(EXTRACT16(SHR32(x[j],SIG_SHIFT+4)),EXTRACT16(SHR32(x[j],SIG_SHIFT+4)));
193          mdct_forward(mdct_lookup, x, tmp);
194          /* Interleaving the sub-frames */
195          for (j=0;j<N;j++)
196             out[C*B*j+C*i+c] = tmp[j];
197       }
198    }
199    RESTORE_STACK;
200    return E;
201 }
202
203 /** Compute the IMDCT and apply window for all sub-frames and all channels in a frame */
204 static void compute_inv_mdcts(const mdct_lookup *mdct_lookup, celt_word16_t *window, celt_sig_t *X, celt_sig_t *out_mem, celt_sig_t *mdct_overlap, int N, int overlap, int B, int C)
205 {
206    int i, c, N4;
207    VARDECL(celt_word32_t *x);
208    VARDECL(celt_word32_t *tmp);
209    SAVE_STACK;
210    ALLOC(x, 2*N, celt_word32_t);
211    ALLOC(tmp, N, celt_word32_t);
212    N4 = (N-overlap)/2;
213    for (c=0;c<C;c++)
214    {
215       for (i=0;i<B;i++)
216       {
217          int j;
218          /* De-interleaving the sub-frames */
219          for (j=0;j<N;j++)
220             tmp[j] = X[C*B*j+C*i+c];
221          mdct_backward(mdct_lookup, tmp, x);
222          /* The first and last part would need to be set to zero if we actually
223             wanted to use them. */
224          for (j=0;j<overlap;j++)
225          {
226             x[j+N4] = MULT16_32_Q15(window[j],x[j+N4]);
227             x[2*N-j-N4-1] = MULT16_32_Q15(window[j],x[2*N-j-N4-1]);
228          }
229          for (j=0;j<overlap;j++)
230             out_mem[C*(MAX_PERIOD+(i-B)*N)+C*j+c] = 2*(x[N4+j]+mdct_overlap[C*j+c]);
231          for (j=0;j<2*N4;j++)
232             out_mem[C*(MAX_PERIOD+(i-B)*N)+C*(j+overlap)+c] = 2*x[j+N4+overlap];
233          for (j=0;j<overlap;j++)
234             mdct_overlap[C*j+c] = x[N+N4+j];
235       }
236    }
237    RESTORE_STACK;
238 }
239
240 int celt_encode(CELTEncoder *st, celt_int16_t *pcm, unsigned char *compressed, int nbCompressedBytes)
241 {
242    int i, c, N, B, C, N4;
243    int has_pitch;
244    int pitch_index;
245    celt_word32_t curr_power, pitch_power;
246    VARDECL(celt_sig_t *in);
247    VARDECL(celt_sig_t *freq);
248    VARDECL(celt_norm_t *X);
249    VARDECL(celt_norm_t *P);
250    VARDECL(celt_ener_t *bandE);
251    VARDECL(celt_pgain_t *gains);
252    SAVE_STACK;
253
254    if (check_mode(st->mode) != CELT_OK)
255       return CELT_INVALID_MODE;
256
257    N = st->block_size;
258    B = st->nb_blocks;
259    C = st->mode->nbChannels;
260    N4 = (N-st->overlap)/2;
261    ALLOC(in, (B+1)*C*N-2*N4, celt_sig_t);
262    ALLOC(freq, B*C*N, celt_sig_t); /**< Interleaved signal MDCTs */
263    ALLOC(bandE,st->mode->nbEBands*C, celt_ener_t);
264    ALLOC(gains,st->mode->nbPBands, celt_pgain_t);
265    
266
267    for (c=0;c<C;c++)
268    {
269       for (i=0;i<st->overlap;i++)
270          in[C*i+c] = st->in_mem[C*i+c];
271       for (i=0;i<B*N;i++)
272       {
273          celt_sig_t tmp = SHL32(EXTEND32(pcm[C*i+c]), SIG_SHIFT);
274          in[C*(i+st->overlap)+c] = SUB32(tmp, MULT16_32_Q15(st->preemph,st->preemph_memE[c]));
275          st->preemph_memE[c] = tmp;
276       }
277       for (i=0;i<st->overlap;i++)
278          st->in_mem[C*i+c] = in[C*(N*(B+1)-2*N4-st->overlap+i)+c];
279    }
280    /*for (i=0;i<(B+1)*C*N;i++) printf ("%f(%d) ", in[i], i); printf ("\n");*/
281    /* Compute MDCTs */
282    curr_power = compute_mdcts(&st->mode->mdct, st->mode->window, in, freq, N, st->overlap, B, C);
283
284 #if 0 /* Mask disabled until it can be made to do something useful */
285    compute_mdct_masking(X, mask, B*C*N, st->Fs);
286
287    /* Invert and stretch the mask to length of X 
288       For some reason, I get better results by using the sqrt instead,
289       although there's no valid reason to. Must investigate further */
290    for (i=0;i<B*C*N;i++)
291       mask[i] = 1/(.1+mask[i]);
292 #endif
293    /* Pitch analysis */
294    for (c=0;c<C;c++)
295    {
296       for (i=0;i<st->overlap;i++)
297       {
298          in[C*i+c] = MULT16_32_Q15(st->mode->window[i], in[C*i+c]);
299          in[C*(B*N+N-i-2*N4-1)+c] = MULT16_32_Q15(st->mode->window[i], in[C*(B*N+N-i-2*N4-1)+c]);
300       }
301    }
302    find_spectral_pitch(st->fft, &st->psy, in, st->out_mem, MAX_PERIOD, (B+1)*N-2*N4, C, &pitch_index);
303    
304    /* Deferred allocation after find_spectral_pitch() to reduce the peak memory usage */
305    ALLOC(X, B*C*N, celt_norm_t);         /**< Interleaved normalised MDCTs */
306    ALLOC(P, B*C*N, celt_norm_t);         /**< Interleaved normalised pitch MDCTs*/
307
308    /*printf ("%f %f\n", curr_power, pitch_power);*/
309    /*int j;
310    for (j=0;j<B*N;j++)
311       printf ("%f ", X[j]);
312    for (j=0;j<B*N;j++)
313       printf ("%f ", P[j]);
314    printf ("\n");*/
315
316    /* Band normalisation */
317    compute_band_energies(st->mode, freq, bandE);
318    normalise_bands(st->mode, freq, X, bandE);
319    /*for (i=0;i<st->mode->nbEBands;i++)printf("%f ", bandE[i]);printf("\n");*/
320    /*for (i=0;i<N*B*C;i++)printf("%f ", X[i]);printf("\n");*/
321
322    /* Compute MDCTs of the pitch part */
323    pitch_power = compute_mdcts(&st->mode->mdct, st->mode->window, st->out_mem+pitch_index*C, freq, N, st->overlap, B, C);
324    
325
326    quant_energy(st->mode, bandE, st->oldBandE, nbCompressedBytes*8/3, &st->enc);
327
328    if (C==2)
329    {
330       stereo_mix(st->mode, X, bandE, 1);
331    }
332
333    /* Check if we can safely use the pitch (i.e. effective gain isn't too high) */
334    if (curr_power + 1e5f*(1.f/SHL16(1,8)) < 10.f*pitch_power)
335    {
336       /* Normalise the pitch vector as well (discard the energies) */
337       VARDECL(celt_ener_t *bandEp);
338       ALLOC(bandEp, st->mode->nbEBands*st->mode->nbChannels, celt_ener_t);
339       compute_band_energies(st->mode, freq, bandEp);
340       normalise_bands(st->mode, freq, P, bandEp);
341
342       if (C==2)
343          stereo_mix(st->mode, P, bandE, 1);
344       /* Simulates intensity stereo */
345       /*for (i=30;i<N*B;i++)
346          X[i*C+1] = P[i*C+1] = 0;*/
347
348       /* Pitch prediction */
349       compute_pitch_gain(st->mode, X, P, gains);
350       has_pitch = quant_pitch(gains, st->mode->nbPBands, &st->enc);
351       if (has_pitch)
352          ec_enc_uint(&st->enc, pitch_index, MAX_PERIOD-((B+1)*N-2*N4));
353    } else {
354       /* No pitch, so we just pretend we found a gain of zero */
355       for (i=0;i<st->mode->nbPBands;i++)
356          gains[i] = 0;
357       ec_enc_uint(&st->enc, 0, 128);
358       for (i=0;i<B*C*N;i++)
359          P[i] = 0;
360    }
361    
362
363    pitch_quant_bands(st->mode, P, gains);
364
365    /*for (i=0;i<B*N;i++) printf("%f ",P[i]);printf("\n");*/
366    /* Compute residual that we're going to encode */
367    for (i=0;i<B*C*N;i++)
368       X[i] -= P[i];
369
370    /* Residual quantisation */
371    quant_bands(st->mode, X, P, NULL, nbCompressedBytes*8, &st->enc);
372    
373    if (C==2)
374    {
375       stereo_mix(st->mode, X, bandE, -1);
376       renormalise_bands(st->mode, X);
377    }
378    /* Synthesis */
379    denormalise_bands(st->mode, X, freq, bandE);
380
381
382    CELT_MOVE(st->out_mem, st->out_mem+C*B*N, C*(MAX_PERIOD-B*N));
383
384    compute_inv_mdcts(&st->mode->mdct, st->mode->window, freq, st->out_mem, st->mdct_overlap, N, st->overlap, B, C);
385    /* De-emphasis and put everything back at the right place in the synthesis history */
386    for (c=0;c<C;c++)
387    {
388       for (i=0;i<B;i++)
389       {
390          int j;
391          for (j=0;j<N;j++)
392          {
393             celt_sig_t tmp = ADD32(st->out_mem[C*(MAX_PERIOD+(i-B)*N)+C*j+c],
394                                    MULT16_32_Q15(st->preemph,st->preemph_memD[c]));
395             st->preemph_memD[c] = tmp;
396             pcm[C*i*N+C*j+c] = SIG2INT16(tmp);
397          }
398       }
399    }
400    
401    if (ec_enc_tell(&st->enc, 0) < nbCompressedBytes*8 - 7)
402       celt_warning_int ("many unused bits: ", nbCompressedBytes*8-ec_enc_tell(&st->enc, 0));
403    /*printf ("%d\n", ec_enc_tell(&st->enc, 0)-8*nbCompressedBytes);*/
404    /* Finishing the stream with a 0101... pattern so that the decoder can check is everything's right */
405    {
406       int val = 0;
407       while (ec_enc_tell(&st->enc, 0) < nbCompressedBytes*8)
408       {
409          ec_enc_uint(&st->enc, val, 2);
410          val = 1-val;
411       }
412    }
413    ec_enc_done(&st->enc);
414    {
415       unsigned char *data;
416       int nbBytes = ec_byte_bytes(&st->buf);
417       if (nbBytes > nbCompressedBytes)
418       {
419          celt_warning_int ("got too many bytes:", nbBytes);
420          RESTORE_STACK;
421          return CELT_INTERNAL_ERROR;
422       }
423       /*printf ("%d\n", *nbBytes);*/
424       data = ec_byte_get_buffer(&st->buf);
425       for (i=0;i<nbBytes;i++)
426          compressed[i] = data[i];
427       for (;i<nbCompressedBytes;i++)
428          compressed[i] = 0;
429    }
430    /* Reset the packing for the next encoding */
431    ec_byte_reset(&st->buf);
432    ec_enc_init(&st->enc,&st->buf);
433
434    RESTORE_STACK;
435    return nbCompressedBytes;
436 }
437
438
439 /****************************************************************************/
440 /*                                                                          */
441 /*                                DECODER                                   */
442 /*                                                                          */
443 /****************************************************************************/
444
445
446 /** Decoder state 
447  @brief Decoder state
448  */
449 struct CELTDecoder {
450    const CELTMode *mode;
451    int frame_size;
452    int block_size;
453    int nb_blocks;
454    int overlap;
455
456    ec_byte_buffer buf;
457    ec_enc         enc;
458
459    celt_word16_t preemph;
460    celt_sig_t *preemph_memD;
461
462    celt_sig_t *mdct_overlap;
463    celt_sig_t *out_mem;
464
465    celt_word16_t *oldBandE;
466    
467    int last_pitch_index;
468 };
469
470 CELTDecoder *celt_decoder_create(const CELTMode *mode)
471 {
472    int N, B, C;
473    CELTDecoder *st;
474
475    if (check_mode(mode) != CELT_OK)
476       return NULL;
477
478    N = mode->mdctSize;
479    B = mode->nbMdctBlocks;
480    C = mode->nbChannels;
481    st = celt_alloc(sizeof(CELTDecoder));
482    
483    st->mode = mode;
484    st->frame_size = B*N;
485    st->block_size = N;
486    st->nb_blocks  = B;
487    st->overlap = mode->overlap;
488
489    st->mdct_overlap = celt_alloc(N*C*sizeof(celt_sig_t));
490    st->out_mem = celt_alloc(MAX_PERIOD*C*sizeof(celt_sig_t));
491    
492    st->oldBandE = (celt_word16_t*)celt_alloc(C*mode->nbEBands*sizeof(celt_word16_t));
493
494    st->preemph = QCONST16(0.8f,15);
495    st->preemph_memD = (celt_sig_t*)celt_alloc(C*sizeof(celt_sig_t));;
496
497    st->last_pitch_index = 0;
498    return st;
499 }
500
501 void celt_decoder_destroy(CELTDecoder *st)
502 {
503    if (st == NULL)
504    {
505       celt_warning("NULL passed to celt_encoder_destroy");
506       return;
507    }
508    if (check_mode(st->mode) != CELT_OK)
509       return;
510
511
512    celt_free(st->mdct_overlap);
513    celt_free(st->out_mem);
514    
515    celt_free(st->oldBandE);
516    
517    celt_free(st->preemph_memD);
518
519    celt_free(st);
520 }
521
522 /** Handles lost packets by just copying past data with the same offset as the last
523     pitch period */
524 static void celt_decode_lost(CELTDecoder *st, short *pcm)
525 {
526    int i, c, N, B, C, N4;
527    int pitch_index;
528    VARDECL(celt_sig_t *freq);
529    SAVE_STACK;
530    N = st->block_size;
531    B = st->nb_blocks;
532    C = st->mode->nbChannels;
533    N4 = (N-st->overlap)/2;
534    ALLOC(freq,C*B*N, celt_sig_t);         /**< Interleaved signal MDCTs */
535    
536    pitch_index = st->last_pitch_index;
537    
538    /* Use the pitch MDCT as the "guessed" signal */
539    compute_mdcts(&st->mode->mdct, st->mode->window, st->out_mem+pitch_index*C, freq, N, st->overlap, B, C);
540
541    CELT_MOVE(st->out_mem, st->out_mem+C*B*N, C*(MAX_PERIOD-B*N));
542    /* Compute inverse MDCTs */
543    compute_inv_mdcts(&st->mode->mdct, st->mode->window, freq, st->out_mem, st->mdct_overlap, N, st->overlap, B, C);
544
545    for (c=0;c<C;c++)
546    {
547       for (i=0;i<B;i++)
548       {
549          int j;
550          for (j=0;j<N;j++)
551          {
552             celt_sig_t tmp = ADD32(st->out_mem[C*(MAX_PERIOD+(i-B)*N)+C*j+c],
553                                    MULT16_32_Q15(st->preemph,st->preemph_memD[c]));
554             st->preemph_memD[c] = tmp;
555             pcm[C*i*N+C*j+c] = SIG2INT16(tmp);
556          }
557       }
558    }
559    RESTORE_STACK;
560 }
561
562 int celt_decode(CELTDecoder *st, unsigned char *data, int len, celt_int16_t *pcm)
563 {
564    int i, c, N, B, C, N4;
565    int has_pitch;
566    int pitch_index;
567    ec_dec dec;
568    ec_byte_buffer buf;
569    VARDECL(celt_sig_t *freq);
570    VARDECL(celt_norm_t *X);
571    VARDECL(celt_norm_t *P);
572    VARDECL(celt_ener_t *bandE);
573    VARDECL(celt_pgain_t *gains);
574    SAVE_STACK;
575
576    if (check_mode(st->mode) != CELT_OK)
577       return CELT_INVALID_MODE;
578
579    N = st->block_size;
580    B = st->nb_blocks;
581    C = st->mode->nbChannels;
582    N4 = (N-st->overlap)/2;
583
584    ALLOC(freq, C*B*N, celt_sig_t); /**< Interleaved signal MDCTs */
585    ALLOC(X, C*B*N, celt_norm_t);         /**< Interleaved normalised MDCTs */
586    ALLOC(P, C*B*N, celt_norm_t);         /**< Interleaved normalised pitch MDCTs*/
587    ALLOC(bandE, st->mode->nbEBands*C, celt_ener_t);
588    ALLOC(gains, st->mode->nbPBands, celt_pgain_t);
589    
590    if (check_mode(st->mode) != CELT_OK)
591    {
592       RESTORE_STACK;
593       return CELT_INVALID_MODE;
594    }
595    if (data == NULL)
596    {
597       celt_decode_lost(st, pcm);
598       RESTORE_STACK;
599       return 0;
600    }
601    
602    ec_byte_readinit(&buf,data,len);
603    ec_dec_init(&dec,&buf);
604    
605    /* Get band energies */
606    unquant_energy(st->mode, bandE, st->oldBandE, len*8/3, &dec);
607    
608    /* Get the pitch gains */
609    has_pitch = unquant_pitch(gains, st->mode->nbPBands, &dec);
610    
611    /* Get the pitch index */
612    if (has_pitch)
613    {
614       pitch_index = ec_dec_uint(&dec, MAX_PERIOD-((B+1)*N-2*N4));
615       st->last_pitch_index = pitch_index;
616    } else {
617       /* FIXME: We could be more intelligent here and just not compute the MDCT */
618       pitch_index = 0;
619    }
620    
621    /* Pitch MDCT */
622    compute_mdcts(&st->mode->mdct, st->mode->window, st->out_mem+pitch_index*C, freq, N, st->overlap, B, C);
623
624    {
625       VARDECL(celt_ener_t *bandEp);
626       ALLOC(bandEp, st->mode->nbEBands*C, celt_ener_t);
627       compute_band_energies(st->mode, freq, bandEp);
628       normalise_bands(st->mode, freq, P, bandEp);
629    }
630
631    if (C==2)
632       stereo_mix(st->mode, P, bandE, 1);
633
634    /* Apply pitch gains */
635    pitch_quant_bands(st->mode, P, gains);
636
637    /* Decode fixed codebook and merge with pitch */
638    unquant_bands(st->mode, X, P, len*8, &dec);
639
640    if (C==2)
641    {
642       stereo_mix(st->mode, X, bandE, -1);
643       renormalise_bands(st->mode, X);
644    }
645    /* Synthesis */
646    denormalise_bands(st->mode, X, freq, bandE);
647
648
649    CELT_MOVE(st->out_mem, st->out_mem+C*B*N, C*(MAX_PERIOD-B*N));
650    /* Compute inverse MDCTs */
651    compute_inv_mdcts(&st->mode->mdct, st->mode->window, freq, st->out_mem, st->mdct_overlap, N, st->overlap, B, C);
652
653    for (c=0;c<C;c++)
654    {
655       for (i=0;i<B;i++)
656       {
657          int j;
658          for (j=0;j<N;j++)
659          {
660             celt_sig_t tmp = ADD32(st->out_mem[C*(MAX_PERIOD+(i-B)*N)+C*j+c],
661                                    MULT16_32_Q15(st->preemph,st->preemph_memD[c]));
662             st->preemph_memD[c] = tmp;
663             pcm[C*i*N+C*j+c] = SIG2INT16(tmp);
664          }
665       }
666    }
667
668    {
669       unsigned int val = 0;
670       while (ec_dec_tell(&dec, 0) < len*8)
671       {
672          if (ec_dec_uint(&dec, 2) != val)
673          {
674             celt_warning("decode error");
675             RESTORE_STACK;
676             return CELT_CORRUPTED_DATA;
677          }
678          val = 1-val;
679       }
680    }
681
682    RESTORE_STACK;
683    return 0;
684    /*printf ("\n");*/
685 }
686