Pitch predictor now uses a larger range of offset by not being limited to the
[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    ALLOC(in, (B+1)*C*N, celt_sig_t);
261    ALLOC(freq, B*C*N, celt_sig_t); /**< Interleaved signal MDCTs */
262    ALLOC(bandE,st->mode->nbEBands*C, celt_ener_t);
263    ALLOC(gains,st->mode->nbPBands, celt_pgain_t);
264    
265    N4 = (N-st->overlap)/2;
266
267    for (c=0;c<C;c++)
268    {
269       for (i=0;i<N4;i++)
270          in[C*i+c] = 0;
271       for (i=0;i<st->overlap;i++)
272          in[C*(i+N4)+c] = st->in_mem[C*i+c];
273       for (i=0;i<B*N;i++)
274       {
275          celt_sig_t tmp = SHL32(EXTEND32(pcm[C*i+c]), SIG_SHIFT);
276          in[C*(i+st->overlap+N4)+c] = SUB32(tmp, MULT16_32_Q15(st->preemph,st->preemph_memE[c]));
277          st->preemph_memE[c] = tmp;
278       }
279       for (i=N*(B+1)-N4;i<N*(B+1);i++)
280          in[C*i+c] = 0;
281       for (i=0;i<st->overlap;i++)
282          st->in_mem[C*i+c] = in[C*(N*(B+1)-N4-st->overlap+i)+c];
283    }
284    /*for (i=0;i<(B+1)*C*N;i++) printf ("%f(%d) ", in[i], i); printf ("\n");*/
285    /* Compute MDCTs */
286    curr_power = compute_mdcts(&st->mode->mdct, st->mode->window, in+N4, freq, N, st->overlap, B, C);
287
288 #if 0 /* Mask disabled until it can be made to do something useful */
289    compute_mdct_masking(X, mask, B*C*N, st->Fs);
290
291    /* Invert and stretch the mask to length of X 
292       For some reason, I get better results by using the sqrt instead,
293       although there's no valid reason to. Must investigate further */
294    for (i=0;i<B*C*N;i++)
295       mask[i] = 1/(.1+mask[i]);
296 #endif
297    /* Pitch analysis */
298    for (c=0;c<C;c++)
299    {
300       for (i=0;i<N4;i++)
301       {
302          in[C*i+c] = 0;
303          in[C*(B*N+N-i-1)+c] = 0;
304       }
305       for (i=0;i<st->overlap;i++)
306       {
307          in[C*(i+N4)+c] = MULT16_32_Q15(st->mode->window[i], in[C*(i+N4)+c]);
308          in[C*(B*N+N-i-N4-1)+c] = MULT16_32_Q15(st->mode->window[i], in[C*(B*N+N-i-N4-1)+c]);
309       }
310    }
311    find_spectral_pitch(st->fft, &st->psy, in+N4, st->out_mem, MAX_PERIOD, (B+1)*N-2*N4, C, &pitch_index);
312    
313    /* Deferred allocation after find_spectral_pitch() to reduce the peak memory usage */
314    ALLOC(X, B*C*N, celt_norm_t);         /**< Interleaved normalised MDCTs */
315    ALLOC(P, B*C*N, celt_norm_t);         /**< Interleaved normalised pitch MDCTs*/
316
317    /*printf ("%f %f\n", curr_power, pitch_power);*/
318    /*int j;
319    for (j=0;j<B*N;j++)
320       printf ("%f ", X[j]);
321    for (j=0;j<B*N;j++)
322       printf ("%f ", P[j]);
323    printf ("\n");*/
324
325    /* Band normalisation */
326    compute_band_energies(st->mode, freq, bandE);
327    normalise_bands(st->mode, freq, X, bandE);
328    /*for (i=0;i<st->mode->nbEBands;i++)printf("%f ", bandE[i]);printf("\n");*/
329    /*for (i=0;i<N*B*C;i++)printf("%f ", X[i]);printf("\n");*/
330
331    /* Compute MDCTs of the pitch part */
332    pitch_power = compute_mdcts(&st->mode->mdct, st->mode->window, st->out_mem+pitch_index*C, freq, N, st->overlap, B, C);
333    
334
335    quant_energy(st->mode, bandE, st->oldBandE, nbCompressedBytes*8/3, &st->enc);
336
337    if (C==2)
338    {
339       stereo_mix(st->mode, X, bandE, 1);
340    }
341
342    /* Check if we can safely use the pitch (i.e. effective gain isn't too high) */
343    if (curr_power + 1e5f*(1.f/SHL16(1,8)) < 10.f*pitch_power)
344    {
345       /* Normalise the pitch vector as well (discard the energies) */
346       VARDECL(celt_ener_t *bandEp);
347       ALLOC(bandEp, st->mode->nbEBands*st->mode->nbChannels, celt_ener_t);
348       compute_band_energies(st->mode, freq, bandEp);
349       normalise_bands(st->mode, freq, P, bandEp);
350
351       if (C==2)
352          stereo_mix(st->mode, P, bandE, 1);
353       /* Simulates intensity stereo */
354       /*for (i=30;i<N*B;i++)
355          X[i*C+1] = P[i*C+1] = 0;*/
356
357       /* Pitch prediction */
358       compute_pitch_gain(st->mode, X, P, gains);
359       has_pitch = quant_pitch(gains, st->mode->nbPBands, &st->enc);
360       if (has_pitch)
361          ec_enc_uint(&st->enc, pitch_index, MAX_PERIOD-((B+1)*N-2*N4));
362    } else {
363       /* No pitch, so we just pretend we found a gain of zero */
364       for (i=0;i<st->mode->nbPBands;i++)
365          gains[i] = 0;
366       ec_enc_uint(&st->enc, 0, 128);
367       for (i=0;i<B*C*N;i++)
368          P[i] = 0;
369    }
370    
371
372    pitch_quant_bands(st->mode, P, gains);
373
374    /*for (i=0;i<B*N;i++) printf("%f ",P[i]);printf("\n");*/
375    /* Compute residual that we're going to encode */
376    for (i=0;i<B*C*N;i++)
377       X[i] -= P[i];
378
379    /* Residual quantisation */
380    quant_bands(st->mode, X, P, NULL, nbCompressedBytes*8, &st->enc);
381    
382    if (C==2)
383    {
384       stereo_mix(st->mode, X, bandE, -1);
385       renormalise_bands(st->mode, X);
386    }
387    /* Synthesis */
388    denormalise_bands(st->mode, X, freq, bandE);
389
390
391    CELT_MOVE(st->out_mem, st->out_mem+C*B*N, C*(MAX_PERIOD-B*N));
392
393    compute_inv_mdcts(&st->mode->mdct, st->mode->window, freq, st->out_mem, st->mdct_overlap, N, st->overlap, B, C);
394    /* De-emphasis and put everything back at the right place in the synthesis history */
395    for (c=0;c<C;c++)
396    {
397       for (i=0;i<B;i++)
398       {
399          int j;
400          for (j=0;j<N;j++)
401          {
402             celt_sig_t tmp = ADD32(st->out_mem[C*(MAX_PERIOD+(i-B)*N)+C*j+c],
403                                    MULT16_32_Q15(st->preemph,st->preemph_memD[c]));
404             st->preemph_memD[c] = tmp;
405             pcm[C*i*N+C*j+c] = SIG2INT16(tmp);
406          }
407       }
408    }
409    
410    if (ec_enc_tell(&st->enc, 0) < nbCompressedBytes*8 - 7)
411       celt_warning_int ("many unused bits: ", nbCompressedBytes*8-ec_enc_tell(&st->enc, 0));
412    /*printf ("%d\n", ec_enc_tell(&st->enc, 0)-8*nbCompressedBytes);*/
413    /* Finishing the stream with a 0101... pattern so that the decoder can check is everything's right */
414    {
415       int val = 0;
416       while (ec_enc_tell(&st->enc, 0) < nbCompressedBytes*8)
417       {
418          ec_enc_uint(&st->enc, val, 2);
419          val = 1-val;
420       }
421    }
422    ec_enc_done(&st->enc);
423    {
424       unsigned char *data;
425       int nbBytes = ec_byte_bytes(&st->buf);
426       if (nbBytes > nbCompressedBytes)
427       {
428          celt_warning_int ("got too many bytes:", nbBytes);
429          RESTORE_STACK;
430          return CELT_INTERNAL_ERROR;
431       }
432       /*printf ("%d\n", *nbBytes);*/
433       data = ec_byte_get_buffer(&st->buf);
434       for (i=0;i<nbBytes;i++)
435          compressed[i] = data[i];
436       for (;i<nbCompressedBytes;i++)
437          compressed[i] = 0;
438    }
439    /* Reset the packing for the next encoding */
440    ec_byte_reset(&st->buf);
441    ec_enc_init(&st->enc,&st->buf);
442
443    RESTORE_STACK;
444    return nbCompressedBytes;
445 }
446
447
448 /****************************************************************************/
449 /*                                                                          */
450 /*                                DECODER                                   */
451 /*                                                                          */
452 /****************************************************************************/
453
454
455 /** Decoder state 
456  @brief Decoder state
457  */
458 struct CELTDecoder {
459    const CELTMode *mode;
460    int frame_size;
461    int block_size;
462    int nb_blocks;
463    int overlap;
464
465    ec_byte_buffer buf;
466    ec_enc         enc;
467
468    celt_word16_t preemph;
469    celt_sig_t *preemph_memD;
470
471    celt_sig_t *mdct_overlap;
472    celt_sig_t *out_mem;
473
474    celt_word16_t *oldBandE;
475    
476    int last_pitch_index;
477 };
478
479 CELTDecoder *celt_decoder_create(const CELTMode *mode)
480 {
481    int N, B, C;
482    CELTDecoder *st;
483
484    if (check_mode(mode) != CELT_OK)
485       return NULL;
486
487    N = mode->mdctSize;
488    B = mode->nbMdctBlocks;
489    C = mode->nbChannels;
490    st = celt_alloc(sizeof(CELTDecoder));
491    
492    st->mode = mode;
493    st->frame_size = B*N;
494    st->block_size = N;
495    st->nb_blocks  = B;
496    st->overlap = mode->overlap;
497
498    st->mdct_overlap = celt_alloc(N*C*sizeof(celt_sig_t));
499    st->out_mem = celt_alloc(MAX_PERIOD*C*sizeof(celt_sig_t));
500    
501    st->oldBandE = (celt_word16_t*)celt_alloc(C*mode->nbEBands*sizeof(celt_word16_t));
502
503    st->preemph = QCONST16(0.8f,15);
504    st->preemph_memD = (celt_sig_t*)celt_alloc(C*sizeof(celt_sig_t));;
505
506    st->last_pitch_index = 0;
507    return st;
508 }
509
510 void celt_decoder_destroy(CELTDecoder *st)
511 {
512    if (st == NULL)
513    {
514       celt_warning("NULL passed to celt_encoder_destroy");
515       return;
516    }
517    if (check_mode(st->mode) != CELT_OK)
518       return;
519
520
521    celt_free(st->mdct_overlap);
522    celt_free(st->out_mem);
523    
524    celt_free(st->oldBandE);
525    
526    celt_free(st->preemph_memD);
527
528    celt_free(st);
529 }
530
531 /** Handles lost packets by just copying past data with the same offset as the last
532     pitch period */
533 static void celt_decode_lost(CELTDecoder *st, short *pcm)
534 {
535    int i, c, N, B, C, N4;
536    int pitch_index;
537    VARDECL(celt_sig_t *freq);
538    SAVE_STACK;
539    N = st->block_size;
540    B = st->nb_blocks;
541    C = st->mode->nbChannels;
542    N4 = (N-st->overlap)/2;
543    ALLOC(freq,C*B*N, celt_sig_t);         /**< Interleaved signal MDCTs */
544    
545    pitch_index = st->last_pitch_index;
546    
547    /* Use the pitch MDCT as the "guessed" signal */
548    compute_mdcts(&st->mode->mdct, st->mode->window, st->out_mem+pitch_index*C, freq, N, st->overlap, B, C);
549
550    CELT_MOVE(st->out_mem, st->out_mem+C*B*N, C*(MAX_PERIOD-B*N));
551    /* Compute inverse MDCTs */
552    compute_inv_mdcts(&st->mode->mdct, st->mode->window, freq, st->out_mem, st->mdct_overlap, N, st->overlap, B, C);
553
554    for (c=0;c<C;c++)
555    {
556       for (i=0;i<B;i++)
557       {
558          int j;
559          for (j=0;j<N;j++)
560          {
561             celt_sig_t tmp = ADD32(st->out_mem[C*(MAX_PERIOD+(i-B)*N)+C*j+c],
562                                    MULT16_32_Q15(st->preemph,st->preemph_memD[c]));
563             st->preemph_memD[c] = tmp;
564             pcm[C*i*N+C*j+c] = SIG2INT16(tmp);
565          }
566       }
567    }
568    RESTORE_STACK;
569 }
570
571 int celt_decode(CELTDecoder *st, unsigned char *data, int len, celt_int16_t *pcm)
572 {
573    int i, c, N, B, C, N4;
574    int has_pitch;
575    int pitch_index;
576    ec_dec dec;
577    ec_byte_buffer buf;
578    VARDECL(celt_sig_t *freq);
579    VARDECL(celt_norm_t *X);
580    VARDECL(celt_norm_t *P);
581    VARDECL(celt_ener_t *bandE);
582    VARDECL(celt_pgain_t *gains);
583    SAVE_STACK;
584
585    if (check_mode(st->mode) != CELT_OK)
586       return CELT_INVALID_MODE;
587
588    N = st->block_size;
589    B = st->nb_blocks;
590    C = st->mode->nbChannels;
591    N4 = (N-st->overlap)/2;
592
593    ALLOC(freq, C*B*N, celt_sig_t); /**< Interleaved signal MDCTs */
594    ALLOC(X, C*B*N, celt_norm_t);         /**< Interleaved normalised MDCTs */
595    ALLOC(P, C*B*N, celt_norm_t);         /**< Interleaved normalised pitch MDCTs*/
596    ALLOC(bandE, st->mode->nbEBands*C, celt_ener_t);
597    ALLOC(gains, st->mode->nbPBands, celt_pgain_t);
598    
599    if (check_mode(st->mode) != CELT_OK)
600    {
601       RESTORE_STACK;
602       return CELT_INVALID_MODE;
603    }
604    if (data == NULL)
605    {
606       celt_decode_lost(st, pcm);
607       RESTORE_STACK;
608       return 0;
609    }
610    
611    ec_byte_readinit(&buf,data,len);
612    ec_dec_init(&dec,&buf);
613    
614    /* Get band energies */
615    unquant_energy(st->mode, bandE, st->oldBandE, len*8/3, &dec);
616    
617    /* Get the pitch gains */
618    has_pitch = unquant_pitch(gains, st->mode->nbPBands, &dec);
619    
620    /* Get the pitch index */
621    if (has_pitch)
622    {
623       pitch_index = ec_dec_uint(&dec, MAX_PERIOD-((B+1)*N-2*N4));
624       st->last_pitch_index = pitch_index;
625    } else {
626       /* FIXME: We could be more intelligent here and just not compute the MDCT */
627       pitch_index = 0;
628    }
629    
630    /* Pitch MDCT */
631    compute_mdcts(&st->mode->mdct, st->mode->window, st->out_mem+pitch_index*C, freq, N, st->overlap, B, C);
632
633    {
634       VARDECL(celt_ener_t *bandEp);
635       ALLOC(bandEp, st->mode->nbEBands*C, celt_ener_t);
636       compute_band_energies(st->mode, freq, bandEp);
637       normalise_bands(st->mode, freq, P, bandEp);
638    }
639
640    if (C==2)
641       stereo_mix(st->mode, P, bandE, 1);
642
643    /* Apply pitch gains */
644    pitch_quant_bands(st->mode, P, gains);
645
646    /* Decode fixed codebook and merge with pitch */
647    unquant_bands(st->mode, X, P, len*8, &dec);
648
649    if (C==2)
650    {
651       stereo_mix(st->mode, X, bandE, -1);
652       renormalise_bands(st->mode, X);
653    }
654    /* Synthesis */
655    denormalise_bands(st->mode, X, freq, bandE);
656
657
658    CELT_MOVE(st->out_mem, st->out_mem+C*B*N, C*(MAX_PERIOD-B*N));
659    /* Compute inverse MDCTs */
660    compute_inv_mdcts(&st->mode->mdct, st->mode->window, freq, st->out_mem, st->mdct_overlap, N, st->overlap, B, C);
661
662    for (c=0;c<C;c++)
663    {
664       for (i=0;i<B;i++)
665       {
666          int j;
667          for (j=0;j<N;j++)
668          {
669             celt_sig_t tmp = ADD32(st->out_mem[C*(MAX_PERIOD+(i-B)*N)+C*j+c],
670                                    MULT16_32_Q15(st->preemph,st->preemph_memD[c]));
671             st->preemph_memD[c] = tmp;
672             pcm[C*i*N+C*j+c] = SIG2INT16(tmp);
673          }
674       }
675    }
676
677    {
678       unsigned int val = 0;
679       while (ec_dec_tell(&dec, 0) < len*8)
680       {
681          if (ec_dec_uint(&dec, 2) != val)
682          {
683             celt_warning("decode error");
684             RESTORE_STACK;
685             return CELT_CORRUPTED_DATA;
686          }
687          val = 1-val;
688       }
689    }
690
691    RESTORE_STACK;
692    return 0;
693    /*printf ("\n");*/
694 }
695