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