fixed-point: celt_sig_t now a 32-bit value.
[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] = 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] = 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] = 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*(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*x[j+N4+overlap];
216          for (j=0;j<overlap;j++)
217             mdct_overlap[C*j+c] = 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 = SIG_SCALING*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             tmp *= SIG_SCALING_1;
388             if (tmp > 32767) tmp = 32767;
389             if (tmp < -32767) tmp = -32767;
390             pcm[C*i*N+C*j+c] = (short)floor(.5+tmp);
391          }
392       }
393    }
394    
395    if (ec_enc_tell(&st->enc, 0) < nbCompressedBytes*8 - 7)
396       celt_warning_int ("many unused bits: ", nbCompressedBytes*8-ec_enc_tell(&st->enc, 0));
397    /*printf ("%d\n", ec_enc_tell(&st->enc, 0)-8*nbCompressedBytes);*/
398    /* Finishing the stream with a 0101... pattern so that the decoder can check is everything's right */
399    {
400       int val = 0;
401       while (ec_enc_tell(&st->enc, 0) < nbCompressedBytes*8)
402       {
403          ec_enc_uint(&st->enc, val, 2);
404          val = 1-val;
405       }
406    }
407    ec_enc_done(&st->enc);
408    {
409       unsigned char *data;
410       int nbBytes = ec_byte_bytes(&st->buf);
411       if (nbBytes > nbCompressedBytes)
412       {
413          celt_warning_int ("got too many bytes:", nbBytes);
414          return CELT_INTERNAL_ERROR;
415       }
416       /*printf ("%d\n", *nbBytes);*/
417       data = ec_byte_get_buffer(&st->buf);
418       for (i=0;i<nbBytes;i++)
419          compressed[i] = data[i];
420       for (;i<nbCompressedBytes;i++)
421          compressed[i] = 0;
422    }
423    /* Reset the packing for the next encoding */
424    ec_byte_reset(&st->buf);
425    ec_enc_init(&st->enc,&st->buf);
426
427    return nbCompressedBytes;
428 }
429
430
431 /****************************************************************************/
432 /*                                                                          */
433 /*                                DECODER                                   */
434 /*                                                                          */
435 /****************************************************************************/
436
437
438 /** Decoder state 
439  @brief Decoder state
440  */
441 struct CELTDecoder {
442    const CELTMode *mode;
443    int frame_size;
444    int block_size;
445    int nb_blocks;
446    int overlap;
447
448    ec_byte_buffer buf;
449    ec_enc         enc;
450
451    float preemph;
452    float *preemph_memD;
453    
454    mdct_lookup mdct_lookup;
455    
456    float *window;
457    celt_sig_t *mdct_overlap;
458    celt_sig_t *out_mem;
459
460    float *oldBandE;
461    
462    int last_pitch_index;
463 };
464
465 CELTDecoder *celt_decoder_create(const CELTMode *mode)
466 {
467    int i, N, B, C, N4;
468    CELTDecoder *st;
469
470    if (check_mode(mode) != CELT_OK)
471       return NULL;
472
473    N = mode->mdctSize;
474    B = mode->nbMdctBlocks;
475    C = mode->nbChannels;
476    st = celt_alloc(sizeof(CELTDecoder));
477    
478    st->mode = mode;
479    st->frame_size = B*N;
480    st->block_size = N;
481    st->nb_blocks  = B;
482    st->overlap = mode->overlap;
483
484    N4 = (N-st->overlap)/2;
485    
486    mdct_init(&st->mdct_lookup, 2*N);
487    
488    st->window = celt_alloc(2*N*sizeof(float));
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    for (i=0;i<2*N;i++)
493       st->window[i] = 0;
494    for (i=0;i<st->overlap;i++)
495       st->window[N4+i] = st->window[2*N-N4-i-1] 
496             = sin(.5*M_PI* sin(.5*M_PI*(i+.5)/st->overlap) * sin(.5*M_PI*(i+.5)/st->overlap));
497    for (i=0;i<2*N4;i++)
498       st->window[N-N4+i] = 1;
499    
500    st->oldBandE = celt_alloc(C*mode->nbEBands*sizeof(float));
501
502    st->preemph = 0.8;
503    st->preemph_memD = celt_alloc(C*sizeof(float));;
504
505    st->last_pitch_index = 0;
506    return st;
507 }
508
509 void celt_decoder_destroy(CELTDecoder *st)
510 {
511    if (st == NULL)
512    {
513       celt_warning("NULL passed to celt_encoder_destroy");
514       return;
515    }
516    if (check_mode(st->mode) != CELT_OK)
517       return;
518
519    mdct_clear(&st->mdct_lookup);
520
521    celt_free(st->window);
522    celt_free(st->mdct_overlap);
523    celt_free(st->out_mem);
524    
525    celt_free(st->oldBandE);
526    
527    celt_free(st->preemph_memD);
528
529    celt_free(st);
530 }
531
532 /** Handles lost packets by just copying past data with the same offset as the last
533     pitch period */
534 static void celt_decode_lost(CELTDecoder *st, short *pcm)
535 {
536    int i, c, N, B, C;
537    int pitch_index;
538    VARDECL(celt_sig_t *freq);
539    N = st->block_size;
540    B = st->nb_blocks;
541    C = st->mode->nbChannels;
542    ALLOC(freq,C*B*N, celt_sig_t);         /**< Interleaved signal MDCTs */
543    
544    pitch_index = st->last_pitch_index;
545    
546    /* Use the pitch MDCT as the "guessed" signal */
547    compute_mdcts(&st->mdct_lookup, st->window, st->out_mem+pitch_index*C, freq, N, B, C);
548
549    CELT_MOVE(st->out_mem, st->out_mem+C*B*N, C*(MAX_PERIOD-B*N));
550    /* Compute inverse MDCTs */
551    compute_inv_mdcts(&st->mdct_lookup, st->window, freq, st->out_mem, st->mdct_overlap, N, st->overlap, B, C);
552
553    for (c=0;c<C;c++)
554    {
555       for (i=0;i<B;i++)
556       {
557          int j;
558          for (j=0;j<N;j++)
559          {
560             float tmp = st->out_mem[C*(MAX_PERIOD+(i-B)*N)+C*j+c] + st->preemph*st->preemph_memD[c];
561             st->preemph_memD[c] = tmp;
562             tmp *= SIG_SCALING_1;
563             if (tmp > 32767) tmp = 32767;
564             if (tmp < -32767) tmp = -32767;
565             pcm[C*i*N+C*j+c] = (short)floor(.5+tmp);
566          }
567       }
568    }
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;
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(float *bandE);
582    VARDECL(float *gains);
583
584    if (check_mode(st->mode) != CELT_OK)
585       return CELT_INVALID_MODE;
586
587    N = st->block_size;
588    B = st->nb_blocks;
589    C = st->mode->nbChannels;
590    
591    ALLOC(freq, C*B*N, celt_sig_t); /**< Interleaved signal MDCTs */
592    ALLOC(X, C*B*N, celt_norm_t);         /**< Interleaved normalised MDCTs */
593    ALLOC(P, C*B*N, celt_norm_t);         /**< Interleaved normalised pitch MDCTs*/
594    ALLOC(bandE, st->mode->nbEBands*C, float);
595    ALLOC(gains, st->mode->nbPBands, float);
596    
597    if (check_mode(st->mode) != CELT_OK)
598       return CELT_INVALID_MODE;
599    if (data == NULL)
600    {
601       celt_decode_lost(st, pcm);
602       return 0;
603    }
604    
605    ec_byte_readinit(&buf,data,len);
606    ec_dec_init(&dec,&buf);
607    
608    /* Get band energies */
609    unquant_energy(st->mode, bandE, st->oldBandE, len*8/3, &dec);
610    
611    /* Get the pitch gains */
612    has_pitch = unquant_pitch(gains, st->mode->nbPBands, &dec);
613    
614    /* Get the pitch index */
615    if (has_pitch)
616    {
617       pitch_index = ec_dec_uint(&dec, MAX_PERIOD-(B+1)*N);
618       st->last_pitch_index = pitch_index;
619    } else {
620       /* FIXME: We could be more intelligent here and just not compute the MDCT */
621       pitch_index = 0;
622    }
623    
624    /* Pitch MDCT */
625    compute_mdcts(&st->mdct_lookup, st->window, st->out_mem+pitch_index*C, freq, N, B, C);
626
627    {
628       VARDECL(float *bandEp);
629       ALLOC(bandEp, st->mode->nbEBands*C, float);
630       compute_band_energies(st->mode, freq, bandEp);
631       normalise_bands(st->mode, freq, P, bandEp);
632    }
633
634    if (C==2)
635       stereo_mix(st->mode, P, bandE, 1);
636
637    /* Apply pitch gains */
638    pitch_quant_bands(st->mode, X, P, gains);
639
640    /* Decode fixed codebook and merge with pitch */
641    unquant_bands(st->mode, X, P, len*8, &dec);
642
643    if (C==2)
644    {
645       stereo_mix(st->mode, X, bandE, -1);
646       renormalise_bands(st->mode, X);
647    }
648    /* Synthesis */
649    denormalise_bands(st->mode, X, freq, bandE);
650
651
652    CELT_MOVE(st->out_mem, st->out_mem+C*B*N, C*(MAX_PERIOD-B*N));
653    /* Compute inverse MDCTs */
654    compute_inv_mdcts(&st->mdct_lookup, st->window, freq, st->out_mem, st->mdct_overlap, N, st->overlap, B, C);
655
656    for (c=0;c<C;c++)
657    {
658       for (i=0;i<B;i++)
659       {
660          int j;
661          for (j=0;j<N;j++)
662          {
663             float tmp = st->out_mem[C*(MAX_PERIOD+(i-B)*N)+C*j+c] + st->preemph*st->preemph_memD[c];
664             st->preemph_memD[c] = tmp;
665             tmp *= SIG_SCALING_1;
666             if (tmp > 32767) tmp = 32767;
667             if (tmp < -32767) tmp = -32767;
668             pcm[C*i*N+C*j+c] = (short)floor(.5+tmp);
669          }
670       }
671    }
672
673    {
674       int val = 0;
675       while (ec_dec_tell(&dec, 0) < len*8)
676       {
677          if (ec_dec_uint(&dec, 2) != val)
678          {
679             celt_warning("decode error");
680             return CELT_CORRUPTED_DATA;
681          }
682          val = 1-val;
683       }
684    }
685
686    return 0;
687    /*printf ("\n");*/
688 }
689