Implemented two pre-echo avoidance techniques: time-domain pre-emphasis and
[opus.git] / libcelt / celt.c
1 /* (C) 2007-2008 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 #include "stack_alloc.h"
52
53 static const celt_word16_t preemph = QCONST16(0.8f,15);
54
55 static const float gainWindow[16] = {
56    0.0085135, 0.0337639, 0.0748914, 0.1304955, 0.1986827, 0.2771308, 0.3631685, 0.4538658,
57    0.5461342, 0.6368315, 0.7228692, 0.8013173, 0.8695045, 0.9251086, 0.9662361, 0.9914865};
58    
59 /** Encoder state 
60  @brief Encoder state
61  */
62 struct CELTEncoder {
63    const CELTMode *mode;     /**< Mode used by the encoder */
64    int frame_size;
65    int block_size;
66    int overlap;
67    int channels;
68    
69    ec_byte_buffer buf;
70    ec_enc         enc;
71
72    celt_word16_t * restrict preemph_memE; /* Input is 16-bit, so why bother with 32 */
73    celt_sig_t    * restrict preemph_memD;
74
75    celt_sig_t *in_mem;
76    celt_sig_t *out_mem;
77
78    celt_word16_t *oldBandE;
79 #ifdef EXP_PSY
80    celt_word16_t *psy_mem;
81    struct PsyDecay psy;
82 #endif
83 };
84
85 CELTEncoder *celt_encoder_create(const CELTMode *mode)
86 {
87    int N, C;
88    CELTEncoder *st;
89
90    if (check_mode(mode) != CELT_OK)
91       return NULL;
92
93    N = mode->mdctSize;
94    C = mode->nbChannels;
95    st = celt_alloc(sizeof(CELTEncoder));
96    
97    st->mode = mode;
98    st->frame_size = N;
99    st->block_size = N;
100    st->overlap = mode->overlap;
101
102    ec_byte_writeinit(&st->buf);
103    ec_enc_init(&st->enc,&st->buf);
104
105    st->in_mem = celt_alloc(st->overlap*C*sizeof(celt_sig_t));
106    st->out_mem = celt_alloc((MAX_PERIOD+st->overlap)*C*sizeof(celt_sig_t));
107
108    st->oldBandE = (celt_word16_t*)celt_alloc(C*mode->nbEBands*sizeof(celt_word16_t));
109
110    st->preemph_memE = (celt_word16_t*)celt_alloc(C*sizeof(celt_word16_t));;
111    st->preemph_memD = (celt_sig_t*)celt_alloc(C*sizeof(celt_sig_t));;
112
113 #ifdef EXP_PSY
114    st->psy_mem = celt_alloc(MAX_PERIOD*sizeof(celt_word16_t));
115    psydecay_init(&st->psy, MAX_PERIOD/2, st->mode->Fs);
116 #endif
117
118    return st;
119 }
120
121 void celt_encoder_destroy(CELTEncoder *st)
122 {
123    if (st == NULL)
124    {
125       celt_warning("NULL passed to celt_encoder_destroy");
126       return;
127    }
128    if (check_mode(st->mode) != CELT_OK)
129       return;
130
131    ec_byte_writeclear(&st->buf);
132
133    celt_free(st->in_mem);
134    celt_free(st->out_mem);
135    
136    celt_free(st->oldBandE);
137    
138    celt_free(st->preemph_memE);
139    celt_free(st->preemph_memD);
140    
141 #ifdef EXP_PSY
142    celt_free (st->psy_mem);
143    psydecay_clear(&st->psy);
144 #endif
145    
146    celt_free(st);
147 }
148
149 static inline celt_int16_t SIG2INT16(celt_sig_t x)
150 {
151    x = PSHR32(x, SIG_SHIFT);
152    x = MAX32(x, -32768);
153    x = MIN32(x, 32767);
154 #ifdef FIXED_POINT
155    return EXTRACT16(x);
156 #else
157    return (celt_int16_t)floor(.5+x);
158 #endif
159 }
160
161 /** Apply window and compute the MDCT for all sub-frames and all channels in a frame */
162 static void compute_mdcts(const CELTMode *mode, const celt_word16_t * restrict window, celt_sig_t * restrict in, celt_sig_t * restrict out)
163 {
164    const mdct_lookup *lookup = MDCT(mode);
165    const int N = FRAMESIZE(mode);
166    const int C = CHANNELS(mode);
167    const int overlap = OVERLAP(mode);
168    if (C==1)
169    {
170       mdct_forward(lookup, in, out, window, overlap);
171    } else {
172       int c;
173       VARDECL(celt_word32_t, x);
174       VARDECL(celt_word32_t, tmp);
175       SAVE_STACK;
176       ALLOC(x, N+overlap, celt_word32_t);
177       ALLOC(tmp, N, celt_word32_t);
178       for (c=0;c<C;c++)
179       {
180          int j;
181          for (j=0;j<N+overlap;j++)
182             x[j] = in[C*j+c];
183          mdct_forward(lookup, x, tmp, window, overlap);
184          /* Interleaving the sub-frames */
185          for (j=0;j<N;j++)
186             out[C*j+c] = tmp[j];
187       }
188       RESTORE_STACK;
189    }
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(const CELTMode *mode, const celt_word16_t * restrict window, celt_sig_t *X, int transient_time, float transient_gain, celt_sig_t * restrict out_mem)
194 {
195    int c, N4;
196    const int C = CHANNELS(mode);
197    const mdct_lookup *lookup = MDCT(mode);
198    const int N = FRAMESIZE(mode);
199    const int overlap = OVERLAP(mode);
200    N4 = (N-overlap)>>1;
201    for (c=0;c<C;c++)
202    {
203       int j;
204       if (transient_time<0 && C==1) {
205          mdct_backward(lookup, X, out_mem+C*(MAX_PERIOD-N-N4), window, overlap);
206       } else {
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          /* De-interleaving the sub-frames */
213          for (j=0;j<N;j++)
214             tmp[j] = X[C*j+c];
215          /* Prevents problems from the imdct doing the overlap-add */
216          CELT_MEMSET(x+N4, 0, overlap);
217          mdct_backward(lookup, tmp, x, window, overlap);
218          if (transient_time >= 0)
219          {
220             for (j=0;j<16;j++)
221                x[N4+transient_time+j-16] *= 1+gainWindow[j]*(transient_gain-1);
222             for (j=transient_time;j<N+overlap;j++)
223                x[N4+j] *= transient_gain;
224          }
225          /* The first and last part would need to be set to zero if we actually
226          wanted to use them. */
227          for (j=0;j<overlap;j++)
228             out_mem[C*(MAX_PERIOD-N)+C*j+c] += x[j+N4];
229          for (j=0;j<overlap;j++)
230             out_mem[C*(MAX_PERIOD)+C*(overlap-j-1)+c] = x[2*N-j-N4-1];
231          for (j=0;j<2*N4;j++)
232             out_mem[C*(MAX_PERIOD-N)+C*(j+overlap)+c] = x[j+N4+overlap];
233          RESTORE_STACK;
234       }
235    }
236 }
237
238 int celt_encode(CELTEncoder * restrict st, celt_int16_t * restrict pcm, unsigned char *compressed, int nbCompressedBytes)
239 {
240    int i, c, N, N4;
241    int has_pitch;
242    int pitch_index;
243    celt_word32_t curr_power, pitch_power;
244    VARDECL(celt_sig_t, in);
245    VARDECL(celt_sig_t, freq);
246    VARDECL(celt_norm_t, X);
247    VARDECL(celt_norm_t, P);
248    VARDECL(celt_ener_t, bandE);
249    VARDECL(celt_pgain_t, gains);
250    VARDECL(int, stereo_mode);
251 #ifdef EXP_PSY
252    VARDECL(celt_word32_t, mask);
253 #endif
254    int time_domain=0;
255    int transient_time;
256    float transient_gain;
257    const int C = CHANNELS(st->mode);
258    SAVE_STACK;
259
260    if (check_mode(st->mode) != CELT_OK)
261       return CELT_INVALID_MODE;
262
263    N = st->block_size;
264    N4 = (N-st->overlap)>>1;
265    ALLOC(in, 2*C*N-2*C*N4, celt_sig_t);
266
267    CELT_COPY(in, st->in_mem, C*st->overlap);
268    for (c=0;c<C;c++)
269    {
270       const celt_int16_t * restrict pcmp = pcm+c;
271       celt_sig_t * restrict inp = in+C*st->overlap+c;
272       for (i=0;i<N;i++)
273       {
274          /* Apply pre-emphasis */
275          celt_sig_t tmp = SHL32(EXTEND32(*pcmp), SIG_SHIFT);
276          *inp = SUB32(tmp, SHR32(MULT16_16(preemph,st->preemph_memE[c]),1));
277          st->preemph_memE[c] = *pcmp;
278          inp += C;
279          pcmp += C;
280       }
281    }
282    CELT_COPY(st->in_mem, in+C*(2*N-2*N4-st->overlap), C*st->overlap);
283    if (1) {
284       int len = N+st->overlap;
285       float maxR, maxD;
286       float begin[C*len], end[C*len];
287       begin[0] = in[0]*in[0];
288       for (i=1;i<len*C;i++)
289          begin[i] = begin[i-1]+in[i]*in[i];
290       end[len-1] = in[len-1]*in[len-1];
291       for (i=C*len-2;i>=0;i--)
292          end[i] = end[i+1] + in[i]*in[i];
293       maxD = 0;
294       maxR = 1;
295       transient_time = -1;
296       for (i=8*C;i<C*(len-8);i++)
297       {
298          float diff = sqrt(sqrt(end[i]/(C*len-i)))-sqrt(sqrt(begin[i]/(i)));
299          float ratio = ((1000+end[i])*i)/((1000+begin[i])*(C*len-i));
300          if (diff > maxD)
301          {
302             maxD = diff;
303             maxR = ratio;
304             transient_time = i;
305          }
306       }
307       transient_time /= C;
308       if (transient_time<32)
309       {
310          transient_time = -1;
311          maxR = 0;
312       }
313       if (maxR > 20)
314       {
315          float gain_1;
316          ec_enc_bits(&st->enc, 1, 1);
317          if (maxR < 30)
318          {
319             transient_gain = 1;
320             ec_enc_bits(&st->enc, 0, 2);
321          } else if (maxR < 100)
322          {
323             transient_gain = 2;
324             ec_enc_bits(&st->enc, 1, 2);
325          } else if (maxR < 500)
326          {
327             transient_gain = 4;
328             ec_enc_bits(&st->enc, 2, 2);
329          } else
330          {
331             transient_gain = 8;
332             ec_enc_bits(&st->enc, 3, 2);
333          }
334          ec_enc_uint(&st->enc, transient_time, len);
335          for (c=0;c<C;c++)
336             for (i=0;i<16;i++)
337                in[C*(transient_time+i-16)+c] /= 1+gainWindow[i]*(transient_gain-1);
338          gain_1 = 1./transient_gain;
339          for (c=0;c<C;c++)
340             for (i=transient_time;i<len;i++)
341                in[C*i+c] *= gain_1;
342          time_domain = 1;
343       } else {
344          ec_enc_bits(&st->enc, 0, 1);
345          transient_time = -1;
346          transient_gain = 1;
347          time_domain = 0;
348       }
349    }
350    /* Pitch analysis: we do it early to save on the peak stack space */
351    if (!time_domain)
352       find_spectral_pitch(st->mode, st->mode->fft, &st->mode->psy, in, st->out_mem, st->mode->window, 2*N-2*N4, MAX_PERIOD-(2*N-2*N4), &pitch_index);
353
354    ALLOC(freq, C*N, celt_sig_t); /**< Interleaved signal MDCTs */
355    
356    /*for (i=0;i<(B+1)*C*N;i++) printf ("%f(%d) ", in[i], i); printf ("\n");*/
357    /* Compute MDCTs */
358    compute_mdcts(st->mode, st->mode->window, in, freq);
359
360 #ifdef EXP_PSY
361    CELT_MOVE(st->psy_mem, st->out_mem+N, MAX_PERIOD+st->overlap-N);
362    for (i=0;i<N;i++)
363       st->psy_mem[MAX_PERIOD+st->overlap-N+i] = in[C*(st->overlap+i)];
364    for (c=1;c<C;c++)
365       for (i=0;i<N;i++)
366          st->psy_mem[MAX_PERIOD+st->overlap-N+i] += in[C*(st->overlap+i)+c];
367
368    ALLOC(mask, N, celt_sig_t);
369    compute_mdct_masking(&st->psy, freq, st->psy_mem, mask, C*N);
370
371    /* Invert and stretch the mask to length of X 
372       For some reason, I get better results by using the sqrt instead,
373       although there's no valid reason to. Must investigate further */
374    for (i=0;i<C*N;i++)
375       mask[i] = 1/(.1+mask[i]);
376 #endif
377    
378    /* Deferred allocation after find_spectral_pitch() to reduce the peak memory usage */
379    ALLOC(X, C*N, celt_norm_t);         /**< Interleaved normalised MDCTs */
380    ALLOC(P, C*N, celt_norm_t);         /**< Interleaved normalised pitch MDCTs*/
381    ALLOC(bandE,st->mode->nbEBands*C, celt_ener_t);
382    ALLOC(gains,st->mode->nbPBands, celt_pgain_t);
383
384    /*printf ("%f %f\n", curr_power, pitch_power);*/
385    /*int j;
386    for (j=0;j<B*N;j++)
387       printf ("%f ", X[j]);
388    for (j=0;j<B*N;j++)
389       printf ("%f ", P[j]);
390    printf ("\n");*/
391
392    /* Band normalisation */
393    compute_band_energies(st->mode, freq, bandE);
394    normalise_bands(st->mode, freq, X, bandE);
395    /*for (i=0;i<st->mode->nbEBands;i++)printf("%f ", bandE[i]);printf("\n");*/
396    /*for (i=0;i<N*B*C;i++)printf("%f ", X[i]);printf("\n");*/
397
398    /* Compute MDCTs of the pitch part */
399    if (!time_domain)
400       compute_mdcts(st->mode, st->mode->window, st->out_mem+pitch_index*C, freq);
401
402    {
403       /* Normalise the pitch vector as well (discard the energies) */
404       VARDECL(celt_ener_t, bandEp);
405       ALLOC(bandEp, st->mode->nbEBands*st->mode->nbChannels, celt_ener_t);
406       compute_band_energies(st->mode, freq, bandEp);
407       normalise_bands(st->mode, freq, P, bandEp);
408       pitch_power = bandEp[0]+bandEp[1]+bandEp[2];
409    }
410    curr_power = bandE[0]+bandE[1]+bandE[2];
411    /* Check if we can safely use the pitch (i.e. effective gain isn't too high) */
412    if (!time_domain && (MULT16_32_Q15(QCONST16(.1f, 15),curr_power) + QCONST32(10.f,ENER_SHIFT) < pitch_power))
413    {
414       /* Simulates intensity stereo */
415       /*for (i=30;i<N*B;i++)
416          X[i*C+1] = P[i*C+1] = 0;*/
417
418       /* Pitch prediction */
419       compute_pitch_gain(st->mode, X, P, gains);
420       has_pitch = quant_pitch(gains, st->mode->nbPBands, &st->enc);
421       if (has_pitch)
422          ec_enc_uint(&st->enc, pitch_index, MAX_PERIOD-(2*N-2*N4));
423    } else {
424       /* No pitch, so we just pretend we found a gain of zero */
425       for (i=0;i<st->mode->nbPBands;i++)
426          gains[i] = 0;
427       ec_enc_bits(&st->enc, 0, 7);
428       for (i=0;i<C*N;i++)
429          P[i] = 0;
430    }
431    quant_energy(st->mode, bandE, st->oldBandE, 20*C+nbCompressedBytes*8/5, st->mode->prob, &st->enc);
432
433    ALLOC(stereo_mode, st->mode->nbEBands, int);
434    stereo_decision(st->mode, X, stereo_mode, st->mode->nbEBands);
435
436    pitch_quant_bands(st->mode, P, gains);
437
438    /*for (i=0;i<B*N;i++) printf("%f ",P[i]);printf("\n");*/
439
440    /* Residual quantisation */
441    quant_bands(st->mode, X, P, NULL, bandE, stereo_mode, nbCompressedBytes*8, time_domain, &st->enc);
442    
443    if (C==2)
444    {
445       renormalise_bands(st->mode, X);
446    }
447    /* Synthesis */
448    denormalise_bands(st->mode, X, freq, bandE);
449
450
451    CELT_MOVE(st->out_mem, st->out_mem+C*N, C*(MAX_PERIOD+st->overlap-N));
452
453    compute_inv_mdcts(st->mode, st->mode->window, freq, transient_time, transient_gain, st->out_mem);
454    /* De-emphasis and put everything back at the right place in the synthesis history */
455 #ifndef SHORTCUTS
456    for (c=0;c<C;c++)
457    {
458       int j;
459       celt_sig_t * restrict outp=st->out_mem+C*(MAX_PERIOD-N)+c;
460       celt_int16_t * restrict pcmp = pcm+c;
461       for (j=0;j<N;j++)
462       {
463          celt_sig_t tmp = ADD32(*outp, MULT16_32_Q15(preemph,st->preemph_memD[c]));
464          st->preemph_memD[c] = tmp;
465          *pcmp = SIG2INT16(tmp);
466          pcmp += C;
467          outp += C;
468       }
469    }
470 #endif
471    if (ec_enc_tell(&st->enc, 0) < nbCompressedBytes*8 - 7)
472       celt_warning_int ("many unused bits: ", nbCompressedBytes*8-ec_enc_tell(&st->enc, 0));
473    /*printf ("%d\n", ec_enc_tell(&st->enc, 0)-8*nbCompressedBytes);*/
474    /* Finishing the stream with a 0101... pattern so that the decoder can check is everything's right */
475    {
476       int val = 0;
477       while (ec_enc_tell(&st->enc, 0) < nbCompressedBytes*8)
478       {
479          ec_enc_uint(&st->enc, val, 2);
480          val = 1-val;
481       }
482    }
483    ec_enc_done(&st->enc);
484    {
485       unsigned char *data;
486       int nbBytes = ec_byte_bytes(&st->buf);
487       if (nbBytes > nbCompressedBytes)
488       {
489          celt_warning_int ("got too many bytes:", nbBytes);
490          RESTORE_STACK;
491          return CELT_INTERNAL_ERROR;
492       }
493       /*printf ("%d\n", *nbBytes);*/
494       data = ec_byte_get_buffer(&st->buf);
495       for (i=0;i<nbBytes;i++)
496          compressed[i] = data[i];
497       for (;i<nbCompressedBytes;i++)
498          compressed[i] = 0;
499    }
500    /* Reset the packing for the next encoding */
501    ec_byte_reset(&st->buf);
502    ec_enc_init(&st->enc,&st->buf);
503
504    RESTORE_STACK;
505    return nbCompressedBytes;
506 }
507
508
509 /****************************************************************************/
510 /*                                                                          */
511 /*                                DECODER                                   */
512 /*                                                                          */
513 /****************************************************************************/
514
515
516 /** Decoder state 
517  @brief Decoder state
518  */
519 struct CELTDecoder {
520    const CELTMode *mode;
521    int frame_size;
522    int block_size;
523    int overlap;
524
525    ec_byte_buffer buf;
526    ec_enc         enc;
527
528    celt_sig_t * restrict preemph_memD;
529
530    celt_sig_t *out_mem;
531
532    celt_word16_t *oldBandE;
533    
534    int last_pitch_index;
535 };
536
537 CELTDecoder *celt_decoder_create(const CELTMode *mode)
538 {
539    int N, C;
540    CELTDecoder *st;
541
542    if (check_mode(mode) != CELT_OK)
543       return NULL;
544
545    N = mode->mdctSize;
546    C = CHANNELS(mode);
547    st = celt_alloc(sizeof(CELTDecoder));
548    
549    st->mode = mode;
550    st->frame_size = N;
551    st->block_size = N;
552    st->overlap = mode->overlap;
553
554    st->out_mem = celt_alloc((MAX_PERIOD+st->overlap)*C*sizeof(celt_sig_t));
555    
556    st->oldBandE = (celt_word16_t*)celt_alloc(C*mode->nbEBands*sizeof(celt_word16_t));
557
558    st->preemph_memD = (celt_sig_t*)celt_alloc(C*sizeof(celt_sig_t));;
559
560    st->last_pitch_index = 0;
561    return st;
562 }
563
564 void celt_decoder_destroy(CELTDecoder *st)
565 {
566    if (st == NULL)
567    {
568       celt_warning("NULL passed to celt_encoder_destroy");
569       return;
570    }
571    if (check_mode(st->mode) != CELT_OK)
572       return;
573
574
575    celt_free(st->out_mem);
576    
577    celt_free(st->oldBandE);
578    
579    celt_free(st->preemph_memD);
580
581    celt_free(st);
582 }
583
584 /** Handles lost packets by just copying past data with the same offset as the last
585     pitch period */
586 static void celt_decode_lost(CELTDecoder * restrict st, short * restrict pcm)
587 {
588    int c, N;
589    int pitch_index;
590    int i, len;
591    VARDECL(celt_sig_t, freq);
592    const int C = CHANNELS(st->mode);
593    int offset;
594    SAVE_STACK;
595    N = st->block_size;
596    ALLOC(freq,C*N, celt_sig_t);         /**< Interleaved signal MDCTs */
597    
598    len = N+st->mode->overlap;
599 #if 0
600    pitch_index = st->last_pitch_index;
601    
602    /* Use the pitch MDCT as the "guessed" signal */
603    compute_mdcts(st->mode, st->mode->window, st->out_mem+pitch_index*C, freq);
604
605 #else
606    find_spectral_pitch(st->mode, st->mode->fft, &st->mode->psy, st->out_mem+MAX_PERIOD-len, st->out_mem, st->mode->window, len, MAX_PERIOD-len-100, &pitch_index);
607    pitch_index = MAX_PERIOD-len-pitch_index;
608    offset = MAX_PERIOD-pitch_index;
609    while (offset+len >= MAX_PERIOD)
610       offset -= pitch_index;
611    compute_mdcts(st->mode, st->mode->window, st->out_mem+offset*C, freq);
612    for (i=0;i<N;i++)
613       freq[i] = MULT16_32_Q15(QCONST16(.9f,15),freq[i]);
614 #endif
615    
616    
617    
618    CELT_MOVE(st->out_mem, st->out_mem+C*N, C*(MAX_PERIOD+st->mode->overlap-N));
619    /* Compute inverse MDCTs */
620    compute_inv_mdcts(st->mode, st->mode->window, freq, -1, 1, st->out_mem);
621
622    for (c=0;c<C;c++)
623    {
624       int j;
625       for (j=0;j<N;j++)
626       {
627          celt_sig_t tmp = ADD32(st->out_mem[C*(MAX_PERIOD-N)+C*j+c],
628                                 MULT16_32_Q15(preemph,st->preemph_memD[c]));
629          st->preemph_memD[c] = tmp;
630          pcm[C*j+c] = SIG2INT16(tmp);
631       }
632    }
633    RESTORE_STACK;
634 }
635
636 int celt_decode(CELTDecoder * restrict st, unsigned char *data, int len, celt_int16_t * restrict pcm)
637 {
638    int c, N, N4;
639    int has_pitch;
640    int pitch_index;
641    ec_dec dec;
642    ec_byte_buffer buf;
643    VARDECL(celt_sig_t, freq);
644    VARDECL(celt_norm_t, X);
645    VARDECL(celt_norm_t, P);
646    VARDECL(celt_ener_t, bandE);
647    VARDECL(celt_pgain_t, gains);
648    VARDECL(int, stereo_mode);
649    int time_domain;
650    int transient_time;
651    float transient_gain;
652    const int C = CHANNELS(st->mode);
653    SAVE_STACK;
654
655    if (check_mode(st->mode) != CELT_OK)
656       return CELT_INVALID_MODE;
657
658    N = st->block_size;
659    N4 = (N-st->overlap)>>1;
660
661    ALLOC(freq, C*N, celt_sig_t); /**< Interleaved signal MDCTs */
662    ALLOC(X, C*N, celt_norm_t);         /**< Interleaved normalised MDCTs */
663    ALLOC(P, C*N, celt_norm_t);         /**< Interleaved normalised pitch MDCTs*/
664    ALLOC(bandE, st->mode->nbEBands*C, celt_ener_t);
665    ALLOC(gains, st->mode->nbPBands, celt_pgain_t);
666    
667    if (check_mode(st->mode) != CELT_OK)
668    {
669       RESTORE_STACK;
670       return CELT_INVALID_MODE;
671    }
672    if (data == NULL)
673    {
674       celt_decode_lost(st, pcm);
675       RESTORE_STACK;
676       return 0;
677    }
678    
679    ec_byte_readinit(&buf,data,len);
680    ec_dec_init(&dec,&buf);
681    
682    time_domain = ec_dec_bits(&dec, 1);
683    if (time_domain)
684    {
685       int gainid = ec_dec_bits(&dec, 2);
686       switch(gainid) {
687          case 0:
688             transient_gain = 1;
689             break;
690          case 1:
691             transient_gain = 2;
692             break;
693          case 2:
694             transient_gain = 4;
695             break;
696          case 3:
697          default:
698             transient_gain = 8;
699             break;
700       }
701       transient_time = ec_dec_uint(&dec, N+st->mode->overlap);
702    } else {
703       transient_time = -1;
704       transient_gain = 1;
705    }
706    /* Get the pitch gains */
707    has_pitch = unquant_pitch(gains, st->mode->nbPBands, &dec);
708    
709    /* Get the pitch index */
710    if (has_pitch)
711    {
712       pitch_index = ec_dec_uint(&dec, MAX_PERIOD-(2*N-2*N4));
713       st->last_pitch_index = pitch_index;
714    } else {
715       /* FIXME: We could be more intelligent here and just not compute the MDCT */
716       pitch_index = 0;
717    }
718
719    /* Get band energies */
720    unquant_energy(st->mode, bandE, st->oldBandE, 20*C+len*8/5, st->mode->prob, &dec);
721
722    /* Pitch MDCT */
723    compute_mdcts(st->mode, st->mode->window, st->out_mem+pitch_index*C, freq);
724
725    {
726       VARDECL(celt_ener_t, bandEp);
727       ALLOC(bandEp, st->mode->nbEBands*C, celt_ener_t);
728       compute_band_energies(st->mode, freq, bandEp);
729       normalise_bands(st->mode, freq, P, bandEp);
730    }
731
732    ALLOC(stereo_mode, st->mode->nbEBands, int);
733    stereo_decision(st->mode, X, stereo_mode, st->mode->nbEBands);
734    /* Apply pitch gains */
735    pitch_quant_bands(st->mode, P, gains);
736
737    /* Decode fixed codebook and merge with pitch */
738    unquant_bands(st->mode, X, P, bandE, stereo_mode, len*8, time_domain, &dec);
739
740    if (C==2)
741    {
742       renormalise_bands(st->mode, X);
743    }
744    /* Synthesis */
745    denormalise_bands(st->mode, X, freq, bandE);
746
747
748    CELT_MOVE(st->out_mem, st->out_mem+C*N, C*(MAX_PERIOD+st->overlap-N));
749    /* Compute inverse MDCTs */
750    compute_inv_mdcts(st->mode, st->mode->window, freq, transient_time, transient_gain, st->out_mem);
751
752    for (c=0;c<C;c++)
753    {
754       int j;
755       const celt_sig_t * restrict outp=st->out_mem+C*(MAX_PERIOD-N)+c;
756       celt_int16_t * restrict pcmp = pcm+c;
757       for (j=0;j<N;j++)
758       {
759          celt_sig_t tmp = ADD32(*outp, MULT16_32_Q15(preemph,st->preemph_memD[c]));
760          st->preemph_memD[c] = tmp;
761          *pcmp = SIG2INT16(tmp);
762          pcmp += C;
763          outp += C;
764       }
765    }
766
767    {
768       unsigned int val = 0;
769       while (ec_dec_tell(&dec, 0) < len*8)
770       {
771          if (ec_dec_uint(&dec, 2) != val)
772          {
773             celt_warning("decode error");
774             RESTORE_STACK;
775             return CELT_CORRUPTED_DATA;
776          }
777          val = 1-val;
778       }
779    }
780
781    RESTORE_STACK;
782    return 0;
783    /*printf ("\n");*/
784 }
785