fixed-point: initial support for using the fixed-point MDCT (rest is still all
[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    float preemph;
71    float *preemph_memE;
72    float *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 = 0.8;
129    st->preemph_memE = celt_alloc(C*sizeof(float));;
130    st->preemph_memD = celt_alloc(C*sizeof(float));;
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 /** Apply window and compute the MDCT for all sub-frames and all channels in a frame */
165 static float compute_mdcts(mdct_lookup *mdct_lookup, float *window, celt_sig_t *in, celt_sig_t *out, int N, int B, int C)
166 {
167    int i, c;
168    float E = 1e-15;
169    VARDECL(celt_word32_t *x);
170    VARDECL(celt_word32_t *tmp);
171    ALLOC(x, 2*N, celt_word32_t);
172    ALLOC(tmp, N, celt_word32_t);
173    for (c=0;c<C;c++)
174    {
175       for (i=0;i<B;i++)
176       {
177          int j;
178          for (j=0;j<2*N;j++)
179          {
180             x[j] = SIG_SCALING*window[j]*in[C*i*N+C*j+c];
181             E += SIG_SCALING_1*SIG_SCALING_1*x[j]*x[j];
182          }
183          mdct_forward(mdct_lookup, x, tmp);
184          /* Interleaving the sub-frames */
185          for (j=0;j<N;j++)
186             out[C*B*j+C*i+c] = SIG_SCALING_1*tmp[j];
187       }
188    }
189    return E;
190 }
191
192 /** Compute the IMDCT and apply window for all sub-frames and all channels in a frame */
193 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)
194 {
195    int i, c, N4;
196    VARDECL(celt_word32_t *x);
197    VARDECL(celt_word32_t *tmp);
198    ALLOC(x, 2*N, celt_word32_t);
199    ALLOC(tmp, N, celt_word32_t);
200    N4 = (N-overlap)/2;
201    for (c=0;c<C;c++)
202    {
203       for (i=0;i<B;i++)
204       {
205          int j;
206          /* De-interleaving the sub-frames */
207          for (j=0;j<N;j++)
208             tmp[j] = SIG_SCALING*X[C*B*j+C*i+c];
209          mdct_backward(mdct_lookup, tmp, x);
210          for (j=0;j<2*N;j++)
211             x[j] = window[j]*x[j];
212          for (j=0;j<overlap;j++)
213             out_mem[C*(MAX_PERIOD+(i-B)*N)+C*j+c] = 2*(SIG_SCALING_1*x[N4+j]+mdct_overlap[C*j+c]);
214          for (j=0;j<2*N4;j++)
215             out_mem[C*(MAX_PERIOD+(i-B)*N)+C*(j+overlap)+c] = 2*SIG_SCALING_1*x[j+N4+overlap];
216          for (j=0;j<overlap;j++)
217             mdct_overlap[C*j+c] = SIG_SCALING_1*x[N+N4+j];
218       }
219    }
220 }
221
222 int celt_encode(CELTEncoder *st, celt_int16_t *pcm, unsigned char *compressed, int nbCompressedBytes)
223 {
224    int i, c, N, B, C, N4;
225    int has_pitch;
226    int pitch_index;
227    float curr_power, pitch_power;
228    VARDECL(celt_sig_t *in);
229    VARDECL(celt_sig_t *freq);
230    VARDECL(celt_norm_t *X);
231    VARDECL(celt_norm_t *P);
232    VARDECL(float *mask);
233    VARDECL(float *bandE);
234    VARDECL(float *gains);
235
236    if (check_mode(st->mode) != CELT_OK)
237       return CELT_INVALID_MODE;
238
239    N = st->block_size;
240    B = st->nb_blocks;
241    C = st->mode->nbChannels;
242    ALLOC(in, (B+1)*C*N, celt_sig_t);
243    ALLOC(freq, B*C*N, celt_sig_t); /**< Interleaved signal MDCTs */
244    ALLOC(X, B*C*N, celt_norm_t);         /**< Interleaved normalised MDCTs */
245    ALLOC(P, B*C*N, celt_norm_t);         /**< Interleaved normalised pitch MDCTs*/
246    ALLOC(mask, B*C*N, float);      /**< Masking curve */
247    ALLOC(bandE,st->mode->nbEBands*C, float);
248    ALLOC(gains,st->mode->nbPBands, float);
249    
250    N4 = (N-st->overlap)/2;
251
252    for (c=0;c<C;c++)
253    {
254       for (i=0;i<N4;i++)
255          in[C*i+c] = 0;
256       for (i=0;i<st->overlap;i++)
257          in[C*(i+N4)+c] = st->in_mem[C*i+c];
258       for (i=0;i<B*N;i++)
259       {
260          float tmp = pcm[C*i+c];
261          in[C*(i+st->overlap+N4)+c] = tmp - st->preemph*st->preemph_memE[c];
262          st->preemph_memE[c] = tmp;
263       }
264       for (i=N*(B+1)-N4;i<N*(B+1);i++)
265          in[C*i+c] = 0;
266       for (i=0;i<st->overlap;i++)
267          st->in_mem[C*i+c] = in[C*(N*(B+1)-N4-st->overlap+i)+c];
268    }
269    /*for (i=0;i<(B+1)*C*N;i++) printf ("%f(%d) ", in[i], i); printf ("\n");*/
270    /* Compute MDCTs */
271    curr_power = compute_mdcts(&st->mdct_lookup, st->window, in, freq, N, B, C);
272
273 #if 0 /* Mask disabled until it can be made to do something useful */
274    compute_mdct_masking(X, mask, B*C*N, st->Fs);
275
276    /* Invert and stretch the mask to length of X 
277       For some reason, I get better results by using the sqrt instead,
278       although there's no valid reason to. Must investigate further */
279    for (i=0;i<B*C*N;i++)
280       mask[i] = 1/(.1+mask[i]);
281 #else
282    for (i=0;i<B*C*N;i++)
283       mask[i] = 1;
284 #endif
285    /* Pitch analysis */
286    for (c=0;c<C;c++)
287    {
288       for (i=0;i<N;i++)
289       {
290          in[C*i+c] *= st->window[i];
291          in[C*(B*N+i)+c] *= st->window[N+i];
292       }
293    }
294    find_spectral_pitch(st->fft, &st->psy, in, st->out_mem, MAX_PERIOD, (B+1)*N, C, &pitch_index);
295    
296    /*printf ("%f %f\n", curr_power, pitch_power);*/
297    /*int j;
298    for (j=0;j<B*N;j++)
299       printf ("%f ", X[j]);
300    for (j=0;j<B*N;j++)
301       printf ("%f ", P[j]);
302    printf ("\n");*/
303
304    /* Band normalisation */
305    compute_band_energies(st->mode, freq, bandE);
306    normalise_bands(st->mode, freq, X, bandE);
307    /*for (i=0;i<st->mode->nbEBands;i++)printf("%f ", bandE[i]);printf("\n");*/
308    /*for (i=0;i<N*B*C;i++)printf("%f ", X[i]);printf("\n");*/
309
310    /* Compute MDCTs of the pitch part */
311    pitch_power = compute_mdcts(&st->mdct_lookup, st->window, st->out_mem+pitch_index*C, freq, N, B, C);
312    
313
314    quant_energy(st->mode, bandE, st->oldBandE, nbCompressedBytes*8/3, &st->enc);
315
316    if (C==2)
317    {
318       stereo_mix(st->mode, X, bandE, 1);
319    }
320
321    /* Check if we can safely use the pitch (i.e. effective gain isn't too high) */
322    if (curr_power + 1e5f < 10.f*pitch_power)
323    {
324       /* Normalise the pitch vector as well (discard the energies) */
325       VARDECL(float *bandEp);
326       ALLOC(bandEp, st->mode->nbEBands*st->mode->nbChannels, float);
327       compute_band_energies(st->mode, freq, bandEp);
328       normalise_bands(st->mode, freq, P, bandEp);
329
330       if (C==2)
331          stereo_mix(st->mode, P, bandE, 1);
332       /* Simulates intensity stereo */
333       /*for (i=30;i<N*B;i++)
334          X[i*C+1] = P[i*C+1] = 0;*/
335
336       /* Pitch prediction */
337       compute_pitch_gain(st->mode, X, P, gains, bandE);
338       has_pitch = quant_pitch(gains, st->mode->nbPBands, &st->enc);
339       if (has_pitch)
340          ec_enc_uint(&st->enc, pitch_index, MAX_PERIOD-(B+1)*N);
341    } else {
342       /* No pitch, so we just pretend we found a gain of zero */
343       for (i=0;i<st->mode->nbPBands;i++)
344          gains[i] = 0;
345       ec_enc_uint(&st->enc, 0, 128);
346       for (i=0;i<B*C*N;i++)
347          P[i] = 0;
348    }
349    
350
351    pitch_quant_bands(st->mode, X, P, gains);
352
353    /*for (i=0;i<B*N;i++) printf("%f ",P[i]);printf("\n");*/
354    /* Compute residual that we're going to encode */
355    for (i=0;i<B*C*N;i++)
356       X[i] -= P[i];
357
358    /*float sum=0;
359    for (i=0;i<B*N;i++)
360       sum += X[i]*X[i];
361    printf ("%f\n", sum);*/
362    /* Residual quantisation */
363    quant_bands(st->mode, X, P, mask, nbCompressedBytes*8, &st->enc);
364    
365    if (C==2)
366    {
367       stereo_mix(st->mode, X, bandE, -1);
368       renormalise_bands(st->mode, X);
369    }
370    /* Synthesis */
371    denormalise_bands(st->mode, X, freq, bandE);
372
373
374    CELT_MOVE(st->out_mem, st->out_mem+C*B*N, C*(MAX_PERIOD-B*N));
375
376    compute_inv_mdcts(&st->mdct_lookup, st->window, freq, st->out_mem, st->mdct_overlap, N, st->overlap, B, C);
377    /* De-emphasis and put everything back at the right place in the synthesis history */
378    for (c=0;c<C;c++)
379    {
380       for (i=0;i<B;i++)
381       {
382          int j;
383          for (j=0;j<N;j++)
384          {
385             float tmp = st->out_mem[C*(MAX_PERIOD+(i-B)*N)+C*j+c] + st->preemph*st->preemph_memD[c];
386             st->preemph_memD[c] = tmp;
387             if (tmp > 32767) tmp = 32767;
388             if (tmp < -32767) tmp = -32767;
389             pcm[C*i*N+C*j+c] = (short)floor(.5+tmp);
390          }
391       }
392    }
393    
394    if (ec_enc_tell(&st->enc, 0) < nbCompressedBytes*8 - 7)
395       celt_warning_int ("many unused bits: ", nbCompressedBytes*8-ec_enc_tell(&st->enc, 0));
396    /*printf ("%d\n", ec_enc_tell(&st->enc, 0)-8*nbCompressedBytes);*/
397    /* Finishing the stream with a 0101... pattern so that the decoder can check is everything's right */
398    {
399       int val = 0;
400       while (ec_enc_tell(&st->enc, 0) < nbCompressedBytes*8)
401       {
402          ec_enc_uint(&st->enc, val, 2);
403          val = 1-val;
404       }
405    }
406    ec_enc_done(&st->enc);
407    {
408       unsigned char *data;
409       int nbBytes = ec_byte_bytes(&st->buf);
410       if (nbBytes > nbCompressedBytes)
411       {
412          celt_warning_int ("got too many bytes:", nbBytes);
413          return CELT_INTERNAL_ERROR;
414       }
415       /*printf ("%d\n", *nbBytes);*/
416       data = ec_byte_get_buffer(&st->buf);
417       for (i=0;i<nbBytes;i++)
418          compressed[i] = data[i];
419       for (;i<nbCompressedBytes;i++)
420          compressed[i] = 0;
421    }
422    /* Reset the packing for the next encoding */
423    ec_byte_reset(&st->buf);
424    ec_enc_init(&st->enc,&st->buf);
425
426    return nbCompressedBytes;
427 }
428
429
430 /****************************************************************************/
431 /*                                                                          */
432 /*                                DECODER                                   */
433 /*                                                                          */
434 /****************************************************************************/
435
436
437 /** Decoder state 
438  @brief Decoder state
439  */
440 struct CELTDecoder {
441    const CELTMode *mode;
442    int frame_size;
443    int block_size;
444    int nb_blocks;
445    int overlap;
446
447    ec_byte_buffer buf;
448    ec_enc         enc;
449
450    float preemph;
451    float *preemph_memD;
452    
453    mdct_lookup mdct_lookup;
454    
455    float *window;
456    celt_sig_t *mdct_overlap;
457    celt_sig_t *out_mem;
458
459    float *oldBandE;
460    
461    int last_pitch_index;
462 };
463
464 CELTDecoder *celt_decoder_create(const CELTMode *mode)
465 {
466    int i, N, B, C, N4;
467    CELTDecoder *st;
468
469    if (check_mode(mode) != CELT_OK)
470       return NULL;
471
472    N = mode->mdctSize;
473    B = mode->nbMdctBlocks;
474    C = mode->nbChannels;
475    st = celt_alloc(sizeof(CELTDecoder));
476    
477    st->mode = mode;
478    st->frame_size = B*N;
479    st->block_size = N;
480    st->nb_blocks  = B;
481    st->overlap = mode->overlap;
482
483    N4 = (N-st->overlap)/2;
484    
485    mdct_init(&st->mdct_lookup, 2*N);
486    
487    st->window = celt_alloc(2*N*sizeof(float));
488    st->mdct_overlap = celt_alloc(N*C*sizeof(celt_sig_t));
489    st->out_mem = celt_alloc(MAX_PERIOD*C*sizeof(celt_sig_t));
490
491    for (i=0;i<2*N;i++)
492       st->window[i] = 0;
493    for (i=0;i<st->overlap;i++)
494       st->window[N4+i] = st->window[2*N-N4-i-1] 
495             = sin(.5*M_PI* sin(.5*M_PI*(i+.5)/st->overlap) * sin(.5*M_PI*(i+.5)/st->overlap));
496    for (i=0;i<2*N4;i++)
497       st->window[N-N4+i] = 1;
498    
499    st->oldBandE = celt_alloc(C*mode->nbEBands*sizeof(float));
500
501    st->preemph = 0.8;
502    st->preemph_memD = celt_alloc(C*sizeof(float));;
503
504    st->last_pitch_index = 0;
505    return st;
506 }
507
508 void celt_decoder_destroy(CELTDecoder *st)
509 {
510    if (st == NULL)
511    {
512       celt_warning("NULL passed to celt_encoder_destroy");
513       return;
514    }
515    if (check_mode(st->mode) != CELT_OK)
516       return;
517
518    mdct_clear(&st->mdct_lookup);
519
520    celt_free(st->window);
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;
536    int pitch_index;
537    VARDECL(celt_sig_t *freq);
538    N = st->block_size;
539    B = st->nb_blocks;
540    C = st->mode->nbChannels;
541    ALLOC(freq,C*B*N, celt_sig_t);         /**< Interleaved signal MDCTs */
542    
543    pitch_index = st->last_pitch_index;
544    
545    /* Use the pitch MDCT as the "guessed" signal */
546    compute_mdcts(&st->mdct_lookup, st->window, st->out_mem+pitch_index*C, freq, N, B, C);
547
548    CELT_MOVE(st->out_mem, st->out_mem+C*B*N, C*(MAX_PERIOD-B*N));
549    /* Compute inverse MDCTs */
550    compute_inv_mdcts(&st->mdct_lookup, st->window, freq, st->out_mem, st->mdct_overlap, N, st->overlap, B, C);
551
552    for (c=0;c<C;c++)
553    {
554       for (i=0;i<B;i++)
555       {
556          int j;
557          for (j=0;j<N;j++)
558          {
559             float tmp = st->out_mem[C*(MAX_PERIOD+(i-B)*N)+C*j+c] + st->preemph*st->preemph_memD[c];
560             st->preemph_memD[c] = tmp;
561             if (tmp > 32767) tmp = 32767;
562             if (tmp < -32767) tmp = -32767;
563             pcm[C*i*N+C*j+c] = (short)floor(.5+tmp);
564          }
565       }
566    }
567 }
568
569 int celt_decode(CELTDecoder *st, unsigned char *data, int len, celt_int16_t *pcm)
570 {
571    int i, c, N, B, C;
572    int has_pitch;
573    int pitch_index;
574    ec_dec dec;
575    ec_byte_buffer buf;
576    VARDECL(celt_sig_t *freq);
577    VARDECL(celt_norm_t *X);
578    VARDECL(celt_norm_t *P);
579    VARDECL(float *bandE);
580    VARDECL(float *gains);
581
582    if (check_mode(st->mode) != CELT_OK)
583       return CELT_INVALID_MODE;
584
585    N = st->block_size;
586    B = st->nb_blocks;
587    C = st->mode->nbChannels;
588    
589    ALLOC(freq, C*B*N, celt_sig_t); /**< Interleaved signal MDCTs */
590    ALLOC(X, C*B*N, celt_norm_t);         /**< Interleaved normalised MDCTs */
591    ALLOC(P, C*B*N, celt_norm_t);         /**< Interleaved normalised pitch MDCTs*/
592    ALLOC(bandE, st->mode->nbEBands*C, float);
593    ALLOC(gains, st->mode->nbPBands, float);
594    
595    if (check_mode(st->mode) != CELT_OK)
596       return CELT_INVALID_MODE;
597    if (data == NULL)
598    {
599       celt_decode_lost(st, pcm);
600       return 0;
601    }
602    
603    ec_byte_readinit(&buf,data,len);
604    ec_dec_init(&dec,&buf);
605    
606    /* Get band energies */
607    unquant_energy(st->mode, bandE, st->oldBandE, len*8/3, &dec);
608    
609    /* Get the pitch gains */
610    has_pitch = unquant_pitch(gains, st->mode->nbPBands, &dec);
611    
612    /* Get the pitch index */
613    if (has_pitch)
614    {
615       pitch_index = ec_dec_uint(&dec, MAX_PERIOD-(B+1)*N);
616       st->last_pitch_index = pitch_index;
617    } else {
618       /* FIXME: We could be more intelligent here and just not compute the MDCT */
619       pitch_index = 0;
620    }
621    
622    /* Pitch MDCT */
623    compute_mdcts(&st->mdct_lookup, st->window, st->out_mem+pitch_index*C, freq, N, B, C);
624
625    {
626       VARDECL(float *bandEp);
627       ALLOC(bandEp, st->mode->nbEBands*C, float);
628       compute_band_energies(st->mode, freq, bandEp);
629       normalise_bands(st->mode, freq, P, bandEp);
630    }
631
632    if (C==2)
633       stereo_mix(st->mode, P, bandE, 1);
634
635    /* Apply pitch gains */
636    pitch_quant_bands(st->mode, X, P, gains);
637
638    /* Decode fixed codebook and merge with pitch */
639    unquant_bands(st->mode, X, P, len*8, &dec);
640
641    if (C==2)
642    {
643       stereo_mix(st->mode, X, bandE, -1);
644       renormalise_bands(st->mode, X);
645    }
646    /* Synthesis */
647    denormalise_bands(st->mode, X, freq, bandE);
648
649
650    CELT_MOVE(st->out_mem, st->out_mem+C*B*N, C*(MAX_PERIOD-B*N));
651    /* Compute inverse MDCTs */
652    compute_inv_mdcts(&st->mdct_lookup, st->window, freq, st->out_mem, st->mdct_overlap, N, st->overlap, B, C);
653
654    for (c=0;c<C;c++)
655    {
656       for (i=0;i<B;i++)
657       {
658          int j;
659          for (j=0;j<N;j++)
660          {
661             float tmp = st->out_mem[C*(MAX_PERIOD+(i-B)*N)+C*j+c] + st->preemph*st->preemph_memD[c];
662             st->preemph_memD[c] = tmp;
663             if (tmp > 32767) tmp = 32767;
664             if (tmp < -32767) tmp = -32767;
665             pcm[C*i*N+C*j+c] = (short)floor(.5+tmp);
666          }
667       }
668    }
669
670    {
671       int val = 0;
672       while (ec_dec_tell(&dec, 0) < len*8)
673       {
674          if (ec_dec_uint(&dec, 2) != val)
675          {
676             celt_warning("decode error");
677             return CELT_CORRUPTED_DATA;
678          }
679          val = 1-val;
680       }
681    }
682
683    return 0;
684    /*printf ("\n");*/
685 }
686