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