Re-organize spreading/folding constants.
[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
58 #define COMBFILTER_MAXPERIOD 1024
59 #define COMBFILTER_MINPERIOD 16
60
61 /** Encoder state 
62  @brief Encoder state
63  */
64 struct CELTEncoder {
65    const CELTMode *mode;     /**< Mode used by the encoder */
66    int overlap;
67    int channels;
68    
69    int force_intra;
70    int complexity;
71    int start, end;
72
73    celt_int32 vbr_rate_norm; /* Target number of 8th bits per frame */
74    int constrained_vbr;      /* If zero, VBR can do whatever it likes with the rate */
75
76    /* Everything beyond this point gets cleared on a reset */
77 #define ENCODER_RESET_START frame_max
78
79    celt_word32 frame_max;
80    int spread_decision;
81    int delayedIntra;
82    int tonal_average;
83    int lastCodedBands;
84
85    int prefilter_period;
86    celt_word16 prefilter_gain;
87 #ifdef RESYNTH
88    int prefilter_period_old;
89    celt_word16 prefilter_gain_old;
90 #endif
91
92    /* VBR-related parameters */
93    celt_int32 vbr_reservoir;
94    celt_int32 vbr_drift;
95    celt_int32 vbr_offset;
96    celt_int32 vbr_count;
97
98    celt_word32 preemph_memE[2];
99    celt_word32 preemph_memD[2];
100
101 #ifdef RESYNTH
102    celt_sig syn_mem[2][2*MAX_PERIOD];
103 #endif
104
105    celt_sig in_mem[1]; /* Size = channels*mode->overlap */
106    /* celt_sig prefilter_mem[],  Size = channels*COMBFILTER_PERIOD */
107    /* celt_sig overlap_mem[],  Size = channels*mode->overlap */
108    /* celt_word16 oldEBands[], Size = channels*mode->nbEBands */
109 };
110
111 int celt_encoder_get_size(const CELTMode *mode, int channels)
112 {
113    int size = sizeof(struct CELTEncoder)
114          + (2*channels*mode->overlap-1)*sizeof(celt_sig)
115          + channels*COMBFILTER_MAXPERIOD*sizeof(celt_sig)
116          + channels*mode->nbEBands*sizeof(celt_word16);
117    return size;
118 }
119
120 CELTEncoder *celt_encoder_create(const CELTMode *mode, int channels, int *error)
121 {
122    return celt_encoder_init(
123          (CELTEncoder *)celt_alloc(celt_encoder_get_size(mode, channels)),
124          mode, channels, error);
125 }
126
127 CELTEncoder *celt_encoder_init(CELTEncoder *st, const CELTMode *mode, int channels, int *error)
128 {
129    if (channels < 0 || channels > 2)
130    {
131       if (error)
132          *error = CELT_BAD_ARG;
133       return NULL;
134    }
135
136    if (st==NULL)
137    {
138       if (error)
139          *error = CELT_ALLOC_FAIL;
140       return NULL;
141    }
142
143    CELT_MEMSET((char*)st, 0, celt_encoder_get_size(mode, channels));
144    
145    st->mode = mode;
146    st->overlap = mode->overlap;
147    st->channels = channels;
148
149    st->start = 0;
150    st->end = st->mode->effEBands;
151    st->constrained_vbr = 1;
152
153    st->vbr_rate_norm = 0;
154    st->vbr_offset = 0;
155    st->force_intra  = 0;
156    st->delayedIntra = 1;
157    st->tonal_average = 256;
158    st->spread_decision = SPREAD_NORMAL;
159    st->complexity = 5;
160
161    if (error)
162       *error = CELT_OK;
163    return st;
164 }
165
166 void celt_encoder_destroy(CELTEncoder *st)
167 {
168    celt_free(st);
169 }
170
171 static inline celt_int16 FLOAT2INT16(float x)
172 {
173    x = x*CELT_SIG_SCALE;
174    x = MAX32(x, -32768);
175    x = MIN32(x, 32767);
176    return (celt_int16)float2int(x);
177 }
178
179 static inline celt_word16 SIG2WORD16(celt_sig x)
180 {
181 #ifdef FIXED_POINT
182    x = PSHR32(x, SIG_SHIFT);
183    x = MAX32(x, -32768);
184    x = MIN32(x, 32767);
185    return EXTRACT16(x);
186 #else
187    return (celt_word16)x;
188 #endif
189 }
190
191 static int transient_analysis(const celt_word32 * restrict in, int len, int C,
192                               celt_word32 *frame_max, int overlap)
193 {
194    int i;
195    VARDECL(celt_word16, tmp);
196    celt_word32 mem0=0,mem1=0;
197    int is_transient = 0;
198    int block;
199    int N;
200    /* FIXME: Make that smaller */
201    celt_word16 bins[50];
202    SAVE_STACK;
203    ALLOC(tmp, len, celt_word16);
204
205    block = overlap/2;
206    N=len/block;
207    if (C==1)
208    {
209       for (i=0;i<len;i++)
210          tmp[i] = SHR32(in[i],SIG_SHIFT);
211    } else {
212       for (i=0;i<len;i++)
213          tmp[i] = SHR32(ADD32(in[i],in[i+len]), SIG_SHIFT+1);
214    }
215
216    /* High-pass filter: (1 - 2*z^-1 + z^-2) / (1 - z^-1 + .5*z^-2) */
217    for (i=0;i<len;i++)
218    {
219       celt_word32 x,y;
220       x = tmp[i];
221       y = ADD32(mem0, x);
222 #ifdef FIXED_POINT
223       mem0 = mem1 + y - SHL32(x,1);
224       mem1 = x - SHR32(y,1);
225 #else
226       mem0 = mem1 + y - 2*x;
227       mem1 = x - .5f*y;
228 #endif
229       tmp[i] = EXTRACT16(SHR(y,2));
230    }
231    /* First few samples are bad because we don't propagate the memory */
232    for (i=0;i<12;i++)
233       tmp[i] = 0;
234
235    for (i=0;i<N;i++)
236    {
237       int j;
238       float max_abs=0;
239       for (j=0;j<block;j++)
240          max_abs = MAX32(max_abs, tmp[i*block+j]);
241       bins[i] = max_abs;
242    }
243    for (i=0;i<N;i++)
244    {
245       int j;
246       int conseq=0;
247       celt_word16 t1, t2, t3;
248
249       t1 = MULT16_16_Q15(QCONST16(.15f, 15), bins[i]);
250       t2 = MULT16_16_Q15(QCONST16(.4f, 15), bins[i]);
251       t3 = MULT16_16_Q15(QCONST16(.15f, 15), bins[i]);
252       for (j=0;j<i;j++)
253       {
254          if (bins[j] < t1)
255             conseq++;
256          if (bins[j] < t2)
257             conseq++;
258          else
259             conseq = 0;
260       }
261       if (conseq>=3)
262          is_transient=1;
263       conseq = 0;
264       for (j=i+1;j<N;j++)
265       {
266          if (bins[j] < t3)
267             conseq++;
268          else
269             conseq = 0;
270       }
271       if (conseq>=7)
272          is_transient=1;
273    }
274    RESTORE_STACK;
275    return is_transient;
276 }
277
278 /** Apply window and compute the MDCT for all sub-frames and 
279     all channels in a frame */
280 static void compute_mdcts(const CELTMode *mode, int shortBlocks, celt_sig * restrict in, celt_sig * restrict out, int _C, int LM)
281 {
282    const int C = CHANNELS(_C);
283    if (C==1 && !shortBlocks)
284    {
285       const int overlap = OVERLAP(mode);
286       clt_mdct_forward(&mode->mdct, in, out, mode->window, overlap, mode->maxLM-LM);
287    } else {
288       const int overlap = OVERLAP(mode);
289       int N = mode->shortMdctSize<<LM;
290       int B = 1;
291       int b, c;
292       VARDECL(celt_word32, tmp);
293       SAVE_STACK;
294       if (shortBlocks)
295       {
296          /*lookup = &mode->mdct[0];*/
297          N = mode->shortMdctSize;
298          B = shortBlocks;
299       }
300       ALLOC(tmp, N, celt_word32);
301       c=0; do {
302          for (b=0;b<B;b++)
303          {
304             int j;
305             clt_mdct_forward(&mode->mdct, in+c*(B*N+overlap)+b*N, tmp, mode->window, overlap, shortBlocks ? mode->maxLM : mode->maxLM-LM);
306             /* Interleaving the sub-frames */
307             for (j=0;j<N;j++)
308                out[(j*B+b)+c*N*B] = tmp[j];
309          }
310       } while (++c<C);
311       RESTORE_STACK;
312    }
313 }
314
315 /** Compute the IMDCT and apply window for all sub-frames and 
316     all channels in a frame */
317 static void compute_inv_mdcts(const CELTMode *mode, int shortBlocks, celt_sig *X,
318       celt_sig * restrict out_mem[],
319       celt_sig * restrict overlap_mem[], int _C, int LM)
320 {
321    int c;
322    const int C = CHANNELS(_C);
323    const int N = mode->shortMdctSize<<LM;
324    const int overlap = OVERLAP(mode);
325    c=0; do {
326       int j;
327          VARDECL(celt_word32, x);
328          VARDECL(celt_word32, tmp);
329          int b;
330          int N2 = N;
331          int B = 1;
332          SAVE_STACK;
333          
334          ALLOC(x, N+overlap, celt_word32);
335          ALLOC(tmp, N, celt_word32);
336
337          if (shortBlocks)
338          {
339             N2 = mode->shortMdctSize;
340             B = shortBlocks;
341          }
342          /* Prevents problems from the imdct doing the overlap-add */
343          CELT_MEMSET(x, 0, overlap);
344
345          for (b=0;b<B;b++)
346          {
347             /* De-interleaving the sub-frames */
348             for (j=0;j<N2;j++)
349                tmp[j] = X[(j*B+b)+c*N2*B];
350             clt_mdct_backward(&mode->mdct, tmp, x+N2*b, mode->window, overlap, shortBlocks ? mode->maxLM : mode->maxLM-LM);
351          }
352
353          for (j=0;j<overlap;j++)
354             out_mem[c][j] = x[j] + overlap_mem[c][j];
355          for (;j<N;j++)
356             out_mem[c][j] = x[j];
357          for (j=0;j<overlap;j++)
358             overlap_mem[c][j] = x[N+j];
359          RESTORE_STACK;
360    } while (++c<C);
361 }
362
363 static void deemphasis(celt_sig *in[], celt_word16 *pcm, int N, int _C, const celt_word16 *coef, celt_sig *mem)
364 {
365    const int C = CHANNELS(_C);
366    int c;
367    c=0; do {
368       int j;
369       celt_sig * restrict x;
370       celt_word16  * restrict y;
371       celt_sig m = mem[c];
372       x =in[c];
373       y = pcm+c;
374       for (j=0;j<N;j++)
375       {
376          celt_sig tmp = *x + m;
377          m = MULT16_32_Q15(coef[0], tmp)
378            - MULT16_32_Q15(coef[1], *x);
379          tmp = SHL32(MULT16_32_Q15(coef[3], tmp), 2);
380          *y = SCALEOUT(SIG2WORD16(tmp));
381          x++;
382          y+=C;
383       }
384       mem[c] = m;
385    } while (++c<C);
386 }
387
388 #ifdef ENABLE_POSTFILTER
389 /* FIXME: Handle the case where T = maxperiod */
390 static void comb_filter(celt_word32 *y, celt_word32 *x, int T0, int T1, int N,
391       int C, celt_word16 g0, celt_word16 g1, const celt_word16 *window, int overlap)
392 {
393    int i;
394    /* printf ("%d %d %f %f\n", T0, T1, g0, g1); */
395    celt_word16 g00, g01, g02, g10, g11, g12;
396    celt_word16 t0, t1, t2;
397    /* zeros at theta = +/- 5*pi/6 */
398    t0 = QCONST16(.26795f, 15);
399    t1 = QCONST16(.46410f, 15);
400    t2 = QCONST16(.26795f, 15);
401    g00 = MULT16_16_Q15(g0, t0);
402    g01 = MULT16_16_Q15(g0, t1);
403    g02 = MULT16_16_Q15(g0, t2);
404    g10 = MULT16_16_Q15(g1, t0);
405    g11 = MULT16_16_Q15(g1, t1);
406    g12 = MULT16_16_Q15(g1, t2);
407    for (i=0;i<overlap;i++)
408    {
409       celt_word16 f;
410       f = MULT16_16_Q15(window[i],window[i]);
411       y[i] = x[i]
412                + MULT16_32_Q15(MULT16_16_Q15((Q15ONE-f),g01),x[i-T0])
413                + MULT16_32_Q15(MULT16_16_Q15((Q15ONE-f),g00),x[i-T0-1])
414                + MULT16_32_Q15(MULT16_16_Q15((Q15ONE-f),g02),x[i-T0+1])
415                + MULT16_32_Q15(MULT16_16_Q15(f,g11),x[i-T1])
416                + MULT16_32_Q15(MULT16_16_Q15(f,g10),x[i-T1-1])
417                + MULT16_32_Q15(MULT16_16_Q15(f,g12),x[i-T1+1]);
418
419    }
420    for (i=overlap;i<N;i++)
421       y[i] = x[i]
422                + MULT16_32_Q15(g11,x[i-T1])
423                + MULT16_32_Q15(g10,x[i-T1-1])
424                + MULT16_32_Q15(g12,x[i-T1+1]);
425 }
426 #endif /* ENABLE_POSTFILTER */
427
428 static const signed char tf_select_table[4][8] = {
429       {0, -1, 0, -1,    0,-1, 0,-1},
430       {0, -1, 0, -2,    1, 0, 1 -1},
431       {0, -2, 0, -3,    2, 0, 1 -1},
432       {0, -2, 0, -3,    2, 0, 1 -1},
433 };
434
435 static celt_word32 l1_metric(const celt_norm *tmp, int N, int LM, int width)
436 {
437    int i, j;
438    static const celt_word16 sqrtM_1[4] = {Q15ONE, QCONST16(0.70711f,15), QCONST16(0.5f,15), QCONST16(0.35355f,15)};
439    celt_word32 L1;
440    celt_word16 bias;
441    L1=0;
442    for (i=0;i<1<<LM;i++)
443    {
444       celt_word32 L2 = 0;
445       for (j=0;j<N>>LM;j++)
446          L2 = MAC16_16(L2, tmp[(j<<LM)+i], tmp[(j<<LM)+i]);
447       L1 += celt_sqrt(L2);
448    }
449    L1 = MULT16_32_Q15(sqrtM_1[LM], L1);
450    if (width==1)
451       bias = QCONST16(.12f,15)*LM;
452    else if (width==2)
453       bias = QCONST16(.05f,15)*LM;
454    else
455       bias = QCONST16(.02f,15)*LM;
456    L1 = MAC16_32_Q15(L1, bias, L1);
457    return L1;
458 }
459
460 static int tf_analysis(const CELTMode *m, celt_word16 *bandLogE, celt_word16 *oldBandE,
461       int len, int C, int isTransient, int *tf_res, int nbCompressedBytes, celt_norm *X,
462       int N0, int LM, int *tf_sum)
463 {
464    int i;
465    VARDECL(int, metric);
466    int cost0;
467    int cost1;
468    VARDECL(int, path0);
469    VARDECL(int, path1);
470    VARDECL(celt_norm, tmp);
471    int lambda;
472    int tf_select=0;
473    SAVE_STACK;
474
475    /* FIXME: Should check number of bytes *left* */
476    if (nbCompressedBytes<15*C)
477    {
478       *tf_sum = 0;
479       for (i=0;i<len;i++)
480          tf_res[i] = isTransient;
481       return 0;
482    }
483    if (nbCompressedBytes<40)
484       lambda = 12;
485    else if (nbCompressedBytes<60)
486       lambda = 6;
487    else if (nbCompressedBytes<100)
488       lambda = 4;
489    else
490       lambda = 3;
491
492    ALLOC(metric, len, int);
493    ALLOC(tmp, (m->eBands[len]-m->eBands[len-1])<<LM, celt_norm);
494    ALLOC(path0, len, int);
495    ALLOC(path1, len, int);
496
497    *tf_sum = 0;
498    for (i=0;i<len;i++)
499    {
500       int j, k, N;
501       celt_word32 L1, best_L1;
502       int best_level=0;
503       N = (m->eBands[i+1]-m->eBands[i])<<LM;
504       for (j=0;j<N;j++)
505          tmp[j] = X[j+(m->eBands[i]<<LM)];
506       /* FIXME: Do something with the right channel */
507       if (C==2)
508          for (j=0;j<N;j++)
509             tmp[j] = ADD16(tmp[j],X[N0+j+(m->eBands[i]<<LM)]);
510       L1 = l1_metric(tmp, N, isTransient ? LM : 0, N>>LM);
511       best_L1 = L1;
512       /*printf ("%f ", L1);*/
513       for (k=0;k<LM;k++)
514       {
515          int B;
516
517          if (isTransient)
518             B = (LM-k-1);
519          else
520             B = k+1;
521
522          if (isTransient)
523             haar1(tmp, N>>(LM-k), 1<<(LM-k));
524          else
525             haar1(tmp, N>>k, 1<<k);
526
527          L1 = l1_metric(tmp, N, B, N>>LM);
528
529          if (L1 < best_L1)
530          {
531             best_L1 = L1;
532             best_level = k+1;
533          }
534       }
535       /*printf ("%d ", isTransient ? LM-best_level : best_level);*/
536       if (isTransient)
537          metric[i] = best_level;
538       else
539          metric[i] = -best_level;
540       *tf_sum += metric[i];
541    }
542    /*printf("\n");*/
543    /* FIXME: Figure out how to set this */
544    tf_select = 0;
545
546    cost0 = 0;
547    cost1 = isTransient ? 0 : lambda;
548    /* Viterbi forward pass */
549    for (i=1;i<len;i++)
550    {
551       int curr0, curr1;
552       int from0, from1;
553
554       from0 = cost0;
555       from1 = cost1 + lambda;
556       if (from0 < from1)
557       {
558          curr0 = from0;
559          path0[i]= 0;
560       } else {
561          curr0 = from1;
562          path0[i]= 1;
563       }
564
565       from0 = cost0 + lambda;
566       from1 = cost1;
567       if (from0 < from1)
568       {
569          curr1 = from0;
570          path1[i]= 0;
571       } else {
572          curr1 = from1;
573          path1[i]= 1;
574       }
575       cost0 = curr0 + abs(metric[i]-tf_select_table[LM][4*isTransient+2*tf_select+0]);
576       cost1 = curr1 + abs(metric[i]-tf_select_table[LM][4*isTransient+2*tf_select+1]);
577    }
578    tf_res[len-1] = cost0 < cost1 ? 0 : 1;
579    /* Viterbi backward pass to check the decisions */
580    for (i=len-2;i>=0;i--)
581    {
582       if (tf_res[i+1] == 1)
583          tf_res[i] = path1[i+1];
584       else
585          tf_res[i] = path0[i+1];
586    }
587    RESTORE_STACK;
588    return tf_select;
589 }
590
591 static void tf_encode(int start, int end, int isTransient, int *tf_res, int LM, int tf_select, ec_enc *enc)
592 {
593    int curr, i;
594    ec_enc_bit_prob(enc, tf_res[start], isTransient ? 16384 : 4096);
595    curr = tf_res[start];
596    for (i=start+1;i<end;i++)
597    {
598       ec_enc_bit_prob(enc, tf_res[i] ^ curr, isTransient ? 4096 : 2048);
599       curr = tf_res[i];
600    }
601    if (LM!=0)
602       ec_enc_bits(enc, tf_select, 1);
603    for (i=start;i<end;i++)
604       tf_res[i] = tf_select_table[LM][4*isTransient+2*tf_select+tf_res[i]];
605    /*printf("%d %d ", isTransient, tf_select); for(i=0;i<end;i++)printf("%d ", tf_res[i]);printf("\n");*/
606 }
607
608 static void tf_decode(int start, int end, int C, int isTransient, int *tf_res, int LM, ec_dec *dec)
609 {
610    int i, curr, tf_select;
611    tf_res[start] = ec_dec_bit_prob(dec, isTransient ? 16384 : 4096);
612    curr = tf_res[start];
613    for (i=start+1;i<end;i++)
614    {
615       tf_res[i] = ec_dec_bit_prob(dec, isTransient ? 4096 : 2048) ^ curr;
616       curr = tf_res[i];
617    }
618    if (LM!=0)
619       tf_select = ec_dec_bits(dec, 1);
620    else
621       tf_select = 0;
622    for (i=start;i<end;i++)
623       tf_res[i] = tf_select_table[LM][4*isTransient+2*tf_select+tf_res[i]];
624 }
625
626 static int alloc_trim_analysis(const CELTMode *m, const celt_norm *X,
627       const celt_word16 *bandLogE, int nbEBands, int LM, int C, int N0)
628 {
629    int i;
630    celt_word32 diff=0;
631    int c;
632    int trim_index = 5;
633    if (C==2)
634    {
635       celt_word16 sum = 0; /* Q10 */
636       /* Compute inter-channel correlation for low frequencies */
637       for (i=0;i<8;i++)
638       {
639          int j;
640          celt_word32 partial = 0;
641          for (j=m->eBands[i]<<LM;j<m->eBands[i+1]<<LM;j++)
642             partial = MAC16_16(partial, X[j], X[N0+j]);
643          sum = ADD16(sum, EXTRACT16(SHR32(partial, 18)));
644       }
645       sum = MULT16_16_Q15(QCONST16(1.f/8, 15), sum);
646       /*printf ("%f\n", sum);*/
647       if (sum > QCONST16(.995f,10))
648          trim_index-=4;
649       else if (sum > QCONST16(.92f,10))
650          trim_index-=3;
651       else if (sum > QCONST16(.85f,10))
652          trim_index-=2;
653       else if (sum > QCONST16(.8f,10))
654          trim_index-=1;
655    }
656
657    /* Estimate spectral tilt */
658    c=0; do {
659       for (i=0;i<nbEBands-1;i++)
660       {
661          diff += bandLogE[i+c*nbEBands]*(celt_int32)(2+2*i-nbEBands);
662       }
663    } while (++c<0);
664    diff /= C*(nbEBands-1);
665    /*printf("%f\n", diff);*/
666    if (diff > QCONST16(2.f, DB_SHIFT))
667       trim_index--;
668    if (diff > QCONST16(8.f, DB_SHIFT))
669       trim_index--;
670    if (diff < -QCONST16(4.f, DB_SHIFT))
671       trim_index++;
672    if (diff < -QCONST16(10.f, DB_SHIFT))
673       trim_index++;
674
675    if (trim_index<0)
676       trim_index = 0;
677    if (trim_index>10)
678       trim_index = 10;
679    return trim_index;
680 }
681
682 static int stereo_analysis(const CELTMode *m, const celt_norm *X,
683       int nbEBands, int LM, int C, int N0)
684 {
685    int i;
686    int thetas;
687    celt_word32 sumLR = EPSILON, sumMS = EPSILON;
688
689    /* Use the L1 norm to model the entropy of the L/R signal vs the M/S signal */
690    for (i=0;i<13;i++)
691    {
692       int j;
693       for (j=m->eBands[i]<<LM;j<m->eBands[i+1]<<LM;j++)
694       {
695          celt_word16 L, R, M, S;
696          L = X[j];
697          R = X[N0+j];
698          M = L+R;
699          S = L-R;
700          sumLR += EXTEND32(ABS16(L)) + EXTEND32(ABS16(R));
701          sumMS += EXTEND32(ABS16(M)) + EXTEND32(ABS16(S));
702       }
703    }
704    sumMS = MULT16_32_Q15(QCONST16(0.707107f, 15), sumMS);
705    thetas = 13;
706    /* We don't need thetas for lower bands with LM<=1 */
707    if (LM<=1)
708       thetas -= 8;
709    return MULT16_32_Q15((m->eBands[13]<<(LM+1))+thetas, sumMS)
710          > MULT16_32_Q15(m->eBands[13]<<(LM+1), sumLR);
711 }
712
713 #ifdef FIXED_POINT
714 int celt_encode_with_ec(CELTEncoder * restrict st, const celt_int16 * pcm, int frame_size, unsigned char *compressed, int nbCompressedBytes, ec_enc *enc)
715 {
716 #else
717 int celt_encode_with_ec_float(CELTEncoder * restrict st, const celt_sig * pcm, int frame_size, unsigned char *compressed, int nbCompressedBytes, ec_enc *enc)
718 {
719 #endif
720    int i, c, N;
721    int bits;
722    ec_byte_buffer buf;
723    ec_enc         _enc;
724    VARDECL(celt_sig, in);
725    VARDECL(celt_sig, freq);
726    VARDECL(celt_norm, X);
727    VARDECL(celt_ener, bandE);
728    VARDECL(celt_word16, bandLogE);
729    VARDECL(int, fine_quant);
730    VARDECL(celt_word16, error);
731    VARDECL(int, pulses);
732    VARDECL(int, offsets);
733    VARDECL(int, fine_priority);
734    VARDECL(int, tf_res);
735    celt_sig *_overlap_mem;
736    celt_sig *prefilter_mem;
737    celt_word16 *oldBandE;
738    int shortBlocks=0;
739    int isTransient=0;
740    int resynth;
741    const int C = CHANNELS(st->channels);
742    int LM, M;
743    int tf_select;
744    int nbFilledBytes, nbAvailableBytes;
745    int effEnd;
746    int codedBands;
747    int tf_sum;
748    int alloc_trim;
749    int pitch_index=COMBFILTER_MINPERIOD;
750    celt_word16 gain1 = 0;
751    int intensity=0;
752    int dual_stereo=0;
753    int effectiveBytes;
754    SAVE_STACK;
755
756    if (nbCompressedBytes<0 || pcm==NULL)
757      return CELT_BAD_ARG;
758
759    for (LM=0;LM<4;LM++)
760       if (st->mode->shortMdctSize<<LM==frame_size)
761          break;
762    if (LM>=MAX_CONFIG_SIZES)
763       return CELT_BAD_ARG;
764    M=1<<LM;
765
766    prefilter_mem = st->in_mem+C*(st->overlap);
767    _overlap_mem = prefilter_mem+C*COMBFILTER_MAXPERIOD;
768    /*_overlap_mem = st->in_mem+C*(st->overlap);*/
769    oldBandE = (celt_word16*)(st->in_mem+C*(2*st->overlap+COMBFILTER_MAXPERIOD));
770
771    if (enc==NULL)
772    {
773       ec_byte_writeinit_buffer(&buf, compressed, nbCompressedBytes);
774       ec_enc_init(&_enc,&buf);
775       enc = &_enc;
776       nbFilledBytes=0;
777    } else {
778       nbFilledBytes=(ec_enc_tell(enc, 0)+4)>>3;
779    }
780    nbAvailableBytes = nbCompressedBytes - nbFilledBytes;
781
782    if (st->vbr_rate_norm>0)
783       effectiveBytes = st->vbr_rate_norm>>BITRES<<LM>>3;
784    else
785       effectiveBytes = nbCompressedBytes;
786
787    effEnd = st->end;
788    if (effEnd > st->mode->effEBands)
789       effEnd = st->mode->effEBands;
790
791    N = M*st->mode->shortMdctSize;
792    ALLOC(in, C*(N+st->overlap), celt_sig);
793
794    /* Find pitch period and gain */
795    {
796       VARDECL(celt_sig, _pre);
797       celt_sig *pre[2];
798       SAVE_STACK;
799       c = 0;
800       ALLOC(_pre, C*(N+COMBFILTER_MAXPERIOD), celt_sig);
801
802       pre[0] = _pre;
803       pre[1] = _pre + (N+COMBFILTER_MAXPERIOD);
804
805       c=0; do {
806          const celt_word16 * restrict pcmp = pcm+c;
807          celt_sig * restrict inp = in+c*(N+st->overlap)+st->overlap;
808
809          for (i=0;i<N;i++)
810          {
811             /* Apply pre-emphasis */
812             celt_sig tmp = MULT16_16(st->mode->preemph[2], SCALEIN(*pcmp));
813             *inp = tmp + st->preemph_memE[c];
814             st->preemph_memE[c] = MULT16_32_Q15(st->mode->preemph[1], *inp)
815                                    - MULT16_32_Q15(st->mode->preemph[0], tmp);
816             inp++;
817             pcmp+=C;
818          }
819          CELT_COPY(pre[c], prefilter_mem+c*COMBFILTER_MAXPERIOD, COMBFILTER_MAXPERIOD);
820          CELT_COPY(pre[c]+COMBFILTER_MAXPERIOD, in+c*(N+st->overlap)+st->overlap, N);
821       } while (++c<C);
822
823 #ifdef ENABLE_POSTFILTER
824       if (nbAvailableBytes>12*C)
825       {
826          VARDECL(celt_word16, pitch_buf);
827          ALLOC(pitch_buf, (COMBFILTER_MAXPERIOD+N)>>1, celt_word16);
828          celt_word32 tmp=0;
829          celt_word32 mem0[2]={0,0};
830          celt_word16 mem1[2]={0,0};
831
832          pitch_downsample(pre, pitch_buf, COMBFILTER_MAXPERIOD+N, COMBFILTER_MAXPERIOD+N,
833                           C, mem0, mem1);
834          pitch_search(st->mode, pitch_buf+(COMBFILTER_MAXPERIOD>>1), pitch_buf, N,
835                COMBFILTER_MAXPERIOD-COMBFILTER_MINPERIOD, &pitch_index, &tmp, 1<<LM);
836          pitch_index = COMBFILTER_MAXPERIOD-pitch_index;
837
838          gain1 = remove_doubling(pitch_buf, COMBFILTER_MAXPERIOD, COMBFILTER_MINPERIOD,
839                N, &pitch_index, st->prefilter_period, st->prefilter_gain);
840          if (pitch_index > COMBFILTER_MAXPERIOD)
841             pitch_index = COMBFILTER_MAXPERIOD;
842          gain1 = MULT16_16_Q15(QCONST16(.7f,15),gain1);
843          if (gain1 > QCONST16(.6f,15))
844             gain1 = QCONST16(.6f,15);
845          if (ABS16(gain1-st->prefilter_gain)<QCONST16(.1,15))
846             gain1=st->prefilter_gain;
847       } else {
848          gain1 = 0;
849       }
850       if (gain1<QCONST16(.2f,15) || (nbAvailableBytes<30 && gain1<QCONST16(.4f,15)))
851       {
852          ec_enc_bit_prob(enc, 0, 32768);
853          gain1 = 0;
854       } else {
855          int qg;
856          int octave;
857 #ifdef FIXED_POINT
858          qg = ((gain1+2048)>>12)-2;
859 #else
860          qg = floor(.5+gain1*8)-2;
861 #endif
862          ec_enc_bit_prob(enc, 1, 32768);
863          octave = EC_ILOG(pitch_index)-5;
864          ec_enc_uint(enc, octave, 6);
865          ec_enc_bits(enc, pitch_index-(16<<octave), 4+octave);
866          ec_enc_bits(enc, qg, 2);
867          gain1 = QCONST16(.125f,15)*(qg+2);
868       }
869       /*printf("%d %f\n", pitch_index, gain1);*/
870 #else /* ENABLE_POSTFILTER */
871       ec_enc_bit_prob(enc, 0, 32768);
872 #endif /* ENABLE_POSTFILTER */
873
874       c=0; do {
875          st->prefilter_period=IMAX(st->prefilter_period, COMBFILTER_MINPERIOD);
876          CELT_COPY(in+c*(N+st->overlap), st->in_mem+c*(st->overlap), st->overlap);
877 #ifdef ENABLE_POSTFILTER
878          comb_filter(in+c*(N+st->overlap)+st->overlap, pre[c]+COMBFILTER_MAXPERIOD,
879                st->prefilter_period, pitch_index, N, C, -st->prefilter_gain, -gain1, st->mode->window, st->mode->overlap);
880 #endif /* ENABLE_POSTFILTER */
881          CELT_COPY(st->in_mem+c*(st->overlap), in+c*(N+st->overlap)+N, st->overlap);
882
883 #ifdef ENABLE_POSTFILTER
884          if (N>COMBFILTER_MAXPERIOD)
885          {
886             CELT_MOVE(prefilter_mem+c*COMBFILTER_MAXPERIOD, pre[c]+N, COMBFILTER_MAXPERIOD);
887          } else {
888             CELT_MOVE(prefilter_mem+c*COMBFILTER_MAXPERIOD, prefilter_mem+c*COMBFILTER_MAXPERIOD+N, COMBFILTER_MAXPERIOD-N);
889             CELT_MOVE(prefilter_mem+c*COMBFILTER_MAXPERIOD+COMBFILTER_MAXPERIOD-N, pre[c]+COMBFILTER_MAXPERIOD, N);
890          }
891 #endif /* ENABLE_POSTFILTER */
892       } while (++c<C);
893
894       RESTORE_STACK;
895    }
896
897 #ifdef RESYNTH
898    resynth = 1;
899 #else
900    resynth = 0;
901 #endif
902
903    if (st->complexity > 1 && LM>0)
904    {
905       isTransient = M > 1 &&
906          transient_analysis(in, N+st->overlap, C, &st->frame_max, st->overlap);
907    } else {
908       isTransient = 0;
909    }
910
911    if (isTransient)
912       shortBlocks = M;
913    else
914       shortBlocks = 0;
915
916    ALLOC(freq, C*N, celt_sig); /**< Interleaved signal MDCTs */
917    ALLOC(bandE,st->mode->nbEBands*C, celt_ener);
918    ALLOC(bandLogE,st->mode->nbEBands*C, celt_word16);
919    /* Compute MDCTs */
920    compute_mdcts(st->mode, shortBlocks, in, freq, C, LM);
921
922    ALLOC(X, C*N, celt_norm);         /**< Interleaved normalised MDCTs */
923
924    compute_band_energies(st->mode, freq, bandE, effEnd, C, M);
925
926    amp2Log2(st->mode, effEnd, st->end, bandE, bandLogE, C);
927
928    /* Band normalisation */
929    normalise_bands(st->mode, freq, X, bandE, effEnd, C, M);
930
931    ALLOC(tf_res, st->mode->nbEBands, int);
932    /* Needs to be before coarse energy quantization because otherwise the energy gets modified */
933    tf_select = tf_analysis(st->mode, bandLogE, oldBandE, effEnd, C, isTransient, tf_res, effectiveBytes, X, N, LM, &tf_sum);
934    for (i=effEnd;i<st->end;i++)
935       tf_res[i] = tf_res[effEnd-1];
936
937    ALLOC(error, C*st->mode->nbEBands, celt_word16);
938    quant_coarse_energy(st->mode, st->start, st->end, effEnd, bandLogE,
939          oldBandE, nbCompressedBytes*8, error, enc,
940          C, LM, nbAvailableBytes, st->force_intra,
941          &st->delayedIntra, st->complexity >= 4);
942
943    if (LM > 0)
944       ec_enc_bit_prob(enc, shortBlocks!=0, 8192);
945
946    tf_encode(st->start, st->end, isTransient, tf_res, LM, tf_select, enc);
947
948    if (shortBlocks || st->complexity < 3 || nbAvailableBytes < 10*C)
949    {
950       if (st->complexity == 0)
951       {
952          st->spread_decision = SPREAD_NONE;
953       } else {
954          st->spread_decision = SPREAD_NORMAL;
955       }
956    } else {
957       st->spread_decision = spreading_decision(st->mode, X, &st->tonal_average, st->spread_decision, effEnd, C, M);
958    }
959    /* Probs: NONE: 21.875%, LIGHT: 6.25%, NORMAL: 65.625%, AGGRESSIVE: 6.25% */
960    ec_enc_bit_prob(enc, st->spread_decision>>1, 18432);
961    ec_enc_bit_prob(enc, st->spread_decision&1,
962          (st->spread_decision>>1) ? 5699 : 14564);
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_bit_prob(dec, 18432)<<1;
1827    spread_decision |= ec_dec_bit_prob(dec, (spread_decision>>1) ? 5699 : 14564);
1828
1829    ALLOC(pulses, st->mode->nbEBands, int);
1830    ALLOC(offsets, st->mode->nbEBands, int);
1831    ALLOC(fine_priority, st->mode->nbEBands, int);
1832
1833    for (i=0;i<st->mode->nbEBands;i++)
1834       offsets[i] = 0;
1835    for (i=0;i<st->mode->nbEBands;i++)
1836    {
1837       if (ec_dec_bit_prob(dec, 1024))
1838       {
1839          while (ec_dec_bit_prob(dec, 32768))
1840             offsets[i]++;
1841          offsets[i]++;
1842          offsets[i] *= (6<<BITRES);
1843       }
1844    }
1845
1846    ALLOC(fine_quant, st->mode->nbEBands, int);
1847    {
1848       int fl;
1849       alloc_trim = 0;
1850       fl = ec_decode_bin(dec, 7);
1851       while (trim_cdf[alloc_trim+1] <= fl)
1852          alloc_trim++;
1853       ec_dec_update(dec, trim_cdf[alloc_trim], trim_cdf[alloc_trim+1], 128);
1854    }
1855
1856    if (C==2)
1857    {
1858       dual_stereo = ec_dec_bit_prob(dec, 32768);
1859       intensity = ec_dec_uint(dec, 1+st->end-st->start);
1860    }
1861
1862    bits = (len*8<<BITRES) - ec_dec_tell(dec, BITRES) - 1;
1863    codedBands = compute_allocation(st->mode, st->start, st->end, offsets,
1864          alloc_trim, bits, pulses, fine_quant, fine_priority, C, LM, dec, 0, 0);
1865    
1866    unquant_fine_energy(st->mode, st->start, st->end, bandE, oldBandE, fine_quant, dec, C);
1867
1868    /* Decode fixed codebook */
1869    quant_all_bands(0, st->mode, st->start, st->end, X, C==2 ? X+N : NULL,
1870          NULL, pulses, shortBlocks, spread_decision, dual_stereo, intensity, tf_res, 1,
1871          len*8, dec, LM, codedBands);
1872
1873    unquant_energy_finalise(st->mode, st->start, st->end, bandE, oldBandE,
1874          fine_quant, fine_priority, len*8-ec_dec_tell(dec, 0), dec, C);
1875
1876    log2Amp(st->mode, st->start, st->end, bandE, oldBandE, C);
1877
1878    /* Synthesis */
1879    denormalise_bands(st->mode, X, freq, bandE, effEnd, C, M);
1880
1881    CELT_MOVE(decode_mem[0], decode_mem[0]+N, DECODE_BUFFER_SIZE-N);
1882    if (C==2)
1883       CELT_MOVE(decode_mem[1], decode_mem[1]+N, DECODE_BUFFER_SIZE-N);
1884
1885    c=0; do
1886       for (i=0;i<M*st->mode->eBands[st->start];i++)
1887          freq[c*N+i] = 0;
1888    while (++c<C);
1889    c=0; do
1890       for (i=M*st->mode->eBands[effEnd];i<N;i++)
1891          freq[c*N+i] = 0;
1892    while (++c<C);
1893
1894    out_syn[0] = out_mem[0]+MAX_PERIOD-N;
1895    if (C==2)
1896       out_syn[1] = out_mem[1]+MAX_PERIOD-N;
1897
1898    /* Compute inverse MDCTs */
1899    compute_inv_mdcts(st->mode, shortBlocks, freq, out_syn, overlap_mem, C, LM);
1900
1901 #ifdef ENABLE_POSTFILTER
1902    c=0; do {
1903       st->postfilter_period=IMAX(st->postfilter_period, COMBFILTER_MINPERIOD);
1904       st->postfilter_period_old=IMAX(st->postfilter_period_old, COMBFILTER_MINPERIOD);
1905       if (LM!=0)
1906       {
1907          comb_filter(out_syn[c], out_syn[c], st->postfilter_period, st->postfilter_period, st->overlap, C,
1908                st->postfilter_gain, st->postfilter_gain, NULL, 0);
1909          comb_filter(out_syn[c]+st->overlap, out_syn[c]+st->overlap, st->postfilter_period, postfilter_pitch, N-st->overlap, C,
1910                st->postfilter_gain, postfilter_gain, st->mode->window, st->mode->overlap);
1911       } else {
1912          comb_filter(out_syn[c], out_syn[c], st->postfilter_period_old, st->postfilter_period, N-st->overlap, C,
1913                st->postfilter_gain_old, st->postfilter_gain, st->mode->window, st->mode->overlap);
1914       }
1915    } while (++c<C);
1916    st->postfilter_period_old = st->postfilter_period;
1917    st->postfilter_gain_old = st->postfilter_gain;
1918    st->postfilter_period = postfilter_pitch;
1919    st->postfilter_gain = postfilter_gain;
1920 #endif /* ENABLE_POSTFILTER */
1921
1922    deemphasis(out_syn, pcm, N, C, st->mode->preemph, st->preemph_memD);
1923    st->loss_count = 0;
1924    RESTORE_STACK;
1925    if (ec_dec_get_error(dec))
1926       return CELT_CORRUPTED_DATA;
1927    else
1928       return CELT_OK;
1929 }
1930
1931 #ifdef FIXED_POINT
1932 #ifndef DISABLE_FLOAT_API
1933 int celt_decode_with_ec_float(CELTDecoder * restrict st, const unsigned char *data, int len, float * restrict pcm, int frame_size, ec_dec *dec)
1934 {
1935    int j, ret, C, N, LM, M;
1936    VARDECL(celt_int16, out);
1937    SAVE_STACK;
1938
1939    if (pcm==NULL)
1940       return CELT_BAD_ARG;
1941
1942    for (LM=0;LM<4;LM++)
1943       if (st->mode->shortMdctSize<<LM==frame_size)
1944          break;
1945    if (LM>=MAX_CONFIG_SIZES)
1946       return CELT_BAD_ARG;
1947    M=1<<LM;
1948
1949    C = CHANNELS(st->channels);
1950    N = M*st->mode->shortMdctSize;
1951    
1952    ALLOC(out, C*N, celt_int16);
1953    ret=celt_decode_with_ec(st, data, len, out, frame_size, dec);
1954    if (ret==0)
1955       for (j=0;j<C*N;j++)
1956          pcm[j]=out[j]*(1.f/32768.f);
1957      
1958    RESTORE_STACK;
1959    return ret;
1960 }
1961 #endif /*DISABLE_FLOAT_API*/
1962 #else
1963 int celt_decode_with_ec(CELTDecoder * restrict st, const unsigned char *data, int len, celt_int16 * restrict pcm, int frame_size, ec_dec *dec)
1964 {
1965    int j, ret, C, N, LM, M;
1966    VARDECL(celt_sig, out);
1967    SAVE_STACK;
1968
1969    if (pcm==NULL)
1970       return CELT_BAD_ARG;
1971
1972    for (LM=0;LM<4;LM++)
1973       if (st->mode->shortMdctSize<<LM==frame_size)
1974          break;
1975    if (LM>=MAX_CONFIG_SIZES)
1976       return CELT_BAD_ARG;
1977    M=1<<LM;
1978
1979    C = CHANNELS(st->channels);
1980    N = M*st->mode->shortMdctSize;
1981    ALLOC(out, C*N, celt_sig);
1982
1983    ret=celt_decode_with_ec_float(st, data, len, out, frame_size, dec);
1984
1985    if (ret==0)
1986       for (j=0;j<C*N;j++)
1987          pcm[j] = FLOAT2INT16 (out[j]);
1988    
1989    RESTORE_STACK;
1990    return ret;
1991 }
1992 #endif
1993
1994 int celt_decode(CELTDecoder * restrict st, const unsigned char *data, int len, celt_int16 * restrict pcm, int frame_size)
1995 {
1996    return celt_decode_with_ec(st, data, len, pcm, frame_size, NULL);
1997 }
1998
1999 #ifndef DISABLE_FLOAT_API
2000 int celt_decode_float(CELTDecoder * restrict st, const unsigned char *data, int len, float * restrict pcm, int frame_size)
2001 {
2002    return celt_decode_with_ec_float(st, data, len, pcm, frame_size, NULL);
2003 }
2004 #endif /* DISABLE_FLOAT_API */
2005
2006 int celt_decoder_ctl(CELTDecoder * restrict st, int request, ...)
2007 {
2008    va_list ap;
2009
2010    va_start(ap, request);
2011    switch (request)
2012    {
2013       case CELT_GET_MODE_REQUEST:
2014       {
2015          const CELTMode ** value = va_arg(ap, const CELTMode**);
2016          if (value==0)
2017             goto bad_arg;
2018          *value=st->mode;
2019       }
2020       break;
2021       case CELT_SET_START_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->start = value;
2027       }
2028       break;
2029       case CELT_SET_END_BAND_REQUEST:
2030       {
2031          celt_int32 value = va_arg(ap, celt_int32);
2032          if (value<0 || value>=st->mode->nbEBands)
2033             goto bad_arg;
2034          st->end = value;
2035       }
2036       break;
2037       case CELT_RESET_STATE:
2038       {
2039          CELT_MEMSET((char*)&st->DECODER_RESET_START, 0,
2040                celt_decoder_get_size(st->mode, st->channels)-
2041                ((char*)&st->DECODER_RESET_START - (char*)st));
2042       }
2043       break;
2044       default:
2045          goto bad_request;
2046    }
2047    va_end(ap);
2048    return CELT_OK;
2049 bad_arg:
2050    va_end(ap);
2051    return CELT_BAD_ARG;
2052 bad_request:
2053       va_end(ap);
2054   return CELT_UNIMPLEMENTED;
2055 }
2056
2057 const char *celt_strerror(int error)
2058 {
2059    static const char *error_strings[8] = {
2060       "success",
2061       "invalid argument",
2062       "invalid mode",
2063       "internal error",
2064       "corrupted stream",
2065       "request not implemented",
2066       "invalid state",
2067       "memory allocation failed"
2068    };
2069    if (error > 0 || error < -7)
2070       return "unknown error";
2071    else 
2072       return error_strings[-error];
2073 }
2074