Add a generic CDF decoding routine.
[opus.git] / libcelt / celt.c
1 /* Copyright (c) 2007-2008 CSIRO
2    Copyright (c) 2007-2010 Xiph.Org Foundation
3    Copyright (c) 2008 Gregory Maxwell 
4    Written by Jean-Marc Valin and Gregory Maxwell */
5 /*
6    Redistribution and use in source and binary forms, with or without
7    modification, are permitted provided that the following conditions
8    are met:
9    
10    - Redistributions of source code must retain the above copyright
11    notice, this list of conditions and the following disclaimer.
12    
13    - Redistributions in binary form must reproduce the above copyright
14    notice, this list of conditions and the following disclaimer in the
15    documentation and/or other materials provided with the distribution.
16    
17    - Neither the name of the Xiph.org Foundation nor the names of its
18    contributors may be used to endorse or promote products derived from
19    this software without specific prior written permission.
20    
21    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
22    ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
23    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
24    A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE FOUNDATION OR
25    CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
26    EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
27    PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
28    PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
29    LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
30    NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
31    SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
32 */
33
34 #ifdef HAVE_CONFIG_H
35 #include "config.h"
36 #endif
37
38 #define CELT_C
39
40 #include "os_support.h"
41 #include "mdct.h"
42 #include <math.h>
43 #include "celt.h"
44 #include "pitch.h"
45 #include "bands.h"
46 #include "modes.h"
47 #include "entcode.h"
48 #include "quant_bands.h"
49 #include "rate.h"
50 #include "stack_alloc.h"
51 #include "mathops.h"
52 #include "float_cast.h"
53 #include <stdarg.h>
54 #include "plc.h"
55
56 static const int trim_cdf[12] = {0, 2, 4, 9, 19, 41, 87, 109, 119, 124, 126, 128};
57 static const int spread_cdf[5] = {0, 7, 9, 30, 32};
58
59 #define COMBFILTER_MAXPERIOD 1024
60 #define COMBFILTER_MINPERIOD 16
61
62 /** Encoder state 
63  @brief Encoder state
64  */
65 struct CELTEncoder {
66    const CELTMode *mode;     /**< Mode used by the encoder */
67    int overlap;
68    int channels;
69    
70    int force_intra;
71    int complexity;
72    int start, end;
73
74    celt_int32 vbr_rate_norm; /* Target number of 8th bits per frame */
75    int constrained_vbr;      /* If zero, VBR can do whatever it likes with the rate */
76
77    /* Everything beyond this point gets cleared on a reset */
78 #define ENCODER_RESET_START frame_max
79
80    celt_word32 frame_max;
81    int spread_decision;
82    int delayedIntra;
83    int tonal_average;
84    int lastCodedBands;
85
86    int prefilter_period;
87    celt_word16 prefilter_gain;
88 #ifdef RESYNTH
89    int prefilter_period_old;
90    celt_word16 prefilter_gain_old;
91 #endif
92
93    /* VBR-related parameters */
94    celt_int32 vbr_reservoir;
95    celt_int32 vbr_drift;
96    celt_int32 vbr_offset;
97    celt_int32 vbr_count;
98
99    celt_word32 preemph_memE[2];
100    celt_word32 preemph_memD[2];
101
102 #ifdef RESYNTH
103    celt_sig syn_mem[2][2*MAX_PERIOD];
104 #endif
105
106    celt_sig in_mem[1]; /* Size = channels*mode->overlap */
107    /* celt_sig prefilter_mem[],  Size = channels*COMBFILTER_PERIOD */
108    /* celt_sig overlap_mem[],  Size = channels*mode->overlap */
109    /* celt_word16 oldEBands[], Size = channels*mode->nbEBands */
110 };
111
112 int celt_encoder_get_size(const CELTMode *mode, int channels)
113 {
114    int size = sizeof(struct CELTEncoder)
115          + (2*channels*mode->overlap-1)*sizeof(celt_sig)
116          + channels*COMBFILTER_MAXPERIOD*sizeof(celt_sig)
117          + channels*mode->nbEBands*sizeof(celt_word16);
118    return size;
119 }
120
121 CELTEncoder *celt_encoder_create(const CELTMode *mode, int channels, int *error)
122 {
123    return celt_encoder_init(
124          (CELTEncoder *)celt_alloc(celt_encoder_get_size(mode, channels)),
125          mode, channels, error);
126 }
127
128 CELTEncoder *celt_encoder_init(CELTEncoder *st, const CELTMode *mode, int channels, int *error)
129 {
130    if (channels < 0 || channels > 2)
131    {
132       if (error)
133          *error = CELT_BAD_ARG;
134       return NULL;
135    }
136
137    if (st==NULL)
138    {
139       if (error)
140          *error = CELT_ALLOC_FAIL;
141       return NULL;
142    }
143
144    CELT_MEMSET((char*)st, 0, celt_encoder_get_size(mode, channels));
145    
146    st->mode = mode;
147    st->overlap = mode->overlap;
148    st->channels = channels;
149
150    st->start = 0;
151    st->end = st->mode->effEBands;
152    st->constrained_vbr = 1;
153
154    st->vbr_rate_norm = 0;
155    st->vbr_offset = 0;
156    st->force_intra  = 0;
157    st->delayedIntra = 1;
158    st->tonal_average = 256;
159    st->spread_decision = SPREAD_NORMAL;
160    st->complexity = 5;
161
162    if (error)
163       *error = CELT_OK;
164    return st;
165 }
166
167 void celt_encoder_destroy(CELTEncoder *st)
168 {
169    celt_free(st);
170 }
171
172 static inline celt_int16 FLOAT2INT16(float x)
173 {
174    x = x*CELT_SIG_SCALE;
175    x = MAX32(x, -32768);
176    x = MIN32(x, 32767);
177    return (celt_int16)float2int(x);
178 }
179
180 static inline celt_word16 SIG2WORD16(celt_sig x)
181 {
182 #ifdef FIXED_POINT
183    x = PSHR32(x, SIG_SHIFT);
184    x = MAX32(x, -32768);
185    x = MIN32(x, 32767);
186    return EXTRACT16(x);
187 #else
188    return (celt_word16)x;
189 #endif
190 }
191
192 static int transient_analysis(const celt_word32 * restrict in, int len, int C,
193                               celt_word32 *frame_max, int overlap)
194 {
195    int i;
196    VARDECL(celt_word16, tmp);
197    celt_word32 mem0=0,mem1=0;
198    int is_transient = 0;
199    int block;
200    int N;
201    /* FIXME: Make that smaller */
202    celt_word16 bins[50];
203    SAVE_STACK;
204    ALLOC(tmp, len, celt_word16);
205
206    block = overlap/2;
207    N=len/block;
208    if (C==1)
209    {
210       for (i=0;i<len;i++)
211          tmp[i] = SHR32(in[i],SIG_SHIFT);
212    } else {
213       for (i=0;i<len;i++)
214          tmp[i] = SHR32(ADD32(in[i],in[i+len]), SIG_SHIFT+1);
215    }
216
217    /* High-pass filter: (1 - 2*z^-1 + z^-2) / (1 - z^-1 + .5*z^-2) */
218    for (i=0;i<len;i++)
219    {
220       celt_word32 x,y;
221       x = tmp[i];
222       y = ADD32(mem0, x);
223 #ifdef FIXED_POINT
224       mem0 = mem1 + y - SHL32(x,1);
225       mem1 = x - SHR32(y,1);
226 #else
227       mem0 = mem1 + y - 2*x;
228       mem1 = x - .5f*y;
229 #endif
230       tmp[i] = EXTRACT16(SHR(y,2));
231    }
232    /* First few samples are bad because we don't propagate the memory */
233    for (i=0;i<12;i++)
234       tmp[i] = 0;
235
236    for (i=0;i<N;i++)
237    {
238       int j;
239       float max_abs=0;
240       for (j=0;j<block;j++)
241          max_abs = MAX32(max_abs, tmp[i*block+j]);
242       bins[i] = max_abs;
243    }
244    for (i=0;i<N;i++)
245    {
246       int j;
247       int conseq=0;
248       celt_word16 t1, t2, t3;
249
250       t1 = MULT16_16_Q15(QCONST16(.15f, 15), bins[i]);
251       t2 = MULT16_16_Q15(QCONST16(.4f, 15), bins[i]);
252       t3 = MULT16_16_Q15(QCONST16(.15f, 15), bins[i]);
253       for (j=0;j<i;j++)
254       {
255          if (bins[j] < t1)
256             conseq++;
257          if (bins[j] < t2)
258             conseq++;
259          else
260             conseq = 0;
261       }
262       if (conseq>=3)
263          is_transient=1;
264       conseq = 0;
265       for (j=i+1;j<N;j++)
266       {
267          if (bins[j] < t3)
268             conseq++;
269          else
270             conseq = 0;
271       }
272       if (conseq>=7)
273          is_transient=1;
274    }
275    RESTORE_STACK;
276    return is_transient;
277 }
278
279 /** Apply window and compute the MDCT for all sub-frames and 
280     all channels in a frame */
281 static void compute_mdcts(const CELTMode *mode, int shortBlocks, celt_sig * restrict in, celt_sig * restrict out, int _C, int LM)
282 {
283    const int C = CHANNELS(_C);
284    if (C==1 && !shortBlocks)
285    {
286       const int overlap = OVERLAP(mode);
287       clt_mdct_forward(&mode->mdct, in, out, mode->window, overlap, mode->maxLM-LM);
288    } else {
289       const int overlap = OVERLAP(mode);
290       int N = mode->shortMdctSize<<LM;
291       int B = 1;
292       int b, c;
293       VARDECL(celt_word32, tmp);
294       SAVE_STACK;
295       if (shortBlocks)
296       {
297          /*lookup = &mode->mdct[0];*/
298          N = mode->shortMdctSize;
299          B = shortBlocks;
300       }
301       ALLOC(tmp, N, celt_word32);
302       c=0; do {
303          for (b=0;b<B;b++)
304          {
305             int j;
306             clt_mdct_forward(&mode->mdct, in+c*(B*N+overlap)+b*N, tmp, mode->window, overlap, shortBlocks ? mode->maxLM : mode->maxLM-LM);
307             /* Interleaving the sub-frames */
308             for (j=0;j<N;j++)
309                out[(j*B+b)+c*N*B] = tmp[j];
310          }
311       } while (++c<C);
312       RESTORE_STACK;
313    }
314 }
315
316 /** Compute the IMDCT and apply window for all sub-frames and 
317     all channels in a frame */
318 static void compute_inv_mdcts(const CELTMode *mode, int shortBlocks, celt_sig *X,
319       celt_sig * restrict out_mem[],
320       celt_sig * restrict overlap_mem[], int _C, int LM)
321 {
322    int c;
323    const int C = CHANNELS(_C);
324    const int N = mode->shortMdctSize<<LM;
325    const int overlap = OVERLAP(mode);
326    c=0; do {
327       int j;
328          VARDECL(celt_word32, x);
329          VARDECL(celt_word32, tmp);
330          int b;
331          int N2 = N;
332          int B = 1;
333          SAVE_STACK;
334          
335          ALLOC(x, N+overlap, celt_word32);
336          ALLOC(tmp, N, celt_word32);
337
338          if (shortBlocks)
339          {
340             N2 = mode->shortMdctSize;
341             B = shortBlocks;
342          }
343          /* Prevents problems from the imdct doing the overlap-add */
344          CELT_MEMSET(x, 0, overlap);
345
346          for (b=0;b<B;b++)
347          {
348             /* De-interleaving the sub-frames */
349             for (j=0;j<N2;j++)
350                tmp[j] = X[(j*B+b)+c*N2*B];
351             clt_mdct_backward(&mode->mdct, tmp, x+N2*b, mode->window, overlap, shortBlocks ? mode->maxLM : mode->maxLM-LM);
352          }
353
354          for (j=0;j<overlap;j++)
355             out_mem[c][j] = x[j] + overlap_mem[c][j];
356          for (;j<N;j++)
357             out_mem[c][j] = x[j];
358          for (j=0;j<overlap;j++)
359             overlap_mem[c][j] = x[N+j];
360          RESTORE_STACK;
361    } while (++c<C);
362 }
363
364 static void deemphasis(celt_sig *in[], celt_word16 *pcm, int N, int _C, const celt_word16 *coef, celt_sig *mem)
365 {
366    const int C = CHANNELS(_C);
367    int c;
368    c=0; do {
369       int j;
370       celt_sig * restrict x;
371       celt_word16  * restrict y;
372       celt_sig m = mem[c];
373       x =in[c];
374       y = pcm+c;
375       for (j=0;j<N;j++)
376       {
377          celt_sig tmp = *x + m;
378          m = MULT16_32_Q15(coef[0], tmp)
379            - MULT16_32_Q15(coef[1], *x);
380          tmp = SHL32(MULT16_32_Q15(coef[3], tmp), 2);
381          *y = SCALEOUT(SIG2WORD16(tmp));
382          x++;
383          y+=C;
384       }
385       mem[c] = m;
386    } while (++c<C);
387 }
388
389 #ifdef ENABLE_POSTFILTER
390 /* FIXME: Handle the case where T = maxperiod */
391 static void comb_filter(celt_word32 *y, celt_word32 *x, int T0, int T1, int N,
392       int C, celt_word16 g0, celt_word16 g1, const celt_word16 *window, int overlap)
393 {
394    int i;
395    /* printf ("%d %d %f %f\n", T0, T1, g0, g1); */
396    celt_word16 g00, g01, g02, g10, g11, g12;
397    celt_word16 t0, t1, t2;
398    /* zeros at theta = +/- 5*pi/6 */
399    t0 = QCONST16(.26795f, 15);
400    t1 = QCONST16(.46410f, 15);
401    t2 = QCONST16(.26795f, 15);
402    g00 = MULT16_16_Q15(g0, t0);
403    g01 = MULT16_16_Q15(g0, t1);
404    g02 = MULT16_16_Q15(g0, t2);
405    g10 = MULT16_16_Q15(g1, t0);
406    g11 = MULT16_16_Q15(g1, t1);
407    g12 = MULT16_16_Q15(g1, t2);
408    for (i=0;i<overlap;i++)
409    {
410       celt_word16 f;
411       f = MULT16_16_Q15(window[i],window[i]);
412       y[i] = x[i]
413                + MULT16_32_Q15(MULT16_16_Q15((Q15ONE-f),g01),x[i-T0])
414                + MULT16_32_Q15(MULT16_16_Q15((Q15ONE-f),g00),x[i-T0-1])
415                + MULT16_32_Q15(MULT16_16_Q15((Q15ONE-f),g02),x[i-T0+1])
416                + MULT16_32_Q15(MULT16_16_Q15(f,g11),x[i-T1])
417                + MULT16_32_Q15(MULT16_16_Q15(f,g10),x[i-T1-1])
418                + MULT16_32_Q15(MULT16_16_Q15(f,g12),x[i-T1+1]);
419
420    }
421    for (i=overlap;i<N;i++)
422       y[i] = x[i]
423                + MULT16_32_Q15(g11,x[i-T1])
424                + MULT16_32_Q15(g10,x[i-T1-1])
425                + MULT16_32_Q15(g12,x[i-T1+1]);
426 }
427 #endif /* ENABLE_POSTFILTER */
428
429 static const signed char tf_select_table[4][8] = {
430       {0, -1, 0, -1,    0,-1, 0,-1},
431       {0, -1, 0, -2,    1, 0, 1 -1},
432       {0, -2, 0, -3,    2, 0, 1 -1},
433       {0, -2, 0, -3,    2, 0, 1 -1},
434 };
435
436 static celt_word32 l1_metric(const celt_norm *tmp, int N, int LM, int width)
437 {
438    int i, j;
439    static const celt_word16 sqrtM_1[4] = {Q15ONE, QCONST16(0.70711f,15), QCONST16(0.5f,15), QCONST16(0.35355f,15)};
440    celt_word32 L1;
441    celt_word16 bias;
442    L1=0;
443    for (i=0;i<1<<LM;i++)
444    {
445       celt_word32 L2 = 0;
446       for (j=0;j<N>>LM;j++)
447          L2 = MAC16_16(L2, tmp[(j<<LM)+i], tmp[(j<<LM)+i]);
448       L1 += celt_sqrt(L2);
449    }
450    L1 = MULT16_32_Q15(sqrtM_1[LM], L1);
451    if (width==1)
452       bias = QCONST16(.12f,15)*LM;
453    else if (width==2)
454       bias = QCONST16(.05f,15)*LM;
455    else
456       bias = QCONST16(.02f,15)*LM;
457    L1 = MAC16_32_Q15(L1, bias, L1);
458    return L1;
459 }
460
461 static int tf_analysis(const CELTMode *m, celt_word16 *bandLogE, celt_word16 *oldBandE,
462       int len, int C, int isTransient, int *tf_res, int nbCompressedBytes, celt_norm *X,
463       int N0, int LM, int *tf_sum)
464 {
465    int i;
466    VARDECL(int, metric);
467    int cost0;
468    int cost1;
469    VARDECL(int, path0);
470    VARDECL(int, path1);
471    VARDECL(celt_norm, tmp);
472    int lambda;
473    int tf_select=0;
474    SAVE_STACK;
475
476    /* FIXME: Should check number of bytes *left* */
477    if (nbCompressedBytes<15*C)
478    {
479       *tf_sum = 0;
480       for (i=0;i<len;i++)
481          tf_res[i] = isTransient;
482       return 0;
483    }
484    if (nbCompressedBytes<40)
485       lambda = 12;
486    else if (nbCompressedBytes<60)
487       lambda = 6;
488    else if (nbCompressedBytes<100)
489       lambda = 4;
490    else
491       lambda = 3;
492
493    ALLOC(metric, len, int);
494    ALLOC(tmp, (m->eBands[len]-m->eBands[len-1])<<LM, celt_norm);
495    ALLOC(path0, len, int);
496    ALLOC(path1, len, int);
497
498    *tf_sum = 0;
499    for (i=0;i<len;i++)
500    {
501       int j, k, N;
502       celt_word32 L1, best_L1;
503       int best_level=0;
504       N = (m->eBands[i+1]-m->eBands[i])<<LM;
505       for (j=0;j<N;j++)
506          tmp[j] = X[j+(m->eBands[i]<<LM)];
507       /* FIXME: Do something with the right channel */
508       if (C==2)
509          for (j=0;j<N;j++)
510             tmp[j] = ADD16(tmp[j],X[N0+j+(m->eBands[i]<<LM)]);
511       L1 = l1_metric(tmp, N, isTransient ? LM : 0, N>>LM);
512       best_L1 = L1;
513       /*printf ("%f ", L1);*/
514       for (k=0;k<LM;k++)
515       {
516          int B;
517
518          if (isTransient)
519             B = (LM-k-1);
520          else
521             B = k+1;
522
523          if (isTransient)
524             haar1(tmp, N>>(LM-k), 1<<(LM-k));
525          else
526             haar1(tmp, N>>k, 1<<k);
527
528          L1 = l1_metric(tmp, N, B, N>>LM);
529
530          if (L1 < best_L1)
531          {
532             best_L1 = L1;
533             best_level = k+1;
534          }
535       }
536       /*printf ("%d ", isTransient ? LM-best_level : best_level);*/
537       if (isTransient)
538          metric[i] = best_level;
539       else
540          metric[i] = -best_level;
541       *tf_sum += metric[i];
542    }
543    /*printf("\n");*/
544    /* FIXME: Figure out how to set this */
545    tf_select = 0;
546
547    cost0 = 0;
548    cost1 = isTransient ? 0 : lambda;
549    /* Viterbi forward pass */
550    for (i=1;i<len;i++)
551    {
552       int curr0, curr1;
553       int from0, from1;
554
555       from0 = cost0;
556       from1 = cost1 + lambda;
557       if (from0 < from1)
558       {
559          curr0 = from0;
560          path0[i]= 0;
561       } else {
562          curr0 = from1;
563          path0[i]= 1;
564       }
565
566       from0 = cost0 + lambda;
567       from1 = cost1;
568       if (from0 < from1)
569       {
570          curr1 = from0;
571          path1[i]= 0;
572       } else {
573          curr1 = from1;
574          path1[i]= 1;
575       }
576       cost0 = curr0 + abs(metric[i]-tf_select_table[LM][4*isTransient+2*tf_select+0]);
577       cost1 = curr1 + abs(metric[i]-tf_select_table[LM][4*isTransient+2*tf_select+1]);
578    }
579    tf_res[len-1] = cost0 < cost1 ? 0 : 1;
580    /* Viterbi backward pass to check the decisions */
581    for (i=len-2;i>=0;i--)
582    {
583       if (tf_res[i+1] == 1)
584          tf_res[i] = path1[i+1];
585       else
586          tf_res[i] = path0[i+1];
587    }
588    RESTORE_STACK;
589    return tf_select;
590 }
591
592 static void tf_encode(int start, int end, int isTransient, int *tf_res, int LM, int tf_select, ec_enc *enc)
593 {
594    int curr, i;
595    ec_enc_bit_prob(enc, tf_res[start], isTransient ? 16384 : 4096);
596    curr = tf_res[start];
597    for (i=start+1;i<end;i++)
598    {
599       ec_enc_bit_prob(enc, tf_res[i] ^ curr, isTransient ? 4096 : 2048);
600       curr = tf_res[i];
601    }
602    if (LM!=0)
603       ec_enc_bits(enc, tf_select, 1);
604    for (i=start;i<end;i++)
605       tf_res[i] = tf_select_table[LM][4*isTransient+2*tf_select+tf_res[i]];
606    /*printf("%d %d ", isTransient, tf_select); for(i=0;i<end;i++)printf("%d ", tf_res[i]);printf("\n");*/
607 }
608
609 static void tf_decode(int start, int end, int C, int isTransient, int *tf_res, int LM, ec_dec *dec)
610 {
611    int i, curr, tf_select;
612    tf_res[start] = ec_dec_bit_prob(dec, isTransient ? 16384 : 4096);
613    curr = tf_res[start];
614    for (i=start+1;i<end;i++)
615    {
616       tf_res[i] = ec_dec_bit_prob(dec, isTransient ? 4096 : 2048) ^ curr;
617       curr = tf_res[i];
618    }
619    if (LM!=0)
620       tf_select = ec_dec_bits(dec, 1);
621    else
622       tf_select = 0;
623    for (i=start;i<end;i++)
624       tf_res[i] = tf_select_table[LM][4*isTransient+2*tf_select+tf_res[i]];
625 }
626
627 static int alloc_trim_analysis(const CELTMode *m, const celt_norm *X,
628       const celt_word16 *bandLogE, int nbEBands, int LM, int C, int N0)
629 {
630    int i;
631    celt_word32 diff=0;
632    int c;
633    int trim_index = 5;
634    if (C==2)
635    {
636       celt_word16 sum = 0; /* Q10 */
637       /* Compute inter-channel correlation for low frequencies */
638       for (i=0;i<8;i++)
639       {
640          int j;
641          celt_word32 partial = 0;
642          for (j=m->eBands[i]<<LM;j<m->eBands[i+1]<<LM;j++)
643             partial = MAC16_16(partial, X[j], X[N0+j]);
644          sum = ADD16(sum, EXTRACT16(SHR32(partial, 18)));
645       }
646       sum = MULT16_16_Q15(QCONST16(1.f/8, 15), sum);
647       /*printf ("%f\n", sum);*/
648       if (sum > QCONST16(.995f,10))
649          trim_index-=4;
650       else if (sum > QCONST16(.92f,10))
651          trim_index-=3;
652       else if (sum > QCONST16(.85f,10))
653          trim_index-=2;
654       else if (sum > QCONST16(.8f,10))
655          trim_index-=1;
656    }
657
658    /* Estimate spectral tilt */
659    c=0; do {
660       for (i=0;i<nbEBands-1;i++)
661       {
662          diff += bandLogE[i+c*nbEBands]*(celt_int32)(2+2*i-nbEBands);
663       }
664    } while (++c<0);
665    diff /= C*(nbEBands-1);
666    /*printf("%f\n", diff);*/
667    if (diff > QCONST16(2.f, DB_SHIFT))
668       trim_index--;
669    if (diff > QCONST16(8.f, DB_SHIFT))
670       trim_index--;
671    if (diff < -QCONST16(4.f, DB_SHIFT))
672       trim_index++;
673    if (diff < -QCONST16(10.f, DB_SHIFT))
674       trim_index++;
675
676    if (trim_index<0)
677       trim_index = 0;
678    if (trim_index>10)
679       trim_index = 10;
680    return trim_index;
681 }
682
683 static int stereo_analysis(const CELTMode *m, const celt_norm *X,
684       int nbEBands, int LM, int C, int N0)
685 {
686    int i;
687    int thetas;
688    celt_word32 sumLR = EPSILON, sumMS = EPSILON;
689
690    /* Use the L1 norm to model the entropy of the L/R signal vs the M/S signal */
691    for (i=0;i<13;i++)
692    {
693       int j;
694       for (j=m->eBands[i]<<LM;j<m->eBands[i+1]<<LM;j++)
695       {
696          celt_word16 L, R, M, S;
697          L = X[j];
698          R = X[N0+j];
699          M = L+R;
700          S = L-R;
701          sumLR += EXTEND32(ABS16(L)) + EXTEND32(ABS16(R));
702          sumMS += EXTEND32(ABS16(M)) + EXTEND32(ABS16(S));
703       }
704    }
705    sumMS = MULT16_32_Q15(QCONST16(0.707107f, 15), sumMS);
706    thetas = 13;
707    /* We don't need thetas for lower bands with LM<=1 */
708    if (LM<=1)
709       thetas -= 8;
710    return MULT16_32_Q15((m->eBands[13]<<(LM+1))+thetas, sumMS)
711          > MULT16_32_Q15(m->eBands[13]<<(LM+1), sumLR);
712 }
713
714 #ifdef FIXED_POINT
715 int celt_encode_with_ec(CELTEncoder * restrict st, const celt_int16 * pcm, int frame_size, unsigned char *compressed, int nbCompressedBytes, ec_enc *enc)
716 {
717 #else
718 int celt_encode_with_ec_float(CELTEncoder * restrict st, const celt_sig * pcm, int frame_size, unsigned char *compressed, int nbCompressedBytes, ec_enc *enc)
719 {
720 #endif
721    int i, c, N;
722    int bits;
723    ec_byte_buffer buf;
724    ec_enc         _enc;
725    VARDECL(celt_sig, in);
726    VARDECL(celt_sig, freq);
727    VARDECL(celt_norm, X);
728    VARDECL(celt_ener, bandE);
729    VARDECL(celt_word16, bandLogE);
730    VARDECL(int, fine_quant);
731    VARDECL(celt_word16, error);
732    VARDECL(int, pulses);
733    VARDECL(int, offsets);
734    VARDECL(int, fine_priority);
735    VARDECL(int, tf_res);
736    celt_sig *_overlap_mem;
737    celt_sig *prefilter_mem;
738    celt_word16 *oldBandE;
739    int shortBlocks=0;
740    int isTransient=0;
741    int resynth;
742    const int C = CHANNELS(st->channels);
743    int LM, M;
744    int tf_select;
745    int nbFilledBytes, nbAvailableBytes;
746    int effEnd;
747    int codedBands;
748    int tf_sum;
749    int alloc_trim;
750    int pitch_index=COMBFILTER_MINPERIOD;
751    celt_word16 gain1 = 0;
752    int intensity=0;
753    int dual_stereo=0;
754    int effectiveBytes;
755    SAVE_STACK;
756
757    if (nbCompressedBytes<0 || pcm==NULL)
758      return CELT_BAD_ARG;
759
760    for (LM=0;LM<4;LM++)
761       if (st->mode->shortMdctSize<<LM==frame_size)
762          break;
763    if (LM>=MAX_CONFIG_SIZES)
764       return CELT_BAD_ARG;
765    M=1<<LM;
766
767    prefilter_mem = st->in_mem+C*(st->overlap);
768    _overlap_mem = prefilter_mem+C*COMBFILTER_MAXPERIOD;
769    /*_overlap_mem = st->in_mem+C*(st->overlap);*/
770    oldBandE = (celt_word16*)(st->in_mem+C*(2*st->overlap+COMBFILTER_MAXPERIOD));
771
772    if (enc==NULL)
773    {
774       ec_byte_writeinit_buffer(&buf, compressed, nbCompressedBytes);
775       ec_enc_init(&_enc,&buf);
776       enc = &_enc;
777       nbFilledBytes=0;
778    } else {
779       nbFilledBytes=(ec_enc_tell(enc, 0)+4)>>3;
780    }
781    nbAvailableBytes = nbCompressedBytes - nbFilledBytes;
782
783    if (st->vbr_rate_norm>0)
784       effectiveBytes = st->vbr_rate_norm>>BITRES<<LM>>3;
785    else
786       effectiveBytes = nbCompressedBytes;
787
788    effEnd = st->end;
789    if (effEnd > st->mode->effEBands)
790       effEnd = st->mode->effEBands;
791
792    N = M*st->mode->shortMdctSize;
793    ALLOC(in, C*(N+st->overlap), celt_sig);
794
795    /* Find pitch period and gain */
796    {
797       VARDECL(celt_sig, _pre);
798       celt_sig *pre[2];
799       SAVE_STACK;
800       c = 0;
801       ALLOC(_pre, C*(N+COMBFILTER_MAXPERIOD), celt_sig);
802
803       pre[0] = _pre;
804       pre[1] = _pre + (N+COMBFILTER_MAXPERIOD);
805
806       c=0; do {
807          const celt_word16 * restrict pcmp = pcm+c;
808          celt_sig * restrict inp = in+c*(N+st->overlap)+st->overlap;
809
810          for (i=0;i<N;i++)
811          {
812             /* Apply pre-emphasis */
813             celt_sig tmp = MULT16_16(st->mode->preemph[2], SCALEIN(*pcmp));
814             *inp = tmp + st->preemph_memE[c];
815             st->preemph_memE[c] = MULT16_32_Q15(st->mode->preemph[1], *inp)
816                                    - MULT16_32_Q15(st->mode->preemph[0], tmp);
817             inp++;
818             pcmp+=C;
819          }
820          CELT_COPY(pre[c], prefilter_mem+c*COMBFILTER_MAXPERIOD, COMBFILTER_MAXPERIOD);
821          CELT_COPY(pre[c]+COMBFILTER_MAXPERIOD, in+c*(N+st->overlap)+st->overlap, N);
822       } while (++c<C);
823
824 #ifdef ENABLE_POSTFILTER
825       if (nbAvailableBytes>12*C)
826       {
827          VARDECL(celt_word16, pitch_buf);
828          ALLOC(pitch_buf, (COMBFILTER_MAXPERIOD+N)>>1, celt_word16);
829          celt_word32 tmp=0;
830          celt_word32 mem0[2]={0,0};
831          celt_word16 mem1[2]={0,0};
832
833          pitch_downsample(pre, pitch_buf, COMBFILTER_MAXPERIOD+N, COMBFILTER_MAXPERIOD+N,
834                           C, mem0, mem1);
835          pitch_search(st->mode, pitch_buf+(COMBFILTER_MAXPERIOD>>1), pitch_buf, N,
836                COMBFILTER_MAXPERIOD-COMBFILTER_MINPERIOD, &pitch_index, &tmp, 1<<LM);
837          pitch_index = COMBFILTER_MAXPERIOD-pitch_index;
838
839          gain1 = remove_doubling(pitch_buf, COMBFILTER_MAXPERIOD, COMBFILTER_MINPERIOD,
840                N, &pitch_index, st->prefilter_period, st->prefilter_gain);
841          if (pitch_index > COMBFILTER_MAXPERIOD)
842             pitch_index = COMBFILTER_MAXPERIOD;
843          gain1 = MULT16_16_Q15(QCONST16(.7f,15),gain1);
844          if (gain1 > QCONST16(.6f,15))
845             gain1 = QCONST16(.6f,15);
846          if (ABS16(gain1-st->prefilter_gain)<QCONST16(.1,15))
847             gain1=st->prefilter_gain;
848       } else {
849          gain1 = 0;
850       }
851       if (gain1<QCONST16(.2f,15) || (nbAvailableBytes<30 && gain1<QCONST16(.4f,15)))
852       {
853          ec_enc_bit_prob(enc, 0, 32768);
854          gain1 = 0;
855       } else {
856          int qg;
857          int octave;
858 #ifdef FIXED_POINT
859          qg = ((gain1+2048)>>12)-2;
860 #else
861          qg = floor(.5+gain1*8)-2;
862 #endif
863          ec_enc_bit_prob(enc, 1, 32768);
864          octave = EC_ILOG(pitch_index)-5;
865          ec_enc_uint(enc, octave, 6);
866          ec_enc_bits(enc, pitch_index-(16<<octave), 4+octave);
867          ec_enc_bits(enc, qg, 2);
868          gain1 = QCONST16(.125f,15)*(qg+2);
869       }
870       /*printf("%d %f\n", pitch_index, gain1);*/
871 #else /* ENABLE_POSTFILTER */
872       ec_enc_bit_prob(enc, 0, 32768);
873 #endif /* ENABLE_POSTFILTER */
874
875       c=0; do {
876          st->prefilter_period=IMAX(st->prefilter_period, COMBFILTER_MINPERIOD);
877          CELT_COPY(in+c*(N+st->overlap), st->in_mem+c*(st->overlap), st->overlap);
878 #ifdef ENABLE_POSTFILTER
879          comb_filter(in+c*(N+st->overlap)+st->overlap, pre[c]+COMBFILTER_MAXPERIOD,
880                st->prefilter_period, pitch_index, N, C, -st->prefilter_gain, -gain1, st->mode->window, st->mode->overlap);
881 #endif /* ENABLE_POSTFILTER */
882          CELT_COPY(st->in_mem+c*(st->overlap), in+c*(N+st->overlap)+N, st->overlap);
883
884 #ifdef ENABLE_POSTFILTER
885          if (N>COMBFILTER_MAXPERIOD)
886          {
887             CELT_MOVE(prefilter_mem+c*COMBFILTER_MAXPERIOD, pre[c]+N, COMBFILTER_MAXPERIOD);
888          } else {
889             CELT_MOVE(prefilter_mem+c*COMBFILTER_MAXPERIOD, prefilter_mem+c*COMBFILTER_MAXPERIOD+N, COMBFILTER_MAXPERIOD-N);
890             CELT_MOVE(prefilter_mem+c*COMBFILTER_MAXPERIOD+COMBFILTER_MAXPERIOD-N, pre[c]+COMBFILTER_MAXPERIOD, N);
891          }
892 #endif /* ENABLE_POSTFILTER */
893       } while (++c<C);
894
895       RESTORE_STACK;
896    }
897
898 #ifdef RESYNTH
899    resynth = 1;
900 #else
901    resynth = 0;
902 #endif
903
904    if (st->complexity > 1 && LM>0)
905    {
906       isTransient = M > 1 &&
907          transient_analysis(in, N+st->overlap, C, &st->frame_max, st->overlap);
908    } else {
909       isTransient = 0;
910    }
911
912    if (isTransient)
913       shortBlocks = M;
914    else
915       shortBlocks = 0;
916
917    ALLOC(freq, C*N, celt_sig); /**< Interleaved signal MDCTs */
918    ALLOC(bandE,st->mode->nbEBands*C, celt_ener);
919    ALLOC(bandLogE,st->mode->nbEBands*C, celt_word16);
920    /* Compute MDCTs */
921    compute_mdcts(st->mode, shortBlocks, in, freq, C, LM);
922
923    ALLOC(X, C*N, celt_norm);         /**< Interleaved normalised MDCTs */
924
925    compute_band_energies(st->mode, freq, bandE, effEnd, C, M);
926
927    amp2Log2(st->mode, effEnd, st->end, bandE, bandLogE, C);
928
929    /* Band normalisation */
930    normalise_bands(st->mode, freq, X, bandE, effEnd, C, M);
931
932    ALLOC(tf_res, st->mode->nbEBands, int);
933    /* Needs to be before coarse energy quantization because otherwise the energy gets modified */
934    tf_select = tf_analysis(st->mode, bandLogE, oldBandE, effEnd, C, isTransient, tf_res, effectiveBytes, X, N, LM, &tf_sum);
935    for (i=effEnd;i<st->end;i++)
936       tf_res[i] = tf_res[effEnd-1];
937
938    ALLOC(error, C*st->mode->nbEBands, celt_word16);
939    quant_coarse_energy(st->mode, st->start, st->end, effEnd, bandLogE,
940          oldBandE, nbCompressedBytes*8, error, enc,
941          C, LM, nbAvailableBytes, st->force_intra,
942          &st->delayedIntra, st->complexity >= 4);
943
944    if (LM > 0)
945       ec_enc_bit_prob(enc, shortBlocks!=0, 8192);
946
947    tf_encode(st->start, st->end, isTransient, tf_res, LM, tf_select, enc);
948
949    if (shortBlocks || st->complexity < 3 || nbAvailableBytes < 10*C)
950    {
951       if (st->complexity == 0)
952       {
953          st->spread_decision = SPREAD_NONE;
954       } else {
955          st->spread_decision = SPREAD_NORMAL;
956       }
957    } else {
958       st->spread_decision = spreading_decision(st->mode, X, &st->tonal_average, st->spread_decision, effEnd, C, M);
959    }
960    /* Probs: NONE: 21.875%, LIGHT: 6.25%, NORMAL: 65.625%, AGGRESSIVE: 6.25% */
961    ec_encode_bin(enc, spread_cdf[st->spread_decision],
962          spread_cdf[st->spread_decision+1], 5);
963
964    ALLOC(offsets, st->mode->nbEBands, int);
965
966    for (i=0;i<st->mode->nbEBands;i++)
967       offsets[i] = 0;
968    /* Dynamic allocation code */
969    /* Make sure that dynamic allocation can't make us bust the budget */
970    if (effectiveBytes > 50 && LM>=1)
971    {
972       int t1, t2;
973       if (LM <= 1)
974       {
975          t1 = 3;
976          t2 = 5;
977       } else {
978          t1 = 2;
979          t2 = 4;
980       }
981       for (i=1;i<st->mode->nbEBands-1;i++)
982       {
983          celt_word32 d2;
984          d2 = 2*bandLogE[i]-bandLogE[i-1]-bandLogE[i+1];
985          if (C==2)
986             d2 = HALF32(d2 + 2*bandLogE[i+st->mode->nbEBands]-
987                   bandLogE[i-1+st->mode->nbEBands]-bandLogE[i+1+st->mode->nbEBands]);
988          if (d2 > SHL16(t1,DB_SHIFT))
989             offsets[i] += 1;
990          if (d2 > SHL16(t2,DB_SHIFT))
991             offsets[i] += 1;
992       }
993    }
994    for (i=0;i<st->mode->nbEBands;i++)
995    {
996       int j;
997       ec_enc_bit_prob(enc, offsets[i]!=0, 1024);
998       if (offsets[i]!=0)
999       {
1000          for (j=0;j<offsets[i]-1;j++)
1001             ec_enc_bit_prob(enc, 1, 32768);
1002          ec_enc_bit_prob(enc, 0, 32768);
1003       }
1004       offsets[i] *= (6<<BITRES);
1005    }
1006    alloc_trim = alloc_trim_analysis(st->mode, X, bandLogE, st->mode->nbEBands, LM, C, N);
1007    ec_encode_bin(enc, trim_cdf[alloc_trim], trim_cdf[alloc_trim+1], 7);
1008
1009    /* Variable bitrate */
1010    if (st->vbr_rate_norm>0)
1011    {
1012      celt_word16 alpha;
1013      celt_int32 delta, tell;
1014      /* The target rate in 8th bits per frame */
1015      celt_int32 vbr_rate;
1016      celt_int32 target;
1017      celt_int32 vbr_bound, max_allowed, min_allowed;
1018
1019      target = vbr_rate = M*st->vbr_rate_norm;
1020
1021      target = target + st->vbr_offset - ((40*C+20)<<BITRES);
1022
1023      /* Shortblocks get a large boost in bitrate, but since they
1024         are uncommon long blocks are not greatly affected */
1025      if (shortBlocks || tf_sum < -2*(st->end-st->start))
1026         target = 7*target/4;
1027      else if (tf_sum < -(st->end-st->start))
1028         target = 3*target/2;
1029      else if (M > 1)
1030         target-=(target+14)/28;
1031
1032      tell = ec_enc_tell(enc, BITRES);
1033
1034      /* The current offset is removed from the target and the space used
1035         so far is added*/
1036      target=target+tell;
1037      /* By how much did we "miss" the target on that frame */
1038      delta = target - vbr_rate;
1039
1040      /* Computes the max bit-rate allowed in VBR more to avoid violating the target rate and buffering */
1041      vbr_bound = vbr_rate;
1042      if (st->constrained_vbr)
1043         max_allowed = IMIN(vbr_rate+vbr_bound-st->vbr_reservoir>>(BITRES+3),nbAvailableBytes);
1044      else
1045         max_allowed = nbAvailableBytes;
1046      min_allowed = (tell>>(BITRES+3)) + 2 - nbFilledBytes;
1047
1048      /* In VBR mode the frame size must not be reduced so much that it would result in the encoder running out of bits */
1049      nbAvailableBytes = target+(1<<(BITRES+2))>>(BITRES+3);
1050      nbAvailableBytes=IMAX(min_allowed,IMIN(max_allowed,nbAvailableBytes));
1051      target=nbAvailableBytes<<(BITRES+3);
1052
1053      if (st->vbr_count < 970)
1054      {
1055         st->vbr_count++;
1056         alpha = celt_rcp(SHL32(EXTEND32(st->vbr_count+20),16));
1057      } else
1058         alpha = QCONST16(.001f,15);
1059      /* How many bits have we used in excess of what we're allowed */
1060      if (st->constrained_vbr)
1061         st->vbr_reservoir += target - vbr_rate;
1062      /*printf ("%d\n", st->vbr_reservoir);*/
1063
1064      /* Compute the offset we need to apply in order to reach the target */
1065      st->vbr_drift += (celt_int32)MULT16_32_Q15(alpha,delta-st->vbr_offset-st->vbr_drift);
1066      st->vbr_offset = -st->vbr_drift;
1067      /*printf ("%d\n", st->vbr_drift);*/
1068
1069      /* We could use any multiple of vbr_rate as bound (depending on the delay) */
1070      if (st->constrained_vbr && st->vbr_reservoir < 0)
1071      {
1072         /* We're under the min value -- increase rate */
1073         int adjust = (-st->vbr_reservoir)/(8<<BITRES);
1074         nbAvailableBytes += adjust;
1075         st->vbr_reservoir = 0;
1076         /*printf ("+%d\n", adjust);*/
1077      }
1078      nbCompressedBytes = IMIN(nbCompressedBytes,nbAvailableBytes+nbFilledBytes);
1079
1080      /* This moves the raw bits to take into account the new compressed size */
1081      ec_byte_shrink(&buf, nbCompressedBytes);
1082    }
1083
1084    if (C==2)
1085    {
1086       /* Always use MS for 2.5 ms frames until we can do a better analysis */
1087       if (LM==0)
1088          dual_stereo = 0;
1089       else
1090          dual_stereo = stereo_analysis(st->mode, X, st->mode->nbEBands, LM, C, N);
1091       ec_enc_bit_prob(enc, dual_stereo, 32768);
1092    }
1093    if (C==2)
1094    {
1095       int effectiveRate;
1096
1097       /* Account for coarse energy */
1098       effectiveRate = (8*effectiveBytes - 80)>>LM;
1099
1100       /* effectiveRate in kb/s */
1101       effectiveRate = 2*effectiveRate/5;
1102       if (effectiveRate<35)
1103          intensity = 8;
1104       else if (effectiveRate<50)
1105          intensity = 12;
1106       else if (effectiveRate<68)
1107          intensity = 16;
1108       else if (effectiveRate<84)
1109          intensity = 18;
1110       else if (effectiveRate<102)
1111          intensity = 19;
1112       else if (effectiveRate<130)
1113          intensity = 20;
1114       else
1115          intensity = 100;
1116       intensity = IMIN(st->end,IMAX(st->start, intensity));
1117       ec_enc_uint(enc, intensity, 1+st->end-st->start);
1118    }
1119
1120    /* Bit allocation */
1121    ALLOC(fine_quant, st->mode->nbEBands, int);
1122    ALLOC(pulses, st->mode->nbEBands, int);
1123    ALLOC(fine_priority, st->mode->nbEBands, int);
1124
1125    /* bits =   packet size        -       where we are           - safety */
1126    bits = (nbCompressedBytes*8<<BITRES) - ec_enc_tell(enc, BITRES) - 1;
1127    codedBands = compute_allocation(st->mode, st->start, st->end, offsets,
1128          alloc_trim, bits, pulses, fine_quant, fine_priority, C, LM, enc, 1, st->lastCodedBands);
1129    st->lastCodedBands = codedBands;
1130
1131    quant_fine_energy(st->mode, st->start, st->end, bandE, oldBandE, error, fine_quant, enc, C);
1132
1133 #ifdef MEASURE_NORM_MSE
1134    float X0[3000];
1135    float bandE0[60];
1136    c=0; do 
1137       for (i=0;i<N;i++)
1138          X0[i+c*N] = X[i+c*N];
1139    while (++c<C);
1140    for (i=0;i<C*st->mode->nbEBands;i++)
1141       bandE0[i] = bandE[i];
1142 #endif
1143
1144    /* Residual quantisation */
1145    quant_all_bands(1, st->mode, st->start, st->end, X, C==2 ? X+N : NULL,
1146          bandE, pulses, shortBlocks, st->spread_decision, dual_stereo, intensity, tf_res, resynth,
1147          nbCompressedBytes*8, enc, LM, codedBands);
1148
1149    quant_energy_finalise(st->mode, st->start, st->end, bandE, oldBandE, error, fine_quant, fine_priority, nbCompressedBytes*8-ec_enc_tell(enc, 0), enc, C);
1150
1151 #ifdef RESYNTH
1152    /* Re-synthesis of the coded audio if required */
1153    if (resynth)
1154    {
1155       celt_sig *out_mem[2];
1156       celt_sig *overlap_mem[2];
1157
1158       log2Amp(st->mode, st->start, st->end, bandE, oldBandE, C);
1159
1160 #ifdef MEASURE_NORM_MSE
1161       measure_norm_mse(st->mode, X, X0, bandE, bandE0, M, N, C);
1162 #endif
1163
1164       /* Synthesis */
1165       denormalise_bands(st->mode, X, freq, bandE, effEnd, C, M);
1166
1167       CELT_MOVE(st->syn_mem[0], st->syn_mem[0]+N, MAX_PERIOD);
1168       if (C==2)
1169          CELT_MOVE(st->syn_mem[1], st->syn_mem[1]+N, MAX_PERIOD);
1170
1171       c=0; do
1172          for (i=0;i<M*st->mode->eBands[st->start];i++)
1173             freq[c*N+i] = 0;
1174       while (++c<C);
1175       c=0; do
1176          for (i=M*st->mode->eBands[st->end];i<N;i++)
1177             freq[c*N+i] = 0;
1178       while (++c<C);
1179
1180       out_mem[0] = st->syn_mem[0]+MAX_PERIOD;
1181       if (C==2)
1182          out_mem[1] = st->syn_mem[1]+MAX_PERIOD;
1183
1184       c=0; do
1185          overlap_mem[c] = _overlap_mem + c*st->overlap;
1186       while (++c<C);
1187
1188       compute_inv_mdcts(st->mode, shortBlocks, freq, out_mem, overlap_mem, C, LM);
1189
1190 #ifdef ENABLE_POSTFILTER
1191       c=0; do {
1192          st->prefilter_period=IMAX(st->prefilter_period, COMBFILTER_MINPERIOD);
1193          st->prefilter_period_old=IMAX(st->prefilter_period_old, COMBFILTER_MINPERIOD);
1194          if (LM!=0)
1195          {
1196             comb_filter(out_mem[c], out_mem[c], st->prefilter_period, st->prefilter_period, st->overlap, C,
1197                   st->prefilter_gain, st->prefilter_gain, NULL, 0);
1198             comb_filter(out_mem[c]+st->overlap, out_mem[c]+st->overlap, st->prefilter_period, pitch_index, N-st->overlap, C,
1199                   st->prefilter_gain, gain1, st->mode->window, st->mode->overlap);
1200          } else {
1201             comb_filter(out_mem[c], out_mem[c], st->prefilter_period_old, st->prefilter_period, N, C,
1202                   st->prefilter_gain_old, st->prefilter_gain, st->mode->window, st->mode->overlap);
1203          }
1204       } while (++c<C);
1205 #endif /* ENABLE_POSTFILTER */
1206
1207       deemphasis(out_mem, (celt_word16*)pcm, N, C, st->mode->preemph, st->preemph_memD);
1208       st->prefilter_period_old = st->prefilter_period;
1209       st->prefilter_gain_old = st->prefilter_gain;
1210    }
1211 #endif
1212
1213    st->prefilter_period = pitch_index;
1214    st->prefilter_gain = gain1;
1215
1216    /* If there's any room left (can only happen for very high rates),
1217       fill it with zeros */
1218    while (ec_enc_tell(enc,0) + 8 <= nbCompressedBytes*8)
1219       ec_enc_bits(enc, 0, 8);
1220    ec_enc_done(enc);
1221    
1222    RESTORE_STACK;
1223    if (ec_enc_get_error(enc))
1224       return CELT_CORRUPTED_DATA;
1225    else
1226       return nbCompressedBytes;
1227 }
1228
1229 #ifdef FIXED_POINT
1230 #ifndef DISABLE_FLOAT_API
1231 int celt_encode_with_ec_float(CELTEncoder * restrict st, const float * pcm, int frame_size, unsigned char *compressed, int nbCompressedBytes, ec_enc *enc)
1232 {
1233    int j, ret, C, N, LM, M;
1234    VARDECL(celt_int16, in);
1235    SAVE_STACK;
1236
1237    if (pcm==NULL)
1238       return CELT_BAD_ARG;
1239
1240    for (LM=0;LM<4;LM++)
1241       if (st->mode->shortMdctSize<<LM==frame_size)
1242          break;
1243    if (LM>=MAX_CONFIG_SIZES)
1244       return CELT_BAD_ARG;
1245    M=1<<LM;
1246
1247    C = CHANNELS(st->channels);
1248    N = M*st->mode->shortMdctSize;
1249    ALLOC(in, C*N, celt_int16);
1250
1251    for (j=0;j<C*N;j++)
1252      in[j] = FLOAT2INT16(pcm[j]);
1253
1254    ret=celt_encode_with_ec(st,in,frame_size,compressed,nbCompressedBytes, enc);
1255 #ifdef RESYNTH
1256    for (j=0;j<C*N;j++)
1257       ((float*)pcm)[j]=in[j]*(1.f/32768.f);
1258 #endif
1259    RESTORE_STACK;
1260    return ret;
1261
1262 }
1263 #endif /*DISABLE_FLOAT_API*/
1264 #else
1265 int celt_encode_with_ec(CELTEncoder * restrict st, const celt_int16 * pcm, int frame_size, unsigned char *compressed, int nbCompressedBytes, ec_enc *enc)
1266 {
1267    int j, ret, C, N, LM, M;
1268    VARDECL(celt_sig, in);
1269    SAVE_STACK;
1270
1271    if (pcm==NULL)
1272       return CELT_BAD_ARG;
1273
1274    for (LM=0;LM<4;LM++)
1275       if (st->mode->shortMdctSize<<LM==frame_size)
1276          break;
1277    if (LM>=MAX_CONFIG_SIZES)
1278       return CELT_BAD_ARG;
1279    M=1<<LM;
1280
1281    C=CHANNELS(st->channels);
1282    N=M*st->mode->shortMdctSize;
1283    ALLOC(in, C*N, celt_sig);
1284    for (j=0;j<C*N;j++) {
1285      in[j] = SCALEOUT(pcm[j]);
1286    }
1287
1288    ret = celt_encode_with_ec_float(st,in,frame_size,compressed,nbCompressedBytes, enc);
1289 #ifdef RESYNTH
1290    for (j=0;j<C*N;j++)
1291       ((celt_int16*)pcm)[j] = FLOAT2INT16(in[j]);
1292 #endif
1293    RESTORE_STACK;
1294    return ret;
1295 }
1296 #endif
1297
1298 int celt_encode(CELTEncoder * restrict st, const celt_int16 * pcm, int frame_size, unsigned char *compressed, int nbCompressedBytes)
1299 {
1300    return celt_encode_with_ec(st, pcm, frame_size, compressed, nbCompressedBytes, NULL);
1301 }
1302
1303 #ifndef DISABLE_FLOAT_API
1304 int celt_encode_float(CELTEncoder * restrict st, const float * pcm, int frame_size, unsigned char *compressed, int nbCompressedBytes)
1305 {
1306    return celt_encode_with_ec_float(st, pcm, frame_size, compressed, nbCompressedBytes, NULL);
1307 }
1308 #endif /* DISABLE_FLOAT_API */
1309
1310 int celt_encoder_ctl(CELTEncoder * restrict st, int request, ...)
1311 {
1312    va_list ap;
1313    
1314    va_start(ap, request);
1315    switch (request)
1316    {
1317       case CELT_GET_MODE_REQUEST:
1318       {
1319          const CELTMode ** value = va_arg(ap, const CELTMode**);
1320          if (value==0)
1321             goto bad_arg;
1322          *value=st->mode;
1323       }
1324       break;
1325       case CELT_SET_COMPLEXITY_REQUEST:
1326       {
1327          int value = va_arg(ap, celt_int32);
1328          if (value<0 || value>10)
1329             goto bad_arg;
1330          st->complexity = value;
1331       }
1332       break;
1333       case CELT_SET_START_BAND_REQUEST:
1334       {
1335          celt_int32 value = va_arg(ap, celt_int32);
1336          if (value<0 || value>=st->mode->nbEBands)
1337             goto bad_arg;
1338          st->start = value;
1339       }
1340       break;
1341       case CELT_SET_END_BAND_REQUEST:
1342       {
1343          celt_int32 value = va_arg(ap, celt_int32);
1344          if (value<0 || value>=st->mode->nbEBands)
1345             goto bad_arg;
1346          st->end = value;
1347       }
1348       break;
1349       case CELT_SET_PREDICTION_REQUEST:
1350       {
1351          int value = va_arg(ap, celt_int32);
1352          if (value<0 || value>2)
1353             goto bad_arg;
1354          if (value==0)
1355          {
1356             st->force_intra   = 1;
1357          } else if (value==1) {
1358             st->force_intra   = 0;
1359          } else {
1360             st->force_intra   = 0;
1361          }   
1362       }
1363       break;
1364       case CELT_SET_VBR_CONSTRAINT_REQUEST:
1365       {
1366          celt_int32 value = va_arg(ap, celt_int32);
1367          st->constrained_vbr = value;
1368       }
1369       break;
1370       case CELT_SET_VBR_RATE_REQUEST:
1371       {
1372          celt_int32 value = va_arg(ap, celt_int32);
1373          int frame_rate;
1374          int N = st->mode->shortMdctSize;
1375          if (value<0)
1376             goto bad_arg;
1377          if (value>3072000)
1378             value = 3072000;
1379          frame_rate = ((st->mode->Fs<<3)+(N>>1))/N;
1380          st->vbr_rate_norm = ((value<<(BITRES+3))+(frame_rate>>1))/frame_rate;
1381       }
1382       break;
1383       case CELT_RESET_STATE:
1384       {
1385          CELT_MEMSET((char*)&st->ENCODER_RESET_START, 0,
1386                celt_encoder_get_size(st->mode, st->channels)-
1387                ((char*)&st->ENCODER_RESET_START - (char*)st));
1388          st->vbr_offset = 0;
1389          st->delayedIntra = 1;
1390          st->spread_decision = SPREAD_NORMAL;
1391          st->tonal_average = QCONST16(1.f,8);
1392       }
1393       break;
1394       default:
1395          goto bad_request;
1396    }
1397    va_end(ap);
1398    return CELT_OK;
1399 bad_arg:
1400    va_end(ap);
1401    return CELT_BAD_ARG;
1402 bad_request:
1403    va_end(ap);
1404    return CELT_UNIMPLEMENTED;
1405 }
1406
1407 /**********************************************************************/
1408 /*                                                                    */
1409 /*                             DECODER                                */
1410 /*                                                                    */
1411 /**********************************************************************/
1412 #define DECODE_BUFFER_SIZE 2048
1413
1414 /** Decoder state 
1415  @brief Decoder state
1416  */
1417 struct CELTDecoder {
1418    const CELTMode *mode;
1419    int overlap;
1420    int channels;
1421
1422    int start, end;
1423
1424    /* Everything beyond this point gets cleared on a reset */
1425 #define DECODER_RESET_START last_pitch_index
1426
1427    int last_pitch_index;
1428    int loss_count;
1429    int postfilter_period;
1430    int postfilter_period_old;
1431    celt_word16 postfilter_gain;
1432    celt_word16 postfilter_gain_old;
1433
1434    celt_sig preemph_memD[2];
1435    
1436    celt_sig _decode_mem[1]; /* Size = channels*(DECODE_BUFFER_SIZE+mode->overlap) */
1437    /* celt_word16 lpc[],  Size = channels*LPC_ORDER */
1438    /* celt_word16 oldEBands[], Size = channels*mode->nbEBands */
1439 };
1440
1441 int celt_decoder_get_size(const CELTMode *mode, int channels)
1442 {
1443    int size = sizeof(struct CELTDecoder)
1444             + (channels*(DECODE_BUFFER_SIZE+mode->overlap)-1)*sizeof(celt_sig)
1445             + channels*LPC_ORDER*sizeof(celt_word16)
1446             + channels*mode->nbEBands*sizeof(celt_word16);
1447    return size;
1448 }
1449
1450 CELTDecoder *celt_decoder_create(const CELTMode *mode, int channels, int *error)
1451 {
1452    return celt_decoder_init(
1453          (CELTDecoder *)celt_alloc(celt_decoder_get_size(mode, channels)),
1454          mode, channels, error);
1455 }
1456
1457 CELTDecoder *celt_decoder_init(CELTDecoder *st, const CELTMode *mode, int channels, int *error)
1458 {
1459    if (channels < 0 || channels > 2)
1460    {
1461       if (error)
1462          *error = CELT_BAD_ARG;
1463       return NULL;
1464    }
1465
1466    if (st==NULL)
1467    {
1468       if (error)
1469          *error = CELT_ALLOC_FAIL;
1470       return NULL;
1471    }
1472
1473    CELT_MEMSET((char*)st, 0, celt_decoder_get_size(mode, channels));
1474
1475    st->mode = mode;
1476    st->overlap = mode->overlap;
1477    st->channels = channels;
1478
1479    st->start = 0;
1480    st->end = st->mode->effEBands;
1481
1482    st->loss_count = 0;
1483
1484    if (error)
1485       *error = CELT_OK;
1486    return st;
1487 }
1488
1489 void celt_decoder_destroy(CELTDecoder *st)
1490 {
1491    celt_free(st);
1492 }
1493
1494 static void celt_decode_lost(CELTDecoder * restrict st, celt_word16 * restrict pcm, int N, int LM)
1495 {
1496    int c;
1497    int pitch_index;
1498    int overlap = st->mode->overlap;
1499    celt_word16 fade = Q15ONE;
1500    int i, len;
1501    const int C = CHANNELS(st->channels);
1502    int offset;
1503    celt_sig *out_mem[2];
1504    celt_sig *decode_mem[2];
1505    celt_sig *overlap_mem[2];
1506    celt_word16 *lpc;
1507    SAVE_STACK;
1508    
1509    c=0; do {
1510       decode_mem[c] = st->_decode_mem + c*(DECODE_BUFFER_SIZE+st->overlap);
1511       out_mem[c] = decode_mem[c]+DECODE_BUFFER_SIZE-MAX_PERIOD;
1512       overlap_mem[c] = decode_mem[c]+DECODE_BUFFER_SIZE;
1513    } while (++c<C);
1514    lpc = (celt_word16*)(st->_decode_mem+(DECODE_BUFFER_SIZE+st->overlap)*C);
1515
1516    len = N+st->mode->overlap;
1517    
1518    if (st->loss_count == 0)
1519    {
1520       celt_word16 pitch_buf[MAX_PERIOD>>1];
1521       celt_word32 tmp=0;
1522       celt_word32 mem0[2]={0,0};
1523       celt_word16 mem1[2]={0,0};
1524       int len2 = len;
1525       /* FIXME: This is a kludge */
1526       if (len2>MAX_PERIOD>>1)
1527          len2 = MAX_PERIOD>>1;
1528       pitch_downsample(out_mem, pitch_buf, MAX_PERIOD, MAX_PERIOD,
1529                        C, mem0, mem1);
1530       pitch_search(st->mode, pitch_buf+((MAX_PERIOD-len2)>>1), pitch_buf, len2,
1531                    MAX_PERIOD-len2-100, &pitch_index, &tmp, 1<<LM);
1532       pitch_index = MAX_PERIOD-len2-pitch_index;
1533       st->last_pitch_index = pitch_index;
1534    } else {
1535       pitch_index = st->last_pitch_index;
1536       if (st->loss_count < 5)
1537          fade = QCONST16(.8f,15);
1538       else
1539          fade = 0;
1540    }
1541
1542    c=0; do {
1543       /* FIXME: This is more memory than necessary */
1544       celt_word32 e[2*MAX_PERIOD];
1545       celt_word16 exc[2*MAX_PERIOD];
1546       celt_word32 ac[LPC_ORDER+1];
1547       celt_word16 decay = 1;
1548       celt_word32 S1=0;
1549       celt_word16 mem[LPC_ORDER]={0};
1550
1551       offset = MAX_PERIOD-pitch_index;
1552       for (i=0;i<MAX_PERIOD;i++)
1553          exc[i] = ROUND16(out_mem[c][i], SIG_SHIFT);
1554
1555       if (st->loss_count == 0)
1556       {
1557          _celt_autocorr(exc, ac, st->mode->window, st->mode->overlap,
1558                         LPC_ORDER, MAX_PERIOD);
1559
1560          /* Noise floor -40 dB */
1561 #ifdef FIXED_POINT
1562          ac[0] += SHR32(ac[0],13);
1563 #else
1564          ac[0] *= 1.0001f;
1565 #endif
1566          /* Lag windowing */
1567          for (i=1;i<=LPC_ORDER;i++)
1568          {
1569             /*ac[i] *= exp(-.5*(2*M_PI*.002*i)*(2*M_PI*.002*i));*/
1570 #ifdef FIXED_POINT
1571             ac[i] -= MULT16_32_Q15(2*i*i, ac[i]);
1572 #else
1573             ac[i] -= ac[i]*(.008f*i)*(.008f*i);
1574 #endif
1575          }
1576
1577          _celt_lpc(lpc+c*LPC_ORDER, ac, LPC_ORDER);
1578       }
1579       for (i=0;i<LPC_ORDER;i++)
1580          mem[i] = ROUND16(out_mem[c][MAX_PERIOD-1-i], SIG_SHIFT);
1581       fir(exc, lpc+c*LPC_ORDER, exc, MAX_PERIOD, LPC_ORDER, mem);
1582       /*for (i=0;i<MAX_PERIOD;i++)printf("%d ", exc[i]); printf("\n");*/
1583       /* Check if the waveform is decaying (and if so how fast) */
1584       {
1585          celt_word32 E1=1, E2=1;
1586          int period;
1587          if (pitch_index <= MAX_PERIOD/2)
1588             period = pitch_index;
1589          else
1590             period = MAX_PERIOD/2;
1591          for (i=0;i<period;i++)
1592          {
1593             E1 += SHR32(MULT16_16(exc[MAX_PERIOD-period+i],exc[MAX_PERIOD-period+i]),8);
1594             E2 += SHR32(MULT16_16(exc[MAX_PERIOD-2*period+i],exc[MAX_PERIOD-2*period+i]),8);
1595          }
1596          if (E1 > E2)
1597             E1 = E2;
1598          decay = celt_sqrt(frac_div32(SHR(E1,1),E2));
1599       }
1600
1601       /* Copy excitation, taking decay into account */
1602       for (i=0;i<len+st->mode->overlap;i++)
1603       {
1604          celt_word16 tmp;
1605          if (offset+i >= MAX_PERIOD)
1606          {
1607             offset -= pitch_index;
1608             decay = MULT16_16_Q15(decay, decay);
1609          }
1610          e[i] = SHL32(EXTEND32(MULT16_16_Q15(decay, exc[offset+i])), SIG_SHIFT);
1611          tmp = ROUND16(out_mem[c][offset+i],SIG_SHIFT);
1612          S1 += SHR32(MULT16_16(tmp,tmp),8);
1613       }
1614       for (i=0;i<LPC_ORDER;i++)
1615          mem[i] = ROUND16(out_mem[c][MAX_PERIOD-1-i], SIG_SHIFT);
1616       for (i=0;i<len+st->mode->overlap;i++)
1617          e[i] = MULT16_32_Q15(fade, e[i]);
1618       iir(e, lpc+c*LPC_ORDER, e, len+st->mode->overlap, LPC_ORDER, mem);
1619
1620       {
1621          celt_word32 S2=0;
1622          for (i=0;i<len+overlap;i++)
1623          {
1624             celt_word16 tmp = ROUND16(e[i],SIG_SHIFT);
1625             S2 += SHR32(MULT16_16(tmp,tmp),8);
1626          }
1627          /* This checks for an "explosion" in the synthesis */
1628 #ifdef FIXED_POINT
1629          if (!(S1 > SHR32(S2,2)))
1630 #else
1631          /* Float test is written this way to catch NaNs at the same time */
1632          if (!(S1 > 0.2f*S2))
1633 #endif
1634          {
1635             for (i=0;i<len+overlap;i++)
1636                e[i] = 0;
1637          } else if (S1 < S2)
1638          {
1639             celt_word16 ratio = celt_sqrt(frac_div32(SHR32(S1,1)+1,S2+1));
1640             for (i=0;i<len+overlap;i++)
1641                e[i] = MULT16_32_Q15(ratio, e[i]);
1642          }
1643       }
1644
1645 #ifdef ENABLE_POSTFILTER
1646       /* Apply post-filter to the MDCT overlap of the previous frame */
1647       comb_filter(out_mem[c]+MAX_PERIOD, out_mem[c]+MAX_PERIOD, st->postfilter_period, st->postfilter_period, st->overlap, C,
1648                   st->postfilter_gain, st->postfilter_gain, NULL, 0);
1649 #endif /* ENABLE_POSTFILTER */
1650
1651       for (i=0;i<MAX_PERIOD+st->mode->overlap-N;i++)
1652          out_mem[c][i] = out_mem[c][N+i];
1653
1654       /* Apply TDAC to the concealed audio so that it blends with the
1655          previous and next frames */
1656       for (i=0;i<overlap/2;i++)
1657       {
1658          celt_word32 tmp;
1659          tmp = MULT16_32_Q15(st->mode->window[i],           e[N+overlap-1-i]) +
1660                MULT16_32_Q15(st->mode->window[overlap-i-1], e[N+i          ]);
1661          out_mem[c][MAX_PERIOD+i] = MULT16_32_Q15(st->mode->window[overlap-i-1], tmp);
1662          out_mem[c][MAX_PERIOD+overlap-i-1] = MULT16_32_Q15(st->mode->window[i], tmp);
1663       }
1664       for (i=0;i<N;i++)
1665          out_mem[c][MAX_PERIOD-N+i] = e[i];
1666
1667 #ifdef ENABLE_POSTFILTER
1668       /* Apply pre-filter to the MDCT overlap for the next frame (post-filter will be applied then) */
1669       comb_filter(e, out_mem[c]+MAX_PERIOD, st->postfilter_period, st->postfilter_period, st->overlap, C,
1670                   -st->postfilter_gain, -st->postfilter_gain, NULL, 0);
1671 #endif /* ENABLE_POSTFILTER */
1672       for (i=0;i<overlap;i++)
1673          out_mem[c][MAX_PERIOD+i] = e[i];
1674    } while (++c<C);
1675
1676    {
1677       celt_word32 *out_syn[2];
1678       out_syn[0] = out_mem[0]+MAX_PERIOD-N;
1679       if (C==2)
1680          out_syn[1] = out_mem[1]+MAX_PERIOD-N;
1681       deemphasis(out_syn, pcm, N, C, st->mode->preemph, st->preemph_memD);
1682    }
1683    
1684    st->loss_count++;
1685
1686    RESTORE_STACK;
1687 }
1688
1689 #ifdef FIXED_POINT
1690 int celt_decode_with_ec(CELTDecoder * restrict st, const unsigned char *data, int len, celt_int16 * restrict pcm, int frame_size, ec_dec *dec)
1691 {
1692 #else
1693 int celt_decode_with_ec_float(CELTDecoder * restrict st, const unsigned char *data, int len, celt_sig * restrict pcm, int frame_size, ec_dec *dec)
1694 {
1695 #endif
1696    int c, i, N;
1697    int spread_decision;
1698    int bits;
1699    ec_dec _dec;
1700    ec_byte_buffer buf;
1701    VARDECL(celt_sig, freq);
1702    VARDECL(celt_norm, X);
1703    VARDECL(celt_ener, bandE);
1704    VARDECL(int, fine_quant);
1705    VARDECL(int, pulses);
1706    VARDECL(int, offsets);
1707    VARDECL(int, fine_priority);
1708    VARDECL(int, tf_res);
1709    celt_sig *out_mem[2];
1710    celt_sig *decode_mem[2];
1711    celt_sig *overlap_mem[2];
1712    celt_sig *out_syn[2];
1713    celt_word16 *lpc;
1714    celt_word16 *oldBandE;
1715
1716    int shortBlocks;
1717    int isTransient;
1718    int intra_ener;
1719    const int C = CHANNELS(st->channels);
1720    int LM, M;
1721    int nbFilledBytes, nbAvailableBytes;
1722    int effEnd;
1723    int codedBands;
1724    int alloc_trim;
1725    int postfilter_pitch;
1726    celt_word16 postfilter_gain;
1727    int intensity=0;
1728    int dual_stereo=0;
1729    SAVE_STACK;
1730
1731    if (pcm==NULL)
1732       return CELT_BAD_ARG;
1733
1734    for (LM=0;LM<4;LM++)
1735       if (st->mode->shortMdctSize<<LM==frame_size)
1736          break;
1737    if (LM>=MAX_CONFIG_SIZES)
1738       return CELT_BAD_ARG;
1739    M=1<<LM;
1740
1741    c=0; do {
1742       decode_mem[c] = st->_decode_mem + c*(DECODE_BUFFER_SIZE+st->overlap);
1743       out_mem[c] = decode_mem[c]+DECODE_BUFFER_SIZE-MAX_PERIOD;
1744       overlap_mem[c] = decode_mem[c]+DECODE_BUFFER_SIZE;
1745    } while (++c<C);
1746    lpc = (celt_word16*)(st->_decode_mem+(DECODE_BUFFER_SIZE+st->overlap)*C);
1747    oldBandE = lpc+C*LPC_ORDER;
1748
1749    N = M*st->mode->shortMdctSize;
1750
1751    effEnd = st->end;
1752    if (effEnd > st->mode->effEBands)
1753       effEnd = st->mode->effEBands;
1754
1755    ALLOC(freq, C*N, celt_sig); /**< Interleaved signal MDCTs */
1756    ALLOC(X, C*N, celt_norm);   /**< Interleaved normalised MDCTs */
1757    ALLOC(bandE, st->mode->nbEBands*C, celt_ener);
1758    c=0; do
1759       for (i=0;i<M*st->mode->eBands[st->start];i++)
1760          X[c*N+i] = 0;
1761    while (++c<C);
1762    c=0; do   
1763       for (i=M*st->mode->eBands[effEnd];i<N;i++)
1764          X[c*N+i] = 0;
1765    while (++c<C);
1766
1767    if (data == NULL)
1768    {
1769       celt_decode_lost(st, pcm, N, LM);
1770       RESTORE_STACK;
1771       return CELT_OK;
1772    }
1773    if (len<0) {
1774      RESTORE_STACK;
1775      return CELT_BAD_ARG;
1776    }
1777    
1778    if (dec == NULL)
1779    {
1780       ec_byte_readinit(&buf,(unsigned char*)data,len);
1781       ec_dec_init(&_dec,&buf);
1782       dec = &_dec;
1783       nbFilledBytes = 0;
1784    } else {
1785       nbFilledBytes = (ec_dec_tell(dec, 0)+4)>>3;
1786    }
1787    nbAvailableBytes = len-nbFilledBytes;
1788
1789    if (ec_dec_bit_prob(dec, 32768))
1790    {
1791 #ifdef ENABLE_POSTFILTER
1792       int qg, octave;
1793       octave = ec_dec_uint(dec, 6);
1794       postfilter_pitch = (16<<octave)+ec_dec_bits(dec, 4+octave);
1795       qg = ec_dec_bits(dec, 2);
1796       postfilter_gain = QCONST16(.125f,15)*(qg+2);
1797 #else /* ENABLE_POSTFILTER */
1798       RESTORE_STACK;
1799       return CELT_CORRUPTED_DATA;
1800 #endif /* ENABLE_POSTFILTER */
1801
1802    } else {
1803       postfilter_gain = 0;
1804       postfilter_pitch = 0;
1805    }
1806
1807    /* Decode the global flags (first symbols in the stream) */
1808    intra_ener = ec_dec_bit_prob(dec, 8192);
1809    /* Get band energies */
1810    unquant_coarse_energy(st->mode, st->start, st->end, bandE, oldBandE,
1811          intra_ener, dec, C, LM);
1812
1813    if (LM > 0)
1814       isTransient = ec_dec_bit_prob(dec, 8192);
1815    else
1816       isTransient = 0;
1817
1818    if (isTransient)
1819       shortBlocks = M;
1820    else
1821       shortBlocks = 0;
1822
1823    ALLOC(tf_res, st->mode->nbEBands, int);
1824    tf_decode(st->start, st->end, C, isTransient, tf_res, LM, dec);
1825
1826    spread_decision = ec_dec_cdf(dec, spread_cdf, 5);
1827
1828    ALLOC(pulses, st->mode->nbEBands, int);
1829    ALLOC(offsets, st->mode->nbEBands, int);
1830    ALLOC(fine_priority, st->mode->nbEBands, int);
1831
1832    for (i=0;i<st->mode->nbEBands;i++)
1833       offsets[i] = 0;
1834    for (i=0;i<st->mode->nbEBands;i++)
1835    {
1836       if (ec_dec_bit_prob(dec, 1024))
1837       {
1838          while (ec_dec_bit_prob(dec, 32768))
1839             offsets[i]++;
1840          offsets[i]++;
1841          offsets[i] *= (6<<BITRES);
1842       }
1843    }
1844
1845    ALLOC(fine_quant, st->mode->nbEBands, int);
1846    alloc_trim = ec_dec_cdf(dec, trim_cdf, 7);
1847
1848    if (C==2)
1849    {
1850       dual_stereo = ec_dec_bit_prob(dec, 32768);
1851       intensity = ec_dec_uint(dec, 1+st->end-st->start);
1852    }
1853
1854    bits = (len*8<<BITRES) - ec_dec_tell(dec, BITRES) - 1;
1855    codedBands = compute_allocation(st->mode, st->start, st->end, offsets,
1856          alloc_trim, bits, pulses, fine_quant, fine_priority, C, LM, dec, 0, 0);
1857    
1858    unquant_fine_energy(st->mode, st->start, st->end, bandE, oldBandE, fine_quant, dec, C);
1859
1860    /* Decode fixed codebook */
1861    quant_all_bands(0, st->mode, st->start, st->end, X, C==2 ? X+N : NULL,
1862          NULL, pulses, shortBlocks, spread_decision, dual_stereo, intensity, tf_res, 1,
1863          len*8, dec, LM, codedBands);
1864
1865    unquant_energy_finalise(st->mode, st->start, st->end, bandE, oldBandE,
1866          fine_quant, fine_priority, len*8-ec_dec_tell(dec, 0), dec, C);
1867
1868    log2Amp(st->mode, st->start, st->end, bandE, oldBandE, C);
1869
1870    /* Synthesis */
1871    denormalise_bands(st->mode, X, freq, bandE, effEnd, C, M);
1872
1873    CELT_MOVE(decode_mem[0], decode_mem[0]+N, DECODE_BUFFER_SIZE-N);
1874    if (C==2)
1875       CELT_MOVE(decode_mem[1], decode_mem[1]+N, DECODE_BUFFER_SIZE-N);
1876
1877    c=0; do
1878       for (i=0;i<M*st->mode->eBands[st->start];i++)
1879          freq[c*N+i] = 0;
1880    while (++c<C);
1881    c=0; do
1882       for (i=M*st->mode->eBands[effEnd];i<N;i++)
1883          freq[c*N+i] = 0;
1884    while (++c<C);
1885
1886    out_syn[0] = out_mem[0]+MAX_PERIOD-N;
1887    if (C==2)
1888       out_syn[1] = out_mem[1]+MAX_PERIOD-N;
1889
1890    /* Compute inverse MDCTs */
1891    compute_inv_mdcts(st->mode, shortBlocks, freq, out_syn, overlap_mem, C, LM);
1892
1893 #ifdef ENABLE_POSTFILTER
1894    c=0; do {
1895       st->postfilter_period=IMAX(st->postfilter_period, COMBFILTER_MINPERIOD);
1896       st->postfilter_period_old=IMAX(st->postfilter_period_old, COMBFILTER_MINPERIOD);
1897       if (LM!=0)
1898       {
1899          comb_filter(out_syn[c], out_syn[c], st->postfilter_period, st->postfilter_period, st->overlap, C,
1900                st->postfilter_gain, st->postfilter_gain, NULL, 0);
1901          comb_filter(out_syn[c]+st->overlap, out_syn[c]+st->overlap, st->postfilter_period, postfilter_pitch, N-st->overlap, C,
1902                st->postfilter_gain, postfilter_gain, st->mode->window, st->mode->overlap);
1903       } else {
1904          comb_filter(out_syn[c], out_syn[c], st->postfilter_period_old, st->postfilter_period, N-st->overlap, C,
1905                st->postfilter_gain_old, st->postfilter_gain, st->mode->window, st->mode->overlap);
1906       }
1907    } while (++c<C);
1908    st->postfilter_period_old = st->postfilter_period;
1909    st->postfilter_gain_old = st->postfilter_gain;
1910    st->postfilter_period = postfilter_pitch;
1911    st->postfilter_gain = postfilter_gain;
1912 #endif /* ENABLE_POSTFILTER */
1913
1914    deemphasis(out_syn, pcm, N, C, st->mode->preemph, st->preemph_memD);
1915    st->loss_count = 0;
1916    RESTORE_STACK;
1917    if (ec_dec_get_error(dec))
1918       return CELT_CORRUPTED_DATA;
1919    else
1920       return CELT_OK;
1921 }
1922
1923 #ifdef FIXED_POINT
1924 #ifndef DISABLE_FLOAT_API
1925 int celt_decode_with_ec_float(CELTDecoder * restrict st, const unsigned char *data, int len, float * restrict pcm, int frame_size, ec_dec *dec)
1926 {
1927    int j, ret, C, N, LM, M;
1928    VARDECL(celt_int16, out);
1929    SAVE_STACK;
1930
1931    if (pcm==NULL)
1932       return CELT_BAD_ARG;
1933
1934    for (LM=0;LM<4;LM++)
1935       if (st->mode->shortMdctSize<<LM==frame_size)
1936          break;
1937    if (LM>=MAX_CONFIG_SIZES)
1938       return CELT_BAD_ARG;
1939    M=1<<LM;
1940
1941    C = CHANNELS(st->channels);
1942    N = M*st->mode->shortMdctSize;
1943    
1944    ALLOC(out, C*N, celt_int16);
1945    ret=celt_decode_with_ec(st, data, len, out, frame_size, dec);
1946    if (ret==0)
1947       for (j=0;j<C*N;j++)
1948          pcm[j]=out[j]*(1.f/32768.f);
1949      
1950    RESTORE_STACK;
1951    return ret;
1952 }
1953 #endif /*DISABLE_FLOAT_API*/
1954 #else
1955 int celt_decode_with_ec(CELTDecoder * restrict st, const unsigned char *data, int len, celt_int16 * restrict pcm, int frame_size, ec_dec *dec)
1956 {
1957    int j, ret, C, N, LM, M;
1958    VARDECL(celt_sig, out);
1959    SAVE_STACK;
1960
1961    if (pcm==NULL)
1962       return CELT_BAD_ARG;
1963
1964    for (LM=0;LM<4;LM++)
1965       if (st->mode->shortMdctSize<<LM==frame_size)
1966          break;
1967    if (LM>=MAX_CONFIG_SIZES)
1968       return CELT_BAD_ARG;
1969    M=1<<LM;
1970
1971    C = CHANNELS(st->channels);
1972    N = M*st->mode->shortMdctSize;
1973    ALLOC(out, C*N, celt_sig);
1974
1975    ret=celt_decode_with_ec_float(st, data, len, out, frame_size, dec);
1976
1977    if (ret==0)
1978       for (j=0;j<C*N;j++)
1979          pcm[j] = FLOAT2INT16 (out[j]);
1980    
1981    RESTORE_STACK;
1982    return ret;
1983 }
1984 #endif
1985
1986 int celt_decode(CELTDecoder * restrict st, const unsigned char *data, int len, celt_int16 * restrict pcm, int frame_size)
1987 {
1988    return celt_decode_with_ec(st, data, len, pcm, frame_size, NULL);
1989 }
1990
1991 #ifndef DISABLE_FLOAT_API
1992 int celt_decode_float(CELTDecoder * restrict st, const unsigned char *data, int len, float * restrict pcm, int frame_size)
1993 {
1994    return celt_decode_with_ec_float(st, data, len, pcm, frame_size, NULL);
1995 }
1996 #endif /* DISABLE_FLOAT_API */
1997
1998 int celt_decoder_ctl(CELTDecoder * restrict st, int request, ...)
1999 {
2000    va_list ap;
2001
2002    va_start(ap, request);
2003    switch (request)
2004    {
2005       case CELT_GET_MODE_REQUEST:
2006       {
2007          const CELTMode ** value = va_arg(ap, const CELTMode**);
2008          if (value==0)
2009             goto bad_arg;
2010          *value=st->mode;
2011       }
2012       break;
2013       case CELT_SET_START_BAND_REQUEST:
2014       {
2015          celt_int32 value = va_arg(ap, celt_int32);
2016          if (value<0 || value>=st->mode->nbEBands)
2017             goto bad_arg;
2018          st->start = value;
2019       }
2020       break;
2021       case CELT_SET_END_BAND_REQUEST:
2022       {
2023          celt_int32 value = va_arg(ap, celt_int32);
2024          if (value<0 || value>=st->mode->nbEBands)
2025             goto bad_arg;
2026          st->end = value;
2027       }
2028       break;
2029       case CELT_RESET_STATE:
2030       {
2031          CELT_MEMSET((char*)&st->DECODER_RESET_START, 0,
2032                celt_decoder_get_size(st->mode, st->channels)-
2033                ((char*)&st->DECODER_RESET_START - (char*)st));
2034       }
2035       break;
2036       default:
2037          goto bad_request;
2038    }
2039    va_end(ap);
2040    return CELT_OK;
2041 bad_arg:
2042    va_end(ap);
2043    return CELT_BAD_ARG;
2044 bad_request:
2045       va_end(ap);
2046   return CELT_UNIMPLEMENTED;
2047 }
2048
2049 const char *celt_strerror(int error)
2050 {
2051    static const char *error_strings[8] = {
2052       "success",
2053       "invalid argument",
2054       "invalid mode",
2055       "internal error",
2056       "corrupted stream",
2057       "request not implemented",
2058       "invalid state",
2059       "memory allocation failed"
2060    };
2061    if (error > 0 || error < -7)
2062       return "unknown error";
2063    else 
2064       return error_strings[-error];
2065 }
2066