26849ac28142a365fa35dbc01dda20a406c74b68
[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 fold_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->fold_decision = 1;
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    int has_fold=1;
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 (LM != 0 && nbAvailableBytes>10*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          has_fold = 0;
954          st->fold_decision = 3;
955       } else {
956          has_fold = 1;
957          st->fold_decision = 1;
958       }
959    } else {
960       has_fold = folding_decision(st->mode, X, &st->tonal_average, &st->fold_decision, effEnd, C, M);
961    }
962    ec_enc_bit_prob(enc, has_fold>>1, 8192);
963    ec_enc_bit_prob(enc, has_fold&1, (has_fold>>1) ? 32768 : 49152);
964
965    ALLOC(offsets, st->mode->nbEBands, int);
966
967    for (i=0;i<st->mode->nbEBands;i++)
968       offsets[i] = 0;
969    /* Dynamic allocation code */
970    /* Make sure that dynamic allocation can't make us bust the budget */
971    if (effectiveBytes > 50 && LM>=1)
972    {
973       int t1, t2;
974       if (LM <= 1)
975       {
976          t1 = 3;
977          t2 = 5;
978       } else {
979          t1 = 2;
980          t2 = 4;
981       }
982       for (i=1;i<st->mode->nbEBands-1;i++)
983       {
984          celt_word32 d2;
985          d2 = 2*bandLogE[i]-bandLogE[i-1]-bandLogE[i+1];
986          if (C==2)
987             d2 = HALF32(d2 + 2*bandLogE[i+st->mode->nbEBands]-
988                   bandLogE[i-1+st->mode->nbEBands]-bandLogE[i+1+st->mode->nbEBands]);
989          if (d2 > SHL16(t1,DB_SHIFT))
990             offsets[i] += 1;
991          if (d2 > SHL16(t2,DB_SHIFT))
992             offsets[i] += 1;
993       }
994    }
995    for (i=0;i<st->mode->nbEBands;i++)
996    {
997       int j;
998       ec_enc_bit_prob(enc, offsets[i]!=0, 1024);
999       if (offsets[i]!=0)
1000       {
1001          for (j=0;j<offsets[i]-1;j++)
1002             ec_enc_bit_prob(enc, 1, 32768);
1003          ec_enc_bit_prob(enc, 0, 32768);
1004       }
1005       offsets[i] *= (6<<BITRES);
1006    }
1007    alloc_trim = alloc_trim_analysis(st->mode, X, bandLogE, st->mode->nbEBands, LM, C, N);
1008    ec_encode_bin(enc, trim_cdf[alloc_trim], trim_cdf[alloc_trim+1], 7);
1009
1010    /* Variable bitrate */
1011    if (st->vbr_rate_norm>0)
1012    {
1013      celt_word16 alpha;
1014      celt_int32 delta, tell;
1015      /* The target rate in 8th bits per frame */
1016      celt_int32 vbr_rate;
1017      celt_int32 target;
1018      celt_int32 vbr_bound, max_allowed, min_allowed;
1019
1020      target = vbr_rate = M*st->vbr_rate_norm;
1021
1022      target = target + st->vbr_offset - ((40*C+20)<<BITRES);
1023
1024      /* Shortblocks get a large boost in bitrate, but since they
1025         are uncommon long blocks are not greatly affected */
1026      if (shortBlocks || tf_sum < -2*(st->end-st->start))
1027         target = 7*target/4;
1028      else if (tf_sum < -(st->end-st->start))
1029         target = 3*target/2;
1030      else if (M > 1)
1031         target-=(target+14)/28;
1032
1033      tell = ec_enc_tell(enc, BITRES);
1034
1035      /* The current offset is removed from the target and the space used
1036         so far is added*/
1037      target=target+tell;
1038      /* By how much did we "miss" the target on that frame */
1039      delta = target - vbr_rate;
1040
1041      /* Computes the max bit-rate allowed in VBR more to avoid violating the target rate and buffering */
1042      vbr_bound = vbr_rate;
1043      if (st->constrained_vbr)
1044         max_allowed = IMIN(vbr_rate+vbr_bound-st->vbr_reservoir>>(BITRES+3),nbAvailableBytes);
1045      else
1046         max_allowed = nbAvailableBytes;
1047      min_allowed = (tell>>(BITRES+3)) + 2 - nbFilledBytes;
1048
1049      /* In VBR mode the frame size must not be reduced so much that it would result in the encoder running out of bits */
1050      nbAvailableBytes = target+(1<<(BITRES+2))>>(BITRES+3);
1051      nbAvailableBytes=IMAX(min_allowed,IMIN(max_allowed,nbAvailableBytes));
1052      target=nbAvailableBytes<<(BITRES+3);
1053
1054      if (st->vbr_count < 970)
1055      {
1056         st->vbr_count++;
1057         alpha = celt_rcp(SHL32(EXTEND32(st->vbr_count+20),16));
1058      } else
1059         alpha = QCONST16(.001f,15);
1060      /* How many bits have we used in excess of what we're allowed */
1061      if (st->constrained_vbr)
1062         st->vbr_reservoir += target - vbr_rate;
1063      /*printf ("%d\n", st->vbr_reservoir);*/
1064
1065      /* Compute the offset we need to apply in order to reach the target */
1066      st->vbr_drift += (celt_int32)MULT16_32_Q15(alpha,delta-st->vbr_offset-st->vbr_drift);
1067      st->vbr_offset = -st->vbr_drift;
1068      /*printf ("%d\n", st->vbr_drift);*/
1069
1070      /* We could use any multiple of vbr_rate as bound (depending on the delay) */
1071      if (st->constrained_vbr && st->vbr_reservoir < 0)
1072      {
1073         /* We're under the min value -- increase rate */
1074         int adjust = (-st->vbr_reservoir)/(8<<BITRES);
1075         nbAvailableBytes += adjust;
1076         st->vbr_reservoir = 0;
1077         /*printf ("+%d\n", adjust);*/
1078      }
1079      nbCompressedBytes = IMIN(nbCompressedBytes,nbAvailableBytes+nbFilledBytes);
1080
1081      /* This moves the raw bits to take into account the new compressed size */
1082      ec_byte_shrink(&buf, nbCompressedBytes);
1083    }
1084
1085    if (C==2)
1086    {
1087       dual_stereo = stereo_analysis(st->mode, X, st->mode->nbEBands, LM, C, N);
1088       ec_enc_bit_prob(enc, dual_stereo, 32768);
1089    }
1090    if (C==2)
1091    {
1092       int effectiveRate;
1093
1094       /* Account for coarse energy */
1095       effectiveRate = (8*effectiveBytes - 80)>>LM;
1096
1097       /* effectiveRate in kb/s */
1098       effectiveRate = 2*effectiveRate/5;
1099       if (effectiveRate<35)
1100          intensity = 8;
1101       else if (effectiveRate<50)
1102          intensity = 12;
1103       else if (effectiveRate<68)
1104          intensity = 16;
1105       else if (effectiveRate<84)
1106          intensity = 18;
1107       else if (effectiveRate<102)
1108          intensity = 19;
1109       else if (effectiveRate<130)
1110          intensity = 20;
1111       else
1112          intensity = 100;
1113       intensity = IMIN(st->end,IMAX(st->start, intensity));
1114       ec_enc_uint(enc, intensity, 1+st->end-st->start);
1115    }
1116
1117    /* Bit allocation */
1118    ALLOC(fine_quant, st->mode->nbEBands, int);
1119    ALLOC(pulses, st->mode->nbEBands, int);
1120    ALLOC(fine_priority, st->mode->nbEBands, int);
1121
1122    /* bits =   packet size        -       where we are           - safety - skip signalling*/
1123    bits = (nbCompressedBytes*8<<BITRES) - ec_enc_tell(enc, BITRES) - 1    - (1<<BITRES);
1124    codedBands = compute_allocation(st->mode, st->start, st->end, offsets,
1125          alloc_trim, bits, pulses, fine_quant, fine_priority, C, LM, enc, 1, st->lastCodedBands);
1126    st->lastCodedBands = codedBands;
1127
1128    quant_fine_energy(st->mode, st->start, st->end, bandE, oldBandE, error, fine_quant, enc, C);
1129
1130 #ifdef MEASURE_NORM_MSE
1131    float X0[3000];
1132    float bandE0[60];
1133    c=0; do 
1134       for (i=0;i<N;i++)
1135          X0[i+c*N] = X[i+c*N];
1136    while (++c<C);
1137    for (i=0;i<C*st->mode->nbEBands;i++)
1138       bandE0[i] = bandE[i];
1139 #endif
1140
1141    /* Residual quantisation */
1142    quant_all_bands(1, st->mode, st->start, st->end, X, C==2 ? X+N : NULL,
1143          bandE, pulses, shortBlocks, has_fold, dual_stereo, intensity, tf_res, resynth,
1144          nbCompressedBytes*8, enc, LM, codedBands);
1145
1146    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);
1147
1148 #ifdef RESYNTH
1149    /* Re-synthesis of the coded audio if required */
1150    if (resynth)
1151    {
1152       celt_sig *out_mem[2];
1153       celt_sig *overlap_mem[2];
1154
1155       log2Amp(st->mode, st->start, st->end, bandE, oldBandE, C);
1156
1157 #ifdef MEASURE_NORM_MSE
1158       measure_norm_mse(st->mode, X, X0, bandE, bandE0, M, N, C);
1159 #endif
1160
1161       /* Synthesis */
1162       denormalise_bands(st->mode, X, freq, bandE, effEnd, C, M);
1163
1164       CELT_MOVE(st->syn_mem[0], st->syn_mem[0]+N, MAX_PERIOD);
1165       if (C==2)
1166          CELT_MOVE(st->syn_mem[1], st->syn_mem[1]+N, MAX_PERIOD);
1167
1168       c=0; do
1169          for (i=0;i<M*st->mode->eBands[st->start];i++)
1170             freq[c*N+i] = 0;
1171       while (++c<C);
1172       c=0; do
1173          for (i=M*st->mode->eBands[st->end];i<N;i++)
1174             freq[c*N+i] = 0;
1175       while (++c<C);
1176
1177       out_mem[0] = st->syn_mem[0]+MAX_PERIOD;
1178       if (C==2)
1179          out_mem[1] = st->syn_mem[1]+MAX_PERIOD;
1180
1181       c=0; do
1182          overlap_mem[c] = _overlap_mem + c*st->overlap;
1183       while (++c<C);
1184
1185       compute_inv_mdcts(st->mode, shortBlocks, freq, out_mem, overlap_mem, C, LM);
1186
1187 #ifdef ENABLE_POSTFILTER
1188       c=0; do {
1189          st->prefilter_period=IMAX(st->prefilter_period, COMBFILTER_MINPERIOD);
1190          st->prefilter_period_old=IMAX(st->prefilter_period_old, COMBFILTER_MINPERIOD);
1191          if (LM!=0)
1192          {
1193             comb_filter(out_mem[c], out_mem[c], st->prefilter_period, st->prefilter_period, st->overlap, C,
1194                   st->prefilter_gain, st->prefilter_gain, NULL, 0);
1195             comb_filter(out_mem[c]+st->overlap, out_mem[c]+st->overlap, st->prefilter_period, pitch_index, N-st->overlap, C,
1196                   st->prefilter_gain, gain1, st->mode->window, st->mode->overlap);
1197          } else {
1198             comb_filter(out_mem[c], out_mem[c], st->prefilter_period_old, st->prefilter_period, N, C,
1199                   st->prefilter_gain_old, st->prefilter_gain, st->mode->window, st->mode->overlap);
1200          }
1201       } while (++c<C);
1202 #endif /* ENABLE_POSTFILTER */
1203
1204       deemphasis(out_mem, (celt_word16*)pcm, N, C, st->mode->preemph, st->preemph_memD);
1205       st->prefilter_period_old = st->prefilter_period;
1206       st->prefilter_gain_old = st->prefilter_gain;
1207    }
1208 #endif
1209
1210    st->prefilter_period = pitch_index;
1211    st->prefilter_gain = gain1;
1212
1213    /* If there's any room left (can only happen for very high rates),
1214       fill it with zeros */
1215    while (ec_enc_tell(enc,0) + 8 <= nbCompressedBytes*8)
1216       ec_enc_bits(enc, 0, 8);
1217    ec_enc_done(enc);
1218    
1219    RESTORE_STACK;
1220    if (ec_enc_get_error(enc))
1221       return CELT_CORRUPTED_DATA;
1222    else
1223       return nbCompressedBytes;
1224 }
1225
1226 #ifdef FIXED_POINT
1227 #ifndef DISABLE_FLOAT_API
1228 int celt_encode_with_ec_float(CELTEncoder * restrict st, const float * pcm, int frame_size, unsigned char *compressed, int nbCompressedBytes, ec_enc *enc)
1229 {
1230    int j, ret, C, N, LM, M;
1231    VARDECL(celt_int16, in);
1232    SAVE_STACK;
1233
1234    if (pcm==NULL)
1235       return CELT_BAD_ARG;
1236
1237    for (LM=0;LM<4;LM++)
1238       if (st->mode->shortMdctSize<<LM==frame_size)
1239          break;
1240    if (LM>=MAX_CONFIG_SIZES)
1241       return CELT_BAD_ARG;
1242    M=1<<LM;
1243
1244    C = CHANNELS(st->channels);
1245    N = M*st->mode->shortMdctSize;
1246    ALLOC(in, C*N, celt_int16);
1247
1248    for (j=0;j<C*N;j++)
1249      in[j] = FLOAT2INT16(pcm[j]);
1250
1251    ret=celt_encode_with_ec(st,in,frame_size,compressed,nbCompressedBytes, enc);
1252 #ifdef RESYNTH
1253    for (j=0;j<C*N;j++)
1254       ((float*)pcm)[j]=in[j]*(1.f/32768.f);
1255 #endif
1256    RESTORE_STACK;
1257    return ret;
1258
1259 }
1260 #endif /*DISABLE_FLOAT_API*/
1261 #else
1262 int celt_encode_with_ec(CELTEncoder * restrict st, const celt_int16 * pcm, int frame_size, unsigned char *compressed, int nbCompressedBytes, ec_enc *enc)
1263 {
1264    int j, ret, C, N, LM, M;
1265    VARDECL(celt_sig, in);
1266    SAVE_STACK;
1267
1268    if (pcm==NULL)
1269       return CELT_BAD_ARG;
1270
1271    for (LM=0;LM<4;LM++)
1272       if (st->mode->shortMdctSize<<LM==frame_size)
1273          break;
1274    if (LM>=MAX_CONFIG_SIZES)
1275       return CELT_BAD_ARG;
1276    M=1<<LM;
1277
1278    C=CHANNELS(st->channels);
1279    N=M*st->mode->shortMdctSize;
1280    ALLOC(in, C*N, celt_sig);
1281    for (j=0;j<C*N;j++) {
1282      in[j] = SCALEOUT(pcm[j]);
1283    }
1284
1285    ret = celt_encode_with_ec_float(st,in,frame_size,compressed,nbCompressedBytes, enc);
1286 #ifdef RESYNTH
1287    for (j=0;j<C*N;j++)
1288       ((celt_int16*)pcm)[j] = FLOAT2INT16(in[j]);
1289 #endif
1290    RESTORE_STACK;
1291    return ret;
1292 }
1293 #endif
1294
1295 int celt_encode(CELTEncoder * restrict st, const celt_int16 * pcm, int frame_size, unsigned char *compressed, int nbCompressedBytes)
1296 {
1297    return celt_encode_with_ec(st, pcm, frame_size, compressed, nbCompressedBytes, NULL);
1298 }
1299
1300 #ifndef DISABLE_FLOAT_API
1301 int celt_encode_float(CELTEncoder * restrict st, const float * pcm, int frame_size, unsigned char *compressed, int nbCompressedBytes)
1302 {
1303    return celt_encode_with_ec_float(st, pcm, frame_size, compressed, nbCompressedBytes, NULL);
1304 }
1305 #endif /* DISABLE_FLOAT_API */
1306
1307 int celt_encoder_ctl(CELTEncoder * restrict st, int request, ...)
1308 {
1309    va_list ap;
1310    
1311    va_start(ap, request);
1312    switch (request)
1313    {
1314       case CELT_GET_MODE_REQUEST:
1315       {
1316          const CELTMode ** value = va_arg(ap, const CELTMode**);
1317          if (value==0)
1318             goto bad_arg;
1319          *value=st->mode;
1320       }
1321       break;
1322       case CELT_SET_COMPLEXITY_REQUEST:
1323       {
1324          int value = va_arg(ap, celt_int32);
1325          if (value<0 || value>10)
1326             goto bad_arg;
1327          st->complexity = value;
1328       }
1329       break;
1330       case CELT_SET_START_BAND_REQUEST:
1331       {
1332          celt_int32 value = va_arg(ap, celt_int32);
1333          if (value<0 || value>=st->mode->nbEBands)
1334             goto bad_arg;
1335          st->start = value;
1336       }
1337       break;
1338       case CELT_SET_END_BAND_REQUEST:
1339       {
1340          celt_int32 value = va_arg(ap, celt_int32);
1341          if (value<0 || value>=st->mode->nbEBands)
1342             goto bad_arg;
1343          st->end = value;
1344       }
1345       break;
1346       case CELT_SET_PREDICTION_REQUEST:
1347       {
1348          int value = va_arg(ap, celt_int32);
1349          if (value<0 || value>2)
1350             goto bad_arg;
1351          if (value==0)
1352          {
1353             st->force_intra   = 1;
1354          } else if (value==1) {
1355             st->force_intra   = 0;
1356          } else {
1357             st->force_intra   = 0;
1358          }   
1359       }
1360       break;
1361       case CELT_SET_VBR_CONSTRAINT_REQUEST:
1362       {
1363          celt_int32 value = va_arg(ap, celt_int32);
1364          st->constrained_vbr = value;
1365       }
1366       break;
1367       case CELT_SET_VBR_RATE_REQUEST:
1368       {
1369          celt_int32 value = va_arg(ap, celt_int32);
1370          int frame_rate;
1371          int N = st->mode->shortMdctSize;
1372          if (value<0)
1373             goto bad_arg;
1374          if (value>3072000)
1375             value = 3072000;
1376          frame_rate = ((st->mode->Fs<<3)+(N>>1))/N;
1377          st->vbr_rate_norm = ((value<<(BITRES+3))+(frame_rate>>1))/frame_rate;
1378       }
1379       break;
1380       case CELT_RESET_STATE:
1381       {
1382          CELT_MEMSET((char*)&st->ENCODER_RESET_START, 0,
1383                celt_encoder_get_size(st->mode, st->channels)-
1384                ((char*)&st->ENCODER_RESET_START - (char*)st));
1385          st->vbr_offset = 0;
1386          st->delayedIntra = 1;
1387          st->fold_decision = 1;
1388          st->tonal_average = QCONST16(1.f,8);
1389       }
1390       break;
1391       default:
1392          goto bad_request;
1393    }
1394    va_end(ap);
1395    return CELT_OK;
1396 bad_arg:
1397    va_end(ap);
1398    return CELT_BAD_ARG;
1399 bad_request:
1400    va_end(ap);
1401    return CELT_UNIMPLEMENTED;
1402 }
1403
1404 /**********************************************************************/
1405 /*                                                                    */
1406 /*                             DECODER                                */
1407 /*                                                                    */
1408 /**********************************************************************/
1409 #define DECODE_BUFFER_SIZE 2048
1410
1411 /** Decoder state 
1412  @brief Decoder state
1413  */
1414 struct CELTDecoder {
1415    const CELTMode *mode;
1416    int overlap;
1417    int channels;
1418
1419    int start, end;
1420
1421    /* Everything beyond this point gets cleared on a reset */
1422 #define DECODER_RESET_START last_pitch_index
1423
1424    int last_pitch_index;
1425    int loss_count;
1426    int postfilter_period;
1427    int postfilter_period_old;
1428    celt_word16 postfilter_gain;
1429    celt_word16 postfilter_gain_old;
1430
1431    celt_sig preemph_memD[2];
1432    
1433    celt_sig _decode_mem[1]; /* Size = channels*(DECODE_BUFFER_SIZE+mode->overlap) */
1434    /* celt_word16 lpc[],  Size = channels*LPC_ORDER */
1435    /* celt_word16 oldEBands[], Size = channels*mode->nbEBands */
1436 };
1437
1438 int celt_decoder_get_size(const CELTMode *mode, int channels)
1439 {
1440    int size = sizeof(struct CELTDecoder)
1441             + (channels*(DECODE_BUFFER_SIZE+mode->overlap)-1)*sizeof(celt_sig)
1442             + channels*LPC_ORDER*sizeof(celt_word16)
1443             + channels*mode->nbEBands*sizeof(celt_word16);
1444    return size;
1445 }
1446
1447 CELTDecoder *celt_decoder_create(const CELTMode *mode, int channels, int *error)
1448 {
1449    return celt_decoder_init(
1450          (CELTDecoder *)celt_alloc(celt_decoder_get_size(mode, channels)),
1451          mode, channels, error);
1452 }
1453
1454 CELTDecoder *celt_decoder_init(CELTDecoder *st, const CELTMode *mode, int channels, int *error)
1455 {
1456    if (channels < 0 || channels > 2)
1457    {
1458       if (error)
1459          *error = CELT_BAD_ARG;
1460       return NULL;
1461    }
1462
1463    if (st==NULL)
1464    {
1465       if (error)
1466          *error = CELT_ALLOC_FAIL;
1467       return NULL;
1468    }
1469
1470    CELT_MEMSET((char*)st, 0, celt_decoder_get_size(mode, channels));
1471
1472    st->mode = mode;
1473    st->overlap = mode->overlap;
1474    st->channels = channels;
1475
1476    st->start = 0;
1477    st->end = st->mode->effEBands;
1478
1479    st->loss_count = 0;
1480
1481    if (error)
1482       *error = CELT_OK;
1483    return st;
1484 }
1485
1486 void celt_decoder_destroy(CELTDecoder *st)
1487 {
1488    celt_free(st);
1489 }
1490
1491 static void celt_decode_lost(CELTDecoder * restrict st, celt_word16 * restrict pcm, int N, int LM)
1492 {
1493    int c;
1494    int pitch_index;
1495    int overlap = st->mode->overlap;
1496    celt_word16 fade = Q15ONE;
1497    int i, len;
1498    const int C = CHANNELS(st->channels);
1499    int offset;
1500    celt_sig *out_mem[2];
1501    celt_sig *decode_mem[2];
1502    celt_sig *overlap_mem[2];
1503    celt_word16 *lpc;
1504    SAVE_STACK;
1505    
1506    c=0; do {
1507       decode_mem[c] = st->_decode_mem + c*(DECODE_BUFFER_SIZE+st->overlap);
1508       out_mem[c] = decode_mem[c]+DECODE_BUFFER_SIZE-MAX_PERIOD;
1509       overlap_mem[c] = decode_mem[c]+DECODE_BUFFER_SIZE;
1510    } while (++c<C);
1511    lpc = (celt_word16*)(st->_decode_mem+(DECODE_BUFFER_SIZE+st->overlap)*C);
1512
1513    len = N+st->mode->overlap;
1514    
1515    if (st->loss_count == 0)
1516    {
1517       celt_word16 pitch_buf[MAX_PERIOD>>1];
1518       celt_word32 tmp=0;
1519       celt_word32 mem0[2]={0,0};
1520       celt_word16 mem1[2]={0,0};
1521       int len2 = len;
1522       /* FIXME: This is a kludge */
1523       if (len2>MAX_PERIOD>>1)
1524          len2 = MAX_PERIOD>>1;
1525       pitch_downsample(out_mem, pitch_buf, MAX_PERIOD, MAX_PERIOD,
1526                        C, mem0, mem1);
1527       pitch_search(st->mode, pitch_buf+((MAX_PERIOD-len2)>>1), pitch_buf, len2,
1528                    MAX_PERIOD-len2-100, &pitch_index, &tmp, 1<<LM);
1529       pitch_index = MAX_PERIOD-len2-pitch_index;
1530       st->last_pitch_index = pitch_index;
1531    } else {
1532       pitch_index = st->last_pitch_index;
1533       if (st->loss_count < 5)
1534          fade = QCONST16(.8f,15);
1535       else
1536          fade = 0;
1537    }
1538
1539    c=0; do {
1540       /* FIXME: This is more memory than necessary */
1541       celt_word32 e[2*MAX_PERIOD];
1542       celt_word16 exc[2*MAX_PERIOD];
1543       celt_word32 ac[LPC_ORDER+1];
1544       celt_word16 decay = 1;
1545       celt_word32 S1=0;
1546       celt_word16 mem[LPC_ORDER]={0};
1547
1548       offset = MAX_PERIOD-pitch_index;
1549       for (i=0;i<MAX_PERIOD;i++)
1550          exc[i] = ROUND16(out_mem[c][i], SIG_SHIFT);
1551
1552       if (st->loss_count == 0)
1553       {
1554          _celt_autocorr(exc, ac, st->mode->window, st->mode->overlap,
1555                         LPC_ORDER, MAX_PERIOD);
1556
1557          /* Noise floor -40 dB */
1558 #ifdef FIXED_POINT
1559          ac[0] += SHR32(ac[0],13);
1560 #else
1561          ac[0] *= 1.0001f;
1562 #endif
1563          /* Lag windowing */
1564          for (i=1;i<=LPC_ORDER;i++)
1565          {
1566             /*ac[i] *= exp(-.5*(2*M_PI*.002*i)*(2*M_PI*.002*i));*/
1567 #ifdef FIXED_POINT
1568             ac[i] -= MULT16_32_Q15(2*i*i, ac[i]);
1569 #else
1570             ac[i] -= ac[i]*(.008f*i)*(.008f*i);
1571 #endif
1572          }
1573
1574          _celt_lpc(lpc+c*LPC_ORDER, ac, LPC_ORDER);
1575       }
1576       for (i=0;i<LPC_ORDER;i++)
1577          mem[i] = ROUND16(out_mem[c][MAX_PERIOD-1-i], SIG_SHIFT);
1578       fir(exc, lpc+c*LPC_ORDER, exc, MAX_PERIOD, LPC_ORDER, mem);
1579       /*for (i=0;i<MAX_PERIOD;i++)printf("%d ", exc[i]); printf("\n");*/
1580       /* Check if the waveform is decaying (and if so how fast) */
1581       {
1582          celt_word32 E1=1, E2=1;
1583          int period;
1584          if (pitch_index <= MAX_PERIOD/2)
1585             period = pitch_index;
1586          else
1587             period = MAX_PERIOD/2;
1588          for (i=0;i<period;i++)
1589          {
1590             E1 += SHR32(MULT16_16(exc[MAX_PERIOD-period+i],exc[MAX_PERIOD-period+i]),8);
1591             E2 += SHR32(MULT16_16(exc[MAX_PERIOD-2*period+i],exc[MAX_PERIOD-2*period+i]),8);
1592          }
1593          if (E1 > E2)
1594             E1 = E2;
1595          decay = celt_sqrt(frac_div32(SHR(E1,1),E2));
1596       }
1597
1598       /* Copy excitation, taking decay into account */
1599       for (i=0;i<len+st->mode->overlap;i++)
1600       {
1601          celt_word16 tmp;
1602          if (offset+i >= MAX_PERIOD)
1603          {
1604             offset -= pitch_index;
1605             decay = MULT16_16_Q15(decay, decay);
1606          }
1607          e[i] = SHL32(EXTEND32(MULT16_16_Q15(decay, exc[offset+i])), SIG_SHIFT);
1608          tmp = ROUND16(out_mem[c][offset+i],SIG_SHIFT);
1609          S1 += SHR32(MULT16_16(tmp,tmp),8);
1610       }
1611       for (i=0;i<LPC_ORDER;i++)
1612          mem[i] = ROUND16(out_mem[c][MAX_PERIOD-1-i], SIG_SHIFT);
1613       for (i=0;i<len+st->mode->overlap;i++)
1614          e[i] = MULT16_32_Q15(fade, e[i]);
1615       iir(e, lpc+c*LPC_ORDER, e, len+st->mode->overlap, LPC_ORDER, mem);
1616
1617       {
1618          celt_word32 S2=0;
1619          for (i=0;i<len+overlap;i++)
1620          {
1621             celt_word16 tmp = ROUND16(e[i],SIG_SHIFT);
1622             S2 += SHR32(MULT16_16(tmp,tmp),8);
1623          }
1624          /* This checks for an "explosion" in the synthesis */
1625 #ifdef FIXED_POINT
1626          if (!(S1 > SHR32(S2,2)))
1627 #else
1628          /* Float test is written this way to catch NaNs at the same time */
1629          if (!(S1 > 0.2f*S2))
1630 #endif
1631          {
1632             for (i=0;i<len+overlap;i++)
1633                e[i] = 0;
1634          } else if (S1 < S2)
1635          {
1636             celt_word16 ratio = celt_sqrt(frac_div32(SHR32(S1,1)+1,S2+1));
1637             for (i=0;i<len+overlap;i++)
1638                e[i] = MULT16_32_Q15(ratio, e[i]);
1639          }
1640       }
1641
1642 #ifdef ENABLE_POSTFILTER
1643       /* Apply post-filter to the MDCT overlap of the previous frame */
1644       comb_filter(out_mem[c]+MAX_PERIOD, out_mem[c]+MAX_PERIOD, st->postfilter_period, st->postfilter_period, st->overlap, C,
1645                   st->postfilter_gain, st->postfilter_gain, NULL, 0);
1646 #endif /* ENABLE_POSTFILTER */
1647
1648       for (i=0;i<MAX_PERIOD+st->mode->overlap-N;i++)
1649          out_mem[c][i] = out_mem[c][N+i];
1650
1651       /* Apply TDAC to the concealed audio so that it blends with the
1652          previous and next frames */
1653       for (i=0;i<overlap/2;i++)
1654       {
1655          celt_word32 tmp;
1656          tmp = MULT16_32_Q15(st->mode->window[i],           e[N+overlap-1-i]) +
1657                MULT16_32_Q15(st->mode->window[overlap-i-1], e[N+i          ]);
1658          out_mem[c][MAX_PERIOD+i] = MULT16_32_Q15(st->mode->window[overlap-i-1], tmp);
1659          out_mem[c][MAX_PERIOD+overlap-i-1] = MULT16_32_Q15(st->mode->window[i], tmp);
1660       }
1661       for (i=0;i<N;i++)
1662          out_mem[c][MAX_PERIOD-N+i] = e[i];
1663
1664 #ifdef ENABLE_POSTFILTER
1665       /* Apply pre-filter to the MDCT overlap for the next frame (post-filter will be applied then) */
1666       comb_filter(e, out_mem[c]+MAX_PERIOD, st->postfilter_period, st->postfilter_period, st->overlap, C,
1667                   -st->postfilter_gain, -st->postfilter_gain, NULL, 0);
1668 #endif /* ENABLE_POSTFILTER */
1669       for (i=0;i<overlap;i++)
1670          out_mem[c][MAX_PERIOD+i] = e[i];
1671    } while (++c<C);
1672
1673    {
1674       celt_word32 *out_syn[2];
1675       out_syn[0] = out_mem[0]+MAX_PERIOD-N;
1676       if (C==2)
1677          out_syn[1] = out_mem[1]+MAX_PERIOD-N;
1678       deemphasis(out_syn, pcm, N, C, st->mode->preemph, st->preemph_memD);
1679    }
1680    
1681    st->loss_count++;
1682
1683    RESTORE_STACK;
1684 }
1685
1686 #ifdef FIXED_POINT
1687 int celt_decode_with_ec(CELTDecoder * restrict st, const unsigned char *data, int len, celt_int16 * restrict pcm, int frame_size, ec_dec *dec)
1688 {
1689 #else
1690 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)
1691 {
1692 #endif
1693    int c, i, N;
1694    int has_fold;
1695    int bits;
1696    ec_dec _dec;
1697    ec_byte_buffer buf;
1698    VARDECL(celt_sig, freq);
1699    VARDECL(celt_norm, X);
1700    VARDECL(celt_ener, bandE);
1701    VARDECL(int, fine_quant);
1702    VARDECL(int, pulses);
1703    VARDECL(int, offsets);
1704    VARDECL(int, fine_priority);
1705    VARDECL(int, tf_res);
1706    celt_sig *out_mem[2];
1707    celt_sig *decode_mem[2];
1708    celt_sig *overlap_mem[2];
1709    celt_sig *out_syn[2];
1710    celt_word16 *lpc;
1711    celt_word16 *oldBandE;
1712
1713    int shortBlocks;
1714    int isTransient;
1715    int intra_ener;
1716    const int C = CHANNELS(st->channels);
1717    int LM, M;
1718    int nbFilledBytes, nbAvailableBytes;
1719    int effEnd;
1720    int codedBands;
1721    int alloc_trim;
1722    int postfilter_pitch;
1723    celt_word16 postfilter_gain;
1724    int intensity=0;
1725    int dual_stereo=0;
1726    SAVE_STACK;
1727
1728    if (pcm==NULL)
1729       return CELT_BAD_ARG;
1730
1731    for (LM=0;LM<4;LM++)
1732       if (st->mode->shortMdctSize<<LM==frame_size)
1733          break;
1734    if (LM>=MAX_CONFIG_SIZES)
1735       return CELT_BAD_ARG;
1736    M=1<<LM;
1737
1738    c=0; do {
1739       decode_mem[c] = st->_decode_mem + c*(DECODE_BUFFER_SIZE+st->overlap);
1740       out_mem[c] = decode_mem[c]+DECODE_BUFFER_SIZE-MAX_PERIOD;
1741       overlap_mem[c] = decode_mem[c]+DECODE_BUFFER_SIZE;
1742    } while (++c<C);
1743    lpc = (celt_word16*)(st->_decode_mem+(DECODE_BUFFER_SIZE+st->overlap)*C);
1744    oldBandE = lpc+C*LPC_ORDER;
1745
1746    N = M*st->mode->shortMdctSize;
1747
1748    effEnd = st->end;
1749    if (effEnd > st->mode->effEBands)
1750       effEnd = st->mode->effEBands;
1751
1752    ALLOC(freq, C*N, celt_sig); /**< Interleaved signal MDCTs */
1753    ALLOC(X, C*N, celt_norm);   /**< Interleaved normalised MDCTs */
1754    ALLOC(bandE, st->mode->nbEBands*C, celt_ener);
1755    c=0; do
1756       for (i=0;i<M*st->mode->eBands[st->start];i++)
1757          X[c*N+i] = 0;
1758    while (++c<C);
1759    c=0; do   
1760       for (i=M*st->mode->eBands[effEnd];i<N;i++)
1761          X[c*N+i] = 0;
1762    while (++c<C);
1763
1764    if (data == NULL)
1765    {
1766       celt_decode_lost(st, pcm, N, LM);
1767       RESTORE_STACK;
1768       return CELT_OK;
1769    }
1770    if (len<0) {
1771      RESTORE_STACK;
1772      return CELT_BAD_ARG;
1773    }
1774    
1775    if (dec == NULL)
1776    {
1777       ec_byte_readinit(&buf,(unsigned char*)data,len);
1778       ec_dec_init(&_dec,&buf);
1779       dec = &_dec;
1780       nbFilledBytes = 0;
1781    } else {
1782       nbFilledBytes = (ec_dec_tell(dec, 0)+4)>>3;
1783    }
1784    nbAvailableBytes = len-nbFilledBytes;
1785
1786    if (ec_dec_bit_prob(dec, 32768))
1787    {
1788 #ifdef ENABLE_POSTFILTER
1789       int qg, octave;
1790       octave = ec_dec_uint(dec, 6);
1791       postfilter_pitch = (16<<octave)+ec_dec_bits(dec, 4+octave);
1792       qg = ec_dec_bits(dec, 2);
1793       postfilter_gain = QCONST16(.125f,15)*(qg+2);
1794 #else /* ENABLE_POSTFILTER */
1795       RESTORE_STACK;
1796       return CELT_CORRUPTED_DATA;
1797 #endif /* ENABLE_POSTFILTER */
1798
1799    } else {
1800       postfilter_gain = 0;
1801       postfilter_pitch = 0;
1802    }
1803
1804    /* Decode the global flags (first symbols in the stream) */
1805    intra_ener = ec_dec_bit_prob(dec, 8192);
1806    /* Get band energies */
1807    unquant_coarse_energy(st->mode, st->start, st->end, bandE, oldBandE,
1808          intra_ener, dec, C, LM);
1809
1810    if (LM > 0)
1811       isTransient = ec_dec_bit_prob(dec, 8192);
1812    else
1813       isTransient = 0;
1814
1815    if (isTransient)
1816       shortBlocks = M;
1817    else
1818       shortBlocks = 0;
1819
1820    ALLOC(tf_res, st->mode->nbEBands, int);
1821    tf_decode(st->start, st->end, C, isTransient, tf_res, LM, dec);
1822
1823    has_fold = ec_dec_bit_prob(dec, 8192)<<1;
1824    has_fold |= ec_dec_bit_prob(dec, (has_fold>>1) ? 32768 : 49152);
1825
1826    ALLOC(pulses, st->mode->nbEBands, int);
1827    ALLOC(offsets, st->mode->nbEBands, int);
1828    ALLOC(fine_priority, st->mode->nbEBands, int);
1829
1830    for (i=0;i<st->mode->nbEBands;i++)
1831       offsets[i] = 0;
1832    for (i=0;i<st->mode->nbEBands;i++)
1833    {
1834       if (ec_dec_bit_prob(dec, 1024))
1835       {
1836          while (ec_dec_bit_prob(dec, 32768))
1837             offsets[i]++;
1838          offsets[i]++;
1839          offsets[i] *= (6<<BITRES);
1840       }
1841    }
1842
1843    ALLOC(fine_quant, st->mode->nbEBands, int);
1844    {
1845       int fl;
1846       alloc_trim = 0;
1847       fl = ec_decode_bin(dec, 7);
1848       while (trim_cdf[alloc_trim+1] <= fl)
1849          alloc_trim++;
1850       ec_dec_update(dec, trim_cdf[alloc_trim], trim_cdf[alloc_trim+1], 128);
1851    }
1852
1853    if (C==2)
1854    {
1855       dual_stereo = ec_dec_bit_prob(dec, 32768);
1856       intensity = ec_dec_uint(dec, 1+st->end-st->start);
1857    }
1858
1859    bits = (len*8<<BITRES) - ec_dec_tell(dec, BITRES) - 1 - (1<<BITRES);
1860    codedBands = compute_allocation(st->mode, st->start, st->end, offsets,
1861          alloc_trim, bits, pulses, fine_quant, fine_priority, C, LM, dec, 0, 0);
1862    
1863    unquant_fine_energy(st->mode, st->start, st->end, bandE, oldBandE, fine_quant, dec, C);
1864
1865    /* Decode fixed codebook */
1866    quant_all_bands(0, st->mode, st->start, st->end, X, C==2 ? X+N : NULL,
1867          NULL, pulses, shortBlocks, has_fold, dual_stereo, intensity, tf_res, 1,
1868          len*8, dec, LM, codedBands);
1869
1870    unquant_energy_finalise(st->mode, st->start, st->end, bandE, oldBandE,
1871          fine_quant, fine_priority, len*8-ec_dec_tell(dec, 0), dec, C);
1872
1873    log2Amp(st->mode, st->start, st->end, bandE, oldBandE, C);
1874
1875    /* Synthesis */
1876    denormalise_bands(st->mode, X, freq, bandE, effEnd, C, M);
1877
1878    CELT_MOVE(decode_mem[0], decode_mem[0]+N, DECODE_BUFFER_SIZE-N);
1879    if (C==2)
1880       CELT_MOVE(decode_mem[1], decode_mem[1]+N, DECODE_BUFFER_SIZE-N);
1881
1882    c=0; do
1883       for (i=0;i<M*st->mode->eBands[st->start];i++)
1884          freq[c*N+i] = 0;
1885    while (++c<C);
1886    c=0; do
1887       for (i=M*st->mode->eBands[effEnd];i<N;i++)
1888          freq[c*N+i] = 0;
1889    while (++c<C);
1890
1891    out_syn[0] = out_mem[0]+MAX_PERIOD-N;
1892    if (C==2)
1893       out_syn[1] = out_mem[1]+MAX_PERIOD-N;
1894
1895    /* Compute inverse MDCTs */
1896    compute_inv_mdcts(st->mode, shortBlocks, freq, out_syn, overlap_mem, C, LM);
1897
1898 #ifdef ENABLE_POSTFILTER
1899    c=0; do {
1900       st->postfilter_period=IMAX(st->postfilter_period, COMBFILTER_MINPERIOD);
1901       st->postfilter_period_old=IMAX(st->postfilter_period_old, COMBFILTER_MINPERIOD);
1902       if (LM!=0)
1903       {
1904          comb_filter(out_syn[c], out_syn[c], st->postfilter_period, st->postfilter_period, st->overlap, C,
1905                st->postfilter_gain, st->postfilter_gain, NULL, 0);
1906          comb_filter(out_syn[c]+st->overlap, out_syn[c]+st->overlap, st->postfilter_period, postfilter_pitch, N-st->overlap, C,
1907                st->postfilter_gain, postfilter_gain, st->mode->window, st->mode->overlap);
1908       } else {
1909          comb_filter(out_syn[c], out_syn[c], st->postfilter_period_old, st->postfilter_period, N-st->overlap, C,
1910                st->postfilter_gain_old, st->postfilter_gain, st->mode->window, st->mode->overlap);
1911       }
1912    } while (++c<C);
1913    st->postfilter_period_old = st->postfilter_period;
1914    st->postfilter_gain_old = st->postfilter_gain;
1915    st->postfilter_period = postfilter_pitch;
1916    st->postfilter_gain = postfilter_gain;
1917 #endif /* ENABLE_POSTFILTER */
1918
1919    deemphasis(out_syn, pcm, N, C, st->mode->preemph, st->preemph_memD);
1920    st->loss_count = 0;
1921    RESTORE_STACK;
1922    if (ec_dec_get_error(dec))
1923       return CELT_CORRUPTED_DATA;
1924    else
1925       return CELT_OK;
1926 }
1927
1928 #ifdef FIXED_POINT
1929 #ifndef DISABLE_FLOAT_API
1930 int celt_decode_with_ec_float(CELTDecoder * restrict st, const unsigned char *data, int len, float * restrict pcm, int frame_size, ec_dec *dec)
1931 {
1932    int j, ret, C, N, LM, M;
1933    VARDECL(celt_int16, out);
1934    SAVE_STACK;
1935
1936    if (pcm==NULL)
1937       return CELT_BAD_ARG;
1938
1939    for (LM=0;LM<4;LM++)
1940       if (st->mode->shortMdctSize<<LM==frame_size)
1941          break;
1942    if (LM>=MAX_CONFIG_SIZES)
1943       return CELT_BAD_ARG;
1944    M=1<<LM;
1945
1946    C = CHANNELS(st->channels);
1947    N = M*st->mode->shortMdctSize;
1948    
1949    ALLOC(out, C*N, celt_int16);
1950    ret=celt_decode_with_ec(st, data, len, out, frame_size, dec);
1951    if (ret==0)
1952       for (j=0;j<C*N;j++)
1953          pcm[j]=out[j]*(1.f/32768.f);
1954      
1955    RESTORE_STACK;
1956    return ret;
1957 }
1958 #endif /*DISABLE_FLOAT_API*/
1959 #else
1960 int celt_decode_with_ec(CELTDecoder * restrict st, const unsigned char *data, int len, celt_int16 * restrict pcm, int frame_size, ec_dec *dec)
1961 {
1962    int j, ret, C, N, LM, M;
1963    VARDECL(celt_sig, out);
1964    SAVE_STACK;
1965
1966    if (pcm==NULL)
1967       return CELT_BAD_ARG;
1968
1969    for (LM=0;LM<4;LM++)
1970       if (st->mode->shortMdctSize<<LM==frame_size)
1971          break;
1972    if (LM>=MAX_CONFIG_SIZES)
1973       return CELT_BAD_ARG;
1974    M=1<<LM;
1975
1976    C = CHANNELS(st->channels);
1977    N = M*st->mode->shortMdctSize;
1978    ALLOC(out, C*N, celt_sig);
1979
1980    ret=celt_decode_with_ec_float(st, data, len, out, frame_size, dec);
1981
1982    if (ret==0)
1983       for (j=0;j<C*N;j++)
1984          pcm[j] = FLOAT2INT16 (out[j]);
1985    
1986    RESTORE_STACK;
1987    return ret;
1988 }
1989 #endif
1990
1991 int celt_decode(CELTDecoder * restrict st, const unsigned char *data, int len, celt_int16 * restrict pcm, int frame_size)
1992 {
1993    return celt_decode_with_ec(st, data, len, pcm, frame_size, NULL);
1994 }
1995
1996 #ifndef DISABLE_FLOAT_API
1997 int celt_decode_float(CELTDecoder * restrict st, const unsigned char *data, int len, float * restrict pcm, int frame_size)
1998 {
1999    return celt_decode_with_ec_float(st, data, len, pcm, frame_size, NULL);
2000 }
2001 #endif /* DISABLE_FLOAT_API */
2002
2003 int celt_decoder_ctl(CELTDecoder * restrict st, int request, ...)
2004 {
2005    va_list ap;
2006
2007    va_start(ap, request);
2008    switch (request)
2009    {
2010       case CELT_GET_MODE_REQUEST:
2011       {
2012          const CELTMode ** value = va_arg(ap, const CELTMode**);
2013          if (value==0)
2014             goto bad_arg;
2015          *value=st->mode;
2016       }
2017       break;
2018       case CELT_SET_START_BAND_REQUEST:
2019       {
2020          celt_int32 value = va_arg(ap, celt_int32);
2021          if (value<0 || value>=st->mode->nbEBands)
2022             goto bad_arg;
2023          st->start = value;
2024       }
2025       break;
2026       case CELT_SET_END_BAND_REQUEST:
2027       {
2028          celt_int32 value = va_arg(ap, celt_int32);
2029          if (value<0 || value>=st->mode->nbEBands)
2030             goto bad_arg;
2031          st->end = value;
2032       }
2033       break;
2034       case CELT_RESET_STATE:
2035       {
2036          CELT_MEMSET((char*)&st->DECODER_RESET_START, 0,
2037                celt_decoder_get_size(st->mode, st->channels)-
2038                ((char*)&st->DECODER_RESET_START - (char*)st));
2039       }
2040       break;
2041       default:
2042          goto bad_request;
2043    }
2044    va_end(ap);
2045    return CELT_OK;
2046 bad_arg:
2047    va_end(ap);
2048    return CELT_BAD_ARG;
2049 bad_request:
2050       va_end(ap);
2051   return CELT_UNIMPLEMENTED;
2052 }
2053
2054 const char *celt_strerror(int error)
2055 {
2056    static const char *error_strings[8] = {
2057       "success",
2058       "invalid argument",
2059       "invalid mode",
2060       "internal error",
2061       "corrupted stream",
2062       "request not implemented",
2063       "invalid state",
2064       "memory allocation failed"
2065    };
2066    if (error > 0 || error < -7)
2067       return "unknown error";
2068    else 
2069       return error_strings[-error];
2070 }
2071