582d306f01639afabf9e9a63270d6ec48cd96df6
[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 (nbAvailableBytes>12*C)
826       {
827          VARDECL(celt_word16, pitch_buf);
828          ALLOC(pitch_buf, (COMBFILTER_MAXPERIOD+N)>>1, celt_word16);
829          celt_word32 tmp=0;
830          celt_word32 mem0[2]={0,0};
831          celt_word16 mem1[2]={0,0};
832
833          pitch_downsample(pre, pitch_buf, COMBFILTER_MAXPERIOD+N, COMBFILTER_MAXPERIOD+N,
834                           C, mem0, mem1);
835          pitch_search(st->mode, pitch_buf+(COMBFILTER_MAXPERIOD>>1), pitch_buf, N,
836                COMBFILTER_MAXPERIOD-COMBFILTER_MINPERIOD, &pitch_index, &tmp, 1<<LM);
837          pitch_index = COMBFILTER_MAXPERIOD-pitch_index;
838
839          gain1 = remove_doubling(pitch_buf, COMBFILTER_MAXPERIOD, COMBFILTER_MINPERIOD,
840                N, &pitch_index, st->prefilter_period, st->prefilter_gain);
841          if (pitch_index > COMBFILTER_MAXPERIOD)
842             pitch_index = COMBFILTER_MAXPERIOD;
843          gain1 = MULT16_16_Q15(QCONST16(.7f,15),gain1);
844          if (gain1 > QCONST16(.6f,15))
845             gain1 = QCONST16(.6f,15);
846          if (ABS16(gain1-st->prefilter_gain)<QCONST16(.1,15))
847             gain1=st->prefilter_gain;
848       } else {
849          gain1 = 0;
850       }
851       if (gain1<QCONST16(.2f,15) || (nbAvailableBytes<30 && gain1<QCONST16(.4f,15)))
852       {
853          ec_enc_bit_prob(enc, 0, 32768);
854          gain1 = 0;
855       } else {
856          int qg;
857          int octave;
858 #ifdef FIXED_POINT
859          qg = ((gain1+2048)>>12)-2;
860 #else
861          qg = floor(.5+gain1*8)-2;
862 #endif
863          ec_enc_bit_prob(enc, 1, 32768);
864          octave = EC_ILOG(pitch_index)-5;
865          ec_enc_uint(enc, octave, 6);
866          ec_enc_bits(enc, pitch_index-(16<<octave), 4+octave);
867          ec_enc_bits(enc, qg, 2);
868          gain1 = QCONST16(.125f,15)*(qg+2);
869       }
870       /*printf("%d %f\n", pitch_index, gain1);*/
871 #else /* ENABLE_POSTFILTER */
872       ec_enc_bit_prob(enc, 0, 32768);
873 #endif /* ENABLE_POSTFILTER */
874
875       c=0; do {
876          st->prefilter_period=IMAX(st->prefilter_period, COMBFILTER_MINPERIOD);
877          CELT_COPY(in+c*(N+st->overlap), st->in_mem+c*(st->overlap), st->overlap);
878 #ifdef ENABLE_POSTFILTER
879          comb_filter(in+c*(N+st->overlap)+st->overlap, pre[c]+COMBFILTER_MAXPERIOD,
880                st->prefilter_period, pitch_index, N, C, -st->prefilter_gain, -gain1, st->mode->window, st->mode->overlap);
881 #endif /* ENABLE_POSTFILTER */
882          CELT_COPY(st->in_mem+c*(st->overlap), in+c*(N+st->overlap)+N, st->overlap);
883
884 #ifdef ENABLE_POSTFILTER
885          if (N>COMBFILTER_MAXPERIOD)
886          {
887             CELT_MOVE(prefilter_mem+c*COMBFILTER_MAXPERIOD, pre[c]+N, COMBFILTER_MAXPERIOD);
888          } else {
889             CELT_MOVE(prefilter_mem+c*COMBFILTER_MAXPERIOD, prefilter_mem+c*COMBFILTER_MAXPERIOD+N, COMBFILTER_MAXPERIOD-N);
890             CELT_MOVE(prefilter_mem+c*COMBFILTER_MAXPERIOD+COMBFILTER_MAXPERIOD-N, pre[c]+COMBFILTER_MAXPERIOD, N);
891          }
892 #endif /* ENABLE_POSTFILTER */
893       } while (++c<C);
894
895       RESTORE_STACK;
896    }
897
898 #ifdef RESYNTH
899    resynth = 1;
900 #else
901    resynth = 0;
902 #endif
903
904    if (st->complexity > 1 && LM>0)
905    {
906       isTransient = M > 1 &&
907          transient_analysis(in, N+st->overlap, C, &st->frame_max, st->overlap);
908    } else {
909       isTransient = 0;
910    }
911
912    if (isTransient)
913       shortBlocks = M;
914    else
915       shortBlocks = 0;
916
917    ALLOC(freq, C*N, celt_sig); /**< Interleaved signal MDCTs */
918    ALLOC(bandE,st->mode->nbEBands*C, celt_ener);
919    ALLOC(bandLogE,st->mode->nbEBands*C, celt_word16);
920    /* Compute MDCTs */
921    compute_mdcts(st->mode, shortBlocks, in, freq, C, LM);
922
923    ALLOC(X, C*N, celt_norm);         /**< Interleaved normalised MDCTs */
924
925    compute_band_energies(st->mode, freq, bandE, effEnd, C, M);
926
927    amp2Log2(st->mode, effEnd, st->end, bandE, bandLogE, C);
928
929    /* Band normalisation */
930    normalise_bands(st->mode, freq, X, bandE, effEnd, C, M);
931
932    ALLOC(tf_res, st->mode->nbEBands, int);
933    /* Needs to be before coarse energy quantization because otherwise the energy gets modified */
934    tf_select = tf_analysis(st->mode, bandLogE, oldBandE, effEnd, C, isTransient, tf_res, effectiveBytes, X, N, LM, &tf_sum);
935    for (i=effEnd;i<st->end;i++)
936       tf_res[i] = tf_res[effEnd-1];
937
938    ALLOC(error, C*st->mode->nbEBands, celt_word16);
939    quant_coarse_energy(st->mode, st->start, st->end, effEnd, bandLogE,
940          oldBandE, nbCompressedBytes*8, error, enc,
941          C, LM, nbAvailableBytes, st->force_intra,
942          &st->delayedIntra, st->complexity >= 4);
943
944    if (LM > 0)
945       ec_enc_bit_prob(enc, shortBlocks!=0, 8192);
946
947    tf_encode(st->start, st->end, isTransient, tf_res, LM, tf_select, enc);
948
949    if (shortBlocks || st->complexity < 3 || nbAvailableBytes < 10*C)
950    {
951       if (st->complexity == 0)
952       {
953          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       /* Always use MS for 2.5 ms frames until we can do a better analysis */
1088       if (LM==0)
1089          dual_stereo = 0;
1090       else
1091          dual_stereo = stereo_analysis(st->mode, X, st->mode->nbEBands, LM, C, N);
1092       ec_enc_bit_prob(enc, dual_stereo, 32768);
1093    }
1094    if (C==2)
1095    {
1096       int effectiveRate;
1097
1098       /* Account for coarse energy */
1099       effectiveRate = (8*effectiveBytes - 80)>>LM;
1100
1101       /* effectiveRate in kb/s */
1102       effectiveRate = 2*effectiveRate/5;
1103       if (effectiveRate<35)
1104          intensity = 8;
1105       else if (effectiveRate<50)
1106          intensity = 12;
1107       else if (effectiveRate<68)
1108          intensity = 16;
1109       else if (effectiveRate<84)
1110          intensity = 18;
1111       else if (effectiveRate<102)
1112          intensity = 19;
1113       else if (effectiveRate<130)
1114          intensity = 20;
1115       else
1116          intensity = 100;
1117       intensity = IMIN(st->end,IMAX(st->start, intensity));
1118       ec_enc_uint(enc, intensity, 1+st->end-st->start);
1119    }
1120
1121    /* Bit allocation */
1122    ALLOC(fine_quant, st->mode->nbEBands, int);
1123    ALLOC(pulses, st->mode->nbEBands, int);
1124    ALLOC(fine_priority, st->mode->nbEBands, int);
1125
1126    /* bits =   packet size        -       where we are           - safety */
1127    bits = (nbCompressedBytes*8<<BITRES) - ec_enc_tell(enc, BITRES) - 1;
1128    codedBands = compute_allocation(st->mode, st->start, st->end, offsets,
1129          alloc_trim, bits, pulses, fine_quant, fine_priority, C, LM, enc, 1, st->lastCodedBands);
1130    st->lastCodedBands = codedBands;
1131
1132    quant_fine_energy(st->mode, st->start, st->end, bandE, oldBandE, error, fine_quant, enc, C);
1133
1134 #ifdef MEASURE_NORM_MSE
1135    float X0[3000];
1136    float bandE0[60];
1137    c=0; do 
1138       for (i=0;i<N;i++)
1139          X0[i+c*N] = X[i+c*N];
1140    while (++c<C);
1141    for (i=0;i<C*st->mode->nbEBands;i++)
1142       bandE0[i] = bandE[i];
1143 #endif
1144
1145    /* Residual quantisation */
1146    quant_all_bands(1, st->mode, st->start, st->end, X, C==2 ? X+N : NULL,
1147          bandE, pulses, shortBlocks, has_fold, dual_stereo, intensity, tf_res, resynth,
1148          nbCompressedBytes*8, enc, LM, codedBands);
1149
1150    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);
1151
1152 #ifdef RESYNTH
1153    /* Re-synthesis of the coded audio if required */
1154    if (resynth)
1155    {
1156       celt_sig *out_mem[2];
1157       celt_sig *overlap_mem[2];
1158
1159       log2Amp(st->mode, st->start, st->end, bandE, oldBandE, C);
1160
1161 #ifdef MEASURE_NORM_MSE
1162       measure_norm_mse(st->mode, X, X0, bandE, bandE0, M, N, C);
1163 #endif
1164
1165       /* Synthesis */
1166       denormalise_bands(st->mode, X, freq, bandE, effEnd, C, M);
1167
1168       CELT_MOVE(st->syn_mem[0], st->syn_mem[0]+N, MAX_PERIOD);
1169       if (C==2)
1170          CELT_MOVE(st->syn_mem[1], st->syn_mem[1]+N, MAX_PERIOD);
1171
1172       c=0; do
1173          for (i=0;i<M*st->mode->eBands[st->start];i++)
1174             freq[c*N+i] = 0;
1175       while (++c<C);
1176       c=0; do
1177          for (i=M*st->mode->eBands[st->end];i<N;i++)
1178             freq[c*N+i] = 0;
1179       while (++c<C);
1180
1181       out_mem[0] = st->syn_mem[0]+MAX_PERIOD;
1182       if (C==2)
1183          out_mem[1] = st->syn_mem[1]+MAX_PERIOD;
1184
1185       c=0; do
1186          overlap_mem[c] = _overlap_mem + c*st->overlap;
1187       while (++c<C);
1188
1189       compute_inv_mdcts(st->mode, shortBlocks, freq, out_mem, overlap_mem, C, LM);
1190
1191 #ifdef ENABLE_POSTFILTER
1192       c=0; do {
1193          st->prefilter_period=IMAX(st->prefilter_period, COMBFILTER_MINPERIOD);
1194          st->prefilter_period_old=IMAX(st->prefilter_period_old, COMBFILTER_MINPERIOD);
1195          if (LM!=0)
1196          {
1197             comb_filter(out_mem[c], out_mem[c], st->prefilter_period, st->prefilter_period, st->overlap, C,
1198                   st->prefilter_gain, st->prefilter_gain, NULL, 0);
1199             comb_filter(out_mem[c]+st->overlap, out_mem[c]+st->overlap, st->prefilter_period, pitch_index, N-st->overlap, C,
1200                   st->prefilter_gain, gain1, st->mode->window, st->mode->overlap);
1201          } else {
1202             comb_filter(out_mem[c], out_mem[c], st->prefilter_period_old, st->prefilter_period, N, C,
1203                   st->prefilter_gain_old, st->prefilter_gain, st->mode->window, st->mode->overlap);
1204          }
1205       } while (++c<C);
1206 #endif /* ENABLE_POSTFILTER */
1207
1208       deemphasis(out_mem, (celt_word16*)pcm, N, C, st->mode->preemph, st->preemph_memD);
1209       st->prefilter_period_old = st->prefilter_period;
1210       st->prefilter_gain_old = st->prefilter_gain;
1211    }
1212 #endif
1213
1214    st->prefilter_period = pitch_index;
1215    st->prefilter_gain = gain1;
1216
1217    /* If there's any room left (can only happen for very high rates),
1218       fill it with zeros */
1219    while (ec_enc_tell(enc,0) + 8 <= nbCompressedBytes*8)
1220       ec_enc_bits(enc, 0, 8);
1221    ec_enc_done(enc);
1222    
1223    RESTORE_STACK;
1224    if (ec_enc_get_error(enc))
1225       return CELT_CORRUPTED_DATA;
1226    else
1227       return nbCompressedBytes;
1228 }
1229
1230 #ifdef FIXED_POINT
1231 #ifndef DISABLE_FLOAT_API
1232 int celt_encode_with_ec_float(CELTEncoder * restrict st, const float * pcm, int frame_size, unsigned char *compressed, int nbCompressedBytes, ec_enc *enc)
1233 {
1234    int j, ret, C, N, LM, M;
1235    VARDECL(celt_int16, in);
1236    SAVE_STACK;
1237
1238    if (pcm==NULL)
1239       return CELT_BAD_ARG;
1240
1241    for (LM=0;LM<4;LM++)
1242       if (st->mode->shortMdctSize<<LM==frame_size)
1243          break;
1244    if (LM>=MAX_CONFIG_SIZES)
1245       return CELT_BAD_ARG;
1246    M=1<<LM;
1247
1248    C = CHANNELS(st->channels);
1249    N = M*st->mode->shortMdctSize;
1250    ALLOC(in, C*N, celt_int16);
1251
1252    for (j=0;j<C*N;j++)
1253      in[j] = FLOAT2INT16(pcm[j]);
1254
1255    ret=celt_encode_with_ec(st,in,frame_size,compressed,nbCompressedBytes, enc);
1256 #ifdef RESYNTH
1257    for (j=0;j<C*N;j++)
1258       ((float*)pcm)[j]=in[j]*(1.f/32768.f);
1259 #endif
1260    RESTORE_STACK;
1261    return ret;
1262
1263 }
1264 #endif /*DISABLE_FLOAT_API*/
1265 #else
1266 int celt_encode_with_ec(CELTEncoder * restrict st, const celt_int16 * pcm, int frame_size, unsigned char *compressed, int nbCompressedBytes, ec_enc *enc)
1267 {
1268    int j, ret, C, N, LM, M;
1269    VARDECL(celt_sig, in);
1270    SAVE_STACK;
1271
1272    if (pcm==NULL)
1273       return CELT_BAD_ARG;
1274
1275    for (LM=0;LM<4;LM++)
1276       if (st->mode->shortMdctSize<<LM==frame_size)
1277          break;
1278    if (LM>=MAX_CONFIG_SIZES)
1279       return CELT_BAD_ARG;
1280    M=1<<LM;
1281
1282    C=CHANNELS(st->channels);
1283    N=M*st->mode->shortMdctSize;
1284    ALLOC(in, C*N, celt_sig);
1285    for (j=0;j<C*N;j++) {
1286      in[j] = SCALEOUT(pcm[j]);
1287    }
1288
1289    ret = celt_encode_with_ec_float(st,in,frame_size,compressed,nbCompressedBytes, enc);
1290 #ifdef RESYNTH
1291    for (j=0;j<C*N;j++)
1292       ((celt_int16*)pcm)[j] = FLOAT2INT16(in[j]);
1293 #endif
1294    RESTORE_STACK;
1295    return ret;
1296 }
1297 #endif
1298
1299 int celt_encode(CELTEncoder * restrict st, const celt_int16 * pcm, int frame_size, unsigned char *compressed, int nbCompressedBytes)
1300 {
1301    return celt_encode_with_ec(st, pcm, frame_size, compressed, nbCompressedBytes, NULL);
1302 }
1303
1304 #ifndef DISABLE_FLOAT_API
1305 int celt_encode_float(CELTEncoder * restrict st, const float * pcm, int frame_size, unsigned char *compressed, int nbCompressedBytes)
1306 {
1307    return celt_encode_with_ec_float(st, pcm, frame_size, compressed, nbCompressedBytes, NULL);
1308 }
1309 #endif /* DISABLE_FLOAT_API */
1310
1311 int celt_encoder_ctl(CELTEncoder * restrict st, int request, ...)
1312 {
1313    va_list ap;
1314    
1315    va_start(ap, request);
1316    switch (request)
1317    {
1318       case CELT_GET_MODE_REQUEST:
1319       {
1320          const CELTMode ** value = va_arg(ap, const CELTMode**);
1321          if (value==0)
1322             goto bad_arg;
1323          *value=st->mode;
1324       }
1325       break;
1326       case CELT_SET_COMPLEXITY_REQUEST:
1327       {
1328          int value = va_arg(ap, celt_int32);
1329          if (value<0 || value>10)
1330             goto bad_arg;
1331          st->complexity = value;
1332       }
1333       break;
1334       case CELT_SET_START_BAND_REQUEST:
1335       {
1336          celt_int32 value = va_arg(ap, celt_int32);
1337          if (value<0 || value>=st->mode->nbEBands)
1338             goto bad_arg;
1339          st->start = value;
1340       }
1341       break;
1342       case CELT_SET_END_BAND_REQUEST:
1343       {
1344          celt_int32 value = va_arg(ap, celt_int32);
1345          if (value<0 || value>=st->mode->nbEBands)
1346             goto bad_arg;
1347          st->end = value;
1348       }
1349       break;
1350       case CELT_SET_PREDICTION_REQUEST:
1351       {
1352          int value = va_arg(ap, celt_int32);
1353          if (value<0 || value>2)
1354             goto bad_arg;
1355          if (value==0)
1356          {
1357             st->force_intra   = 1;
1358          } else if (value==1) {
1359             st->force_intra   = 0;
1360          } else {
1361             st->force_intra   = 0;
1362          }   
1363       }
1364       break;
1365       case CELT_SET_VBR_CONSTRAINT_REQUEST:
1366       {
1367          celt_int32 value = va_arg(ap, celt_int32);
1368          st->constrained_vbr = value;
1369       }
1370       break;
1371       case CELT_SET_VBR_RATE_REQUEST:
1372       {
1373          celt_int32 value = va_arg(ap, celt_int32);
1374          int frame_rate;
1375          int N = st->mode->shortMdctSize;
1376          if (value<0)
1377             goto bad_arg;
1378          if (value>3072000)
1379             value = 3072000;
1380          frame_rate = ((st->mode->Fs<<3)+(N>>1))/N;
1381          st->vbr_rate_norm = ((value<<(BITRES+3))+(frame_rate>>1))/frame_rate;
1382       }
1383       break;
1384       case CELT_RESET_STATE:
1385       {
1386          CELT_MEMSET((char*)&st->ENCODER_RESET_START, 0,
1387                celt_encoder_get_size(st->mode, st->channels)-
1388                ((char*)&st->ENCODER_RESET_START - (char*)st));
1389          st->vbr_offset = 0;
1390          st->delayedIntra = 1;
1391          st->fold_decision = 1;
1392          st->tonal_average = QCONST16(1.f,8);
1393       }
1394       break;
1395       default:
1396          goto bad_request;
1397    }
1398    va_end(ap);
1399    return CELT_OK;
1400 bad_arg:
1401    va_end(ap);
1402    return CELT_BAD_ARG;
1403 bad_request:
1404    va_end(ap);
1405    return CELT_UNIMPLEMENTED;
1406 }
1407
1408 /**********************************************************************/
1409 /*                                                                    */
1410 /*                             DECODER                                */
1411 /*                                                                    */
1412 /**********************************************************************/
1413 #define DECODE_BUFFER_SIZE 2048
1414
1415 /** Decoder state 
1416  @brief Decoder state
1417  */
1418 struct CELTDecoder {
1419    const CELTMode *mode;
1420    int overlap;
1421    int channels;
1422
1423    int start, end;
1424
1425    /* Everything beyond this point gets cleared on a reset */
1426 #define DECODER_RESET_START last_pitch_index
1427
1428    int last_pitch_index;
1429    int loss_count;
1430    int postfilter_period;
1431    int postfilter_period_old;
1432    celt_word16 postfilter_gain;
1433    celt_word16 postfilter_gain_old;
1434
1435    celt_sig preemph_memD[2];
1436    
1437    celt_sig _decode_mem[1]; /* Size = channels*(DECODE_BUFFER_SIZE+mode->overlap) */
1438    /* celt_word16 lpc[],  Size = channels*LPC_ORDER */
1439    /* celt_word16 oldEBands[], Size = channels*mode->nbEBands */
1440 };
1441
1442 int celt_decoder_get_size(const CELTMode *mode, int channels)
1443 {
1444    int size = sizeof(struct CELTDecoder)
1445             + (channels*(DECODE_BUFFER_SIZE+mode->overlap)-1)*sizeof(celt_sig)
1446             + channels*LPC_ORDER*sizeof(celt_word16)
1447             + channels*mode->nbEBands*sizeof(celt_word16);
1448    return size;
1449 }
1450
1451 CELTDecoder *celt_decoder_create(const CELTMode *mode, int channels, int *error)
1452 {
1453    return celt_decoder_init(
1454          (CELTDecoder *)celt_alloc(celt_decoder_get_size(mode, channels)),
1455          mode, channels, error);
1456 }
1457
1458 CELTDecoder *celt_decoder_init(CELTDecoder *st, const CELTMode *mode, int channels, int *error)
1459 {
1460    if (channels < 0 || channels > 2)
1461    {
1462       if (error)
1463          *error = CELT_BAD_ARG;
1464       return NULL;
1465    }
1466
1467    if (st==NULL)
1468    {
1469       if (error)
1470          *error = CELT_ALLOC_FAIL;
1471       return NULL;
1472    }
1473
1474    CELT_MEMSET((char*)st, 0, celt_decoder_get_size(mode, channels));
1475
1476    st->mode = mode;
1477    st->overlap = mode->overlap;
1478    st->channels = channels;
1479
1480    st->start = 0;
1481    st->end = st->mode->effEBands;
1482
1483    st->loss_count = 0;
1484
1485    if (error)
1486       *error = CELT_OK;
1487    return st;
1488 }
1489
1490 void celt_decoder_destroy(CELTDecoder *st)
1491 {
1492    celt_free(st);
1493 }
1494
1495 static void celt_decode_lost(CELTDecoder * restrict st, celt_word16 * restrict pcm, int N, int LM)
1496 {
1497    int c;
1498    int pitch_index;
1499    int overlap = st->mode->overlap;
1500    celt_word16 fade = Q15ONE;
1501    int i, len;
1502    const int C = CHANNELS(st->channels);
1503    int offset;
1504    celt_sig *out_mem[2];
1505    celt_sig *decode_mem[2];
1506    celt_sig *overlap_mem[2];
1507    celt_word16 *lpc;
1508    SAVE_STACK;
1509    
1510    c=0; do {
1511       decode_mem[c] = st->_decode_mem + c*(DECODE_BUFFER_SIZE+st->overlap);
1512       out_mem[c] = decode_mem[c]+DECODE_BUFFER_SIZE-MAX_PERIOD;
1513       overlap_mem[c] = decode_mem[c]+DECODE_BUFFER_SIZE;
1514    } while (++c<C);
1515    lpc = (celt_word16*)(st->_decode_mem+(DECODE_BUFFER_SIZE+st->overlap)*C);
1516
1517    len = N+st->mode->overlap;
1518    
1519    if (st->loss_count == 0)
1520    {
1521       celt_word16 pitch_buf[MAX_PERIOD>>1];
1522       celt_word32 tmp=0;
1523       celt_word32 mem0[2]={0,0};
1524       celt_word16 mem1[2]={0,0};
1525       int len2 = len;
1526       /* FIXME: This is a kludge */
1527       if (len2>MAX_PERIOD>>1)
1528          len2 = MAX_PERIOD>>1;
1529       pitch_downsample(out_mem, pitch_buf, MAX_PERIOD, MAX_PERIOD,
1530                        C, mem0, mem1);
1531       pitch_search(st->mode, pitch_buf+((MAX_PERIOD-len2)>>1), pitch_buf, len2,
1532                    MAX_PERIOD-len2-100, &pitch_index, &tmp, 1<<LM);
1533       pitch_index = MAX_PERIOD-len2-pitch_index;
1534       st->last_pitch_index = pitch_index;
1535    } else {
1536       pitch_index = st->last_pitch_index;
1537       if (st->loss_count < 5)
1538          fade = QCONST16(.8f,15);
1539       else
1540          fade = 0;
1541    }
1542
1543    c=0; do {
1544       /* FIXME: This is more memory than necessary */
1545       celt_word32 e[2*MAX_PERIOD];
1546       celt_word16 exc[2*MAX_PERIOD];
1547       celt_word32 ac[LPC_ORDER+1];
1548       celt_word16 decay = 1;
1549       celt_word32 S1=0;
1550       celt_word16 mem[LPC_ORDER]={0};
1551
1552       offset = MAX_PERIOD-pitch_index;
1553       for (i=0;i<MAX_PERIOD;i++)
1554          exc[i] = ROUND16(out_mem[c][i], SIG_SHIFT);
1555
1556       if (st->loss_count == 0)
1557       {
1558          _celt_autocorr(exc, ac, st->mode->window, st->mode->overlap,
1559                         LPC_ORDER, MAX_PERIOD);
1560
1561          /* Noise floor -40 dB */
1562 #ifdef FIXED_POINT
1563          ac[0] += SHR32(ac[0],13);
1564 #else
1565          ac[0] *= 1.0001f;
1566 #endif
1567          /* Lag windowing */
1568          for (i=1;i<=LPC_ORDER;i++)
1569          {
1570             /*ac[i] *= exp(-.5*(2*M_PI*.002*i)*(2*M_PI*.002*i));*/
1571 #ifdef FIXED_POINT
1572             ac[i] -= MULT16_32_Q15(2*i*i, ac[i]);
1573 #else
1574             ac[i] -= ac[i]*(.008f*i)*(.008f*i);
1575 #endif
1576          }
1577
1578          _celt_lpc(lpc+c*LPC_ORDER, ac, LPC_ORDER);
1579       }
1580       for (i=0;i<LPC_ORDER;i++)
1581          mem[i] = ROUND16(out_mem[c][MAX_PERIOD-1-i], SIG_SHIFT);
1582       fir(exc, lpc+c*LPC_ORDER, exc, MAX_PERIOD, LPC_ORDER, mem);
1583       /*for (i=0;i<MAX_PERIOD;i++)printf("%d ", exc[i]); printf("\n");*/
1584       /* Check if the waveform is decaying (and if so how fast) */
1585       {
1586          celt_word32 E1=1, E2=1;
1587          int period;
1588          if (pitch_index <= MAX_PERIOD/2)
1589             period = pitch_index;
1590          else
1591             period = MAX_PERIOD/2;
1592          for (i=0;i<period;i++)
1593          {
1594             E1 += SHR32(MULT16_16(exc[MAX_PERIOD-period+i],exc[MAX_PERIOD-period+i]),8);
1595             E2 += SHR32(MULT16_16(exc[MAX_PERIOD-2*period+i],exc[MAX_PERIOD-2*period+i]),8);
1596          }
1597          if (E1 > E2)
1598             E1 = E2;
1599          decay = celt_sqrt(frac_div32(SHR(E1,1),E2));
1600       }
1601
1602       /* Copy excitation, taking decay into account */
1603       for (i=0;i<len+st->mode->overlap;i++)
1604       {
1605          celt_word16 tmp;
1606          if (offset+i >= MAX_PERIOD)
1607          {
1608             offset -= pitch_index;
1609             decay = MULT16_16_Q15(decay, decay);
1610          }
1611          e[i] = SHL32(EXTEND32(MULT16_16_Q15(decay, exc[offset+i])), SIG_SHIFT);
1612          tmp = ROUND16(out_mem[c][offset+i],SIG_SHIFT);
1613          S1 += SHR32(MULT16_16(tmp,tmp),8);
1614       }
1615       for (i=0;i<LPC_ORDER;i++)
1616          mem[i] = ROUND16(out_mem[c][MAX_PERIOD-1-i], SIG_SHIFT);
1617       for (i=0;i<len+st->mode->overlap;i++)
1618          e[i] = MULT16_32_Q15(fade, e[i]);
1619       iir(e, lpc+c*LPC_ORDER, e, len+st->mode->overlap, LPC_ORDER, mem);
1620
1621       {
1622          celt_word32 S2=0;
1623          for (i=0;i<len+overlap;i++)
1624          {
1625             celt_word16 tmp = ROUND16(e[i],SIG_SHIFT);
1626             S2 += SHR32(MULT16_16(tmp,tmp),8);
1627          }
1628          /* This checks for an "explosion" in the synthesis */
1629 #ifdef FIXED_POINT
1630          if (!(S1 > SHR32(S2,2)))
1631 #else
1632          /* Float test is written this way to catch NaNs at the same time */
1633          if (!(S1 > 0.2f*S2))
1634 #endif
1635          {
1636             for (i=0;i<len+overlap;i++)
1637                e[i] = 0;
1638          } else if (S1 < S2)
1639          {
1640             celt_word16 ratio = celt_sqrt(frac_div32(SHR32(S1,1)+1,S2+1));
1641             for (i=0;i<len+overlap;i++)
1642                e[i] = MULT16_32_Q15(ratio, e[i]);
1643          }
1644       }
1645
1646 #ifdef ENABLE_POSTFILTER
1647       /* Apply post-filter to the MDCT overlap of the previous frame */
1648       comb_filter(out_mem[c]+MAX_PERIOD, out_mem[c]+MAX_PERIOD, st->postfilter_period, st->postfilter_period, st->overlap, C,
1649                   st->postfilter_gain, st->postfilter_gain, NULL, 0);
1650 #endif /* ENABLE_POSTFILTER */
1651
1652       for (i=0;i<MAX_PERIOD+st->mode->overlap-N;i++)
1653          out_mem[c][i] = out_mem[c][N+i];
1654
1655       /* Apply TDAC to the concealed audio so that it blends with the
1656          previous and next frames */
1657       for (i=0;i<overlap/2;i++)
1658       {
1659          celt_word32 tmp;
1660          tmp = MULT16_32_Q15(st->mode->window[i],           e[N+overlap-1-i]) +
1661                MULT16_32_Q15(st->mode->window[overlap-i-1], e[N+i          ]);
1662          out_mem[c][MAX_PERIOD+i] = MULT16_32_Q15(st->mode->window[overlap-i-1], tmp);
1663          out_mem[c][MAX_PERIOD+overlap-i-1] = MULT16_32_Q15(st->mode->window[i], tmp);
1664       }
1665       for (i=0;i<N;i++)
1666          out_mem[c][MAX_PERIOD-N+i] = e[i];
1667
1668 #ifdef ENABLE_POSTFILTER
1669       /* Apply pre-filter to the MDCT overlap for the next frame (post-filter will be applied then) */
1670       comb_filter(e, out_mem[c]+MAX_PERIOD, st->postfilter_period, st->postfilter_period, st->overlap, C,
1671                   -st->postfilter_gain, -st->postfilter_gain, NULL, 0);
1672 #endif /* ENABLE_POSTFILTER */
1673       for (i=0;i<overlap;i++)
1674          out_mem[c][MAX_PERIOD+i] = e[i];
1675    } while (++c<C);
1676
1677    {
1678       celt_word32 *out_syn[2];
1679       out_syn[0] = out_mem[0]+MAX_PERIOD-N;
1680       if (C==2)
1681          out_syn[1] = out_mem[1]+MAX_PERIOD-N;
1682       deemphasis(out_syn, pcm, N, C, st->mode->preemph, st->preemph_memD);
1683    }
1684    
1685    st->loss_count++;
1686
1687    RESTORE_STACK;
1688 }
1689
1690 #ifdef FIXED_POINT
1691 int celt_decode_with_ec(CELTDecoder * restrict st, const unsigned char *data, int len, celt_int16 * restrict pcm, int frame_size, ec_dec *dec)
1692 {
1693 #else
1694 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)
1695 {
1696 #endif
1697    int c, i, N;
1698    int has_fold;
1699    int bits;
1700    ec_dec _dec;
1701    ec_byte_buffer buf;
1702    VARDECL(celt_sig, freq);
1703    VARDECL(celt_norm, X);
1704    VARDECL(celt_ener, bandE);
1705    VARDECL(int, fine_quant);
1706    VARDECL(int, pulses);
1707    VARDECL(int, offsets);
1708    VARDECL(int, fine_priority);
1709    VARDECL(int, tf_res);
1710    celt_sig *out_mem[2];
1711    celt_sig *decode_mem[2];
1712    celt_sig *overlap_mem[2];
1713    celt_sig *out_syn[2];
1714    celt_word16 *lpc;
1715    celt_word16 *oldBandE;
1716
1717    int shortBlocks;
1718    int isTransient;
1719    int intra_ener;
1720    const int C = CHANNELS(st->channels);
1721    int LM, M;
1722    int nbFilledBytes, nbAvailableBytes;
1723    int effEnd;
1724    int codedBands;
1725    int alloc_trim;
1726    int postfilter_pitch;
1727    celt_word16 postfilter_gain;
1728    int intensity=0;
1729    int dual_stereo=0;
1730    SAVE_STACK;
1731
1732    if (pcm==NULL)
1733       return CELT_BAD_ARG;
1734
1735    for (LM=0;LM<4;LM++)
1736       if (st->mode->shortMdctSize<<LM==frame_size)
1737          break;
1738    if (LM>=MAX_CONFIG_SIZES)
1739       return CELT_BAD_ARG;
1740    M=1<<LM;
1741
1742    c=0; do {
1743       decode_mem[c] = st->_decode_mem + c*(DECODE_BUFFER_SIZE+st->overlap);
1744       out_mem[c] = decode_mem[c]+DECODE_BUFFER_SIZE-MAX_PERIOD;
1745       overlap_mem[c] = decode_mem[c]+DECODE_BUFFER_SIZE;
1746    } while (++c<C);
1747    lpc = (celt_word16*)(st->_decode_mem+(DECODE_BUFFER_SIZE+st->overlap)*C);
1748    oldBandE = lpc+C*LPC_ORDER;
1749
1750    N = M*st->mode->shortMdctSize;
1751
1752    effEnd = st->end;
1753    if (effEnd > st->mode->effEBands)
1754       effEnd = st->mode->effEBands;
1755
1756    ALLOC(freq, C*N, celt_sig); /**< Interleaved signal MDCTs */
1757    ALLOC(X, C*N, celt_norm);   /**< Interleaved normalised MDCTs */
1758    ALLOC(bandE, st->mode->nbEBands*C, celt_ener);
1759    c=0; do
1760       for (i=0;i<M*st->mode->eBands[st->start];i++)
1761          X[c*N+i] = 0;
1762    while (++c<C);
1763    c=0; do   
1764       for (i=M*st->mode->eBands[effEnd];i<N;i++)
1765          X[c*N+i] = 0;
1766    while (++c<C);
1767
1768    if (data == NULL)
1769    {
1770       celt_decode_lost(st, pcm, N, LM);
1771       RESTORE_STACK;
1772       return CELT_OK;
1773    }
1774    if (len<0) {
1775      RESTORE_STACK;
1776      return CELT_BAD_ARG;
1777    }
1778    
1779    if (dec == NULL)
1780    {
1781       ec_byte_readinit(&buf,(unsigned char*)data,len);
1782       ec_dec_init(&_dec,&buf);
1783       dec = &_dec;
1784       nbFilledBytes = 0;
1785    } else {
1786       nbFilledBytes = (ec_dec_tell(dec, 0)+4)>>3;
1787    }
1788    nbAvailableBytes = len-nbFilledBytes;
1789
1790    if (ec_dec_bit_prob(dec, 32768))
1791    {
1792 #ifdef ENABLE_POSTFILTER
1793       int qg, octave;
1794       octave = ec_dec_uint(dec, 6);
1795       postfilter_pitch = (16<<octave)+ec_dec_bits(dec, 4+octave);
1796       qg = ec_dec_bits(dec, 2);
1797       postfilter_gain = QCONST16(.125f,15)*(qg+2);
1798 #else /* ENABLE_POSTFILTER */
1799       RESTORE_STACK;
1800       return CELT_CORRUPTED_DATA;
1801 #endif /* ENABLE_POSTFILTER */
1802
1803    } else {
1804       postfilter_gain = 0;
1805       postfilter_pitch = 0;
1806    }
1807
1808    /* Decode the global flags (first symbols in the stream) */
1809    intra_ener = ec_dec_bit_prob(dec, 8192);
1810    /* Get band energies */
1811    unquant_coarse_energy(st->mode, st->start, st->end, bandE, oldBandE,
1812          intra_ener, dec, C, LM);
1813
1814    if (LM > 0)
1815       isTransient = ec_dec_bit_prob(dec, 8192);
1816    else
1817       isTransient = 0;
1818
1819    if (isTransient)
1820       shortBlocks = M;
1821    else
1822       shortBlocks = 0;
1823
1824    ALLOC(tf_res, st->mode->nbEBands, int);
1825    tf_decode(st->start, st->end, C, isTransient, tf_res, LM, dec);
1826
1827    has_fold = ec_dec_bit_prob(dec, 8192)<<1;
1828    has_fold |= ec_dec_bit_prob(dec, (has_fold>>1) ? 32768 : 49152);
1829
1830    ALLOC(pulses, st->mode->nbEBands, int);
1831    ALLOC(offsets, st->mode->nbEBands, int);
1832    ALLOC(fine_priority, st->mode->nbEBands, int);
1833
1834    for (i=0;i<st->mode->nbEBands;i++)
1835       offsets[i] = 0;
1836    for (i=0;i<st->mode->nbEBands;i++)
1837    {
1838       if (ec_dec_bit_prob(dec, 1024))
1839       {
1840          while (ec_dec_bit_prob(dec, 32768))
1841             offsets[i]++;
1842          offsets[i]++;
1843          offsets[i] *= (6<<BITRES);
1844       }
1845    }
1846
1847    ALLOC(fine_quant, st->mode->nbEBands, int);
1848    {
1849       int fl;
1850       alloc_trim = 0;
1851       fl = ec_decode_bin(dec, 7);
1852       while (trim_cdf[alloc_trim+1] <= fl)
1853          alloc_trim++;
1854       ec_dec_update(dec, trim_cdf[alloc_trim], trim_cdf[alloc_trim+1], 128);
1855    }
1856
1857    if (C==2)
1858    {
1859       dual_stereo = ec_dec_bit_prob(dec, 32768);
1860       intensity = ec_dec_uint(dec, 1+st->end-st->start);
1861    }
1862
1863    bits = (len*8<<BITRES) - ec_dec_tell(dec, BITRES) - 1;
1864    codedBands = compute_allocation(st->mode, st->start, st->end, offsets,
1865          alloc_trim, bits, pulses, fine_quant, fine_priority, C, LM, dec, 0, 0);
1866    
1867    unquant_fine_energy(st->mode, st->start, st->end, bandE, oldBandE, fine_quant, dec, C);
1868
1869    /* Decode fixed codebook */
1870    quant_all_bands(0, st->mode, st->start, st->end, X, C==2 ? X+N : NULL,
1871          NULL, pulses, shortBlocks, has_fold, dual_stereo, intensity, tf_res, 1,
1872          len*8, dec, LM, codedBands);
1873
1874    unquant_energy_finalise(st->mode, st->start, st->end, bandE, oldBandE,
1875          fine_quant, fine_priority, len*8-ec_dec_tell(dec, 0), dec, C);
1876
1877    log2Amp(st->mode, st->start, st->end, bandE, oldBandE, C);
1878
1879    /* Synthesis */
1880    denormalise_bands(st->mode, X, freq, bandE, effEnd, C, M);
1881
1882    CELT_MOVE(decode_mem[0], decode_mem[0]+N, DECODE_BUFFER_SIZE-N);
1883    if (C==2)
1884       CELT_MOVE(decode_mem[1], decode_mem[1]+N, DECODE_BUFFER_SIZE-N);
1885
1886    c=0; do
1887       for (i=0;i<M*st->mode->eBands[st->start];i++)
1888          freq[c*N+i] = 0;
1889    while (++c<C);
1890    c=0; do
1891       for (i=M*st->mode->eBands[effEnd];i<N;i++)
1892          freq[c*N+i] = 0;
1893    while (++c<C);
1894
1895    out_syn[0] = out_mem[0]+MAX_PERIOD-N;
1896    if (C==2)
1897       out_syn[1] = out_mem[1]+MAX_PERIOD-N;
1898
1899    /* Compute inverse MDCTs */
1900    compute_inv_mdcts(st->mode, shortBlocks, freq, out_syn, overlap_mem, C, LM);
1901
1902 #ifdef ENABLE_POSTFILTER
1903    c=0; do {
1904       st->postfilter_period=IMAX(st->postfilter_period, COMBFILTER_MINPERIOD);
1905       st->postfilter_period_old=IMAX(st->postfilter_period_old, COMBFILTER_MINPERIOD);
1906       if (LM!=0)
1907       {
1908          comb_filter(out_syn[c], out_syn[c], st->postfilter_period, st->postfilter_period, st->overlap, C,
1909                st->postfilter_gain, st->postfilter_gain, NULL, 0);
1910          comb_filter(out_syn[c]+st->overlap, out_syn[c]+st->overlap, st->postfilter_period, postfilter_pitch, N-st->overlap, C,
1911                st->postfilter_gain, postfilter_gain, st->mode->window, st->mode->overlap);
1912       } else {
1913          comb_filter(out_syn[c], out_syn[c], st->postfilter_period_old, st->postfilter_period, N-st->overlap, C,
1914                st->postfilter_gain_old, st->postfilter_gain, st->mode->window, st->mode->overlap);
1915       }
1916    } while (++c<C);
1917    st->postfilter_period_old = st->postfilter_period;
1918    st->postfilter_gain_old = st->postfilter_gain;
1919    st->postfilter_period = postfilter_pitch;
1920    st->postfilter_gain = postfilter_gain;
1921 #endif /* ENABLE_POSTFILTER */
1922
1923    deemphasis(out_syn, pcm, N, C, st->mode->preemph, st->preemph_memD);
1924    st->loss_count = 0;
1925    RESTORE_STACK;
1926    if (ec_dec_get_error(dec))
1927       return CELT_CORRUPTED_DATA;
1928    else
1929       return CELT_OK;
1930 }
1931
1932 #ifdef FIXED_POINT
1933 #ifndef DISABLE_FLOAT_API
1934 int celt_decode_with_ec_float(CELTDecoder * restrict st, const unsigned char *data, int len, float * restrict pcm, int frame_size, ec_dec *dec)
1935 {
1936    int j, ret, C, N, LM, M;
1937    VARDECL(celt_int16, out);
1938    SAVE_STACK;
1939
1940    if (pcm==NULL)
1941       return CELT_BAD_ARG;
1942
1943    for (LM=0;LM<4;LM++)
1944       if (st->mode->shortMdctSize<<LM==frame_size)
1945          break;
1946    if (LM>=MAX_CONFIG_SIZES)
1947       return CELT_BAD_ARG;
1948    M=1<<LM;
1949
1950    C = CHANNELS(st->channels);
1951    N = M*st->mode->shortMdctSize;
1952    
1953    ALLOC(out, C*N, celt_int16);
1954    ret=celt_decode_with_ec(st, data, len, out, frame_size, dec);
1955    if (ret==0)
1956       for (j=0;j<C*N;j++)
1957          pcm[j]=out[j]*(1.f/32768.f);
1958      
1959    RESTORE_STACK;
1960    return ret;
1961 }
1962 #endif /*DISABLE_FLOAT_API*/
1963 #else
1964 int celt_decode_with_ec(CELTDecoder * restrict st, const unsigned char *data, int len, celt_int16 * restrict pcm, int frame_size, ec_dec *dec)
1965 {
1966    int j, ret, C, N, LM, M;
1967    VARDECL(celt_sig, out);
1968    SAVE_STACK;
1969
1970    if (pcm==NULL)
1971       return CELT_BAD_ARG;
1972
1973    for (LM=0;LM<4;LM++)
1974       if (st->mode->shortMdctSize<<LM==frame_size)
1975          break;
1976    if (LM>=MAX_CONFIG_SIZES)
1977       return CELT_BAD_ARG;
1978    M=1<<LM;
1979
1980    C = CHANNELS(st->channels);
1981    N = M*st->mode->shortMdctSize;
1982    ALLOC(out, C*N, celt_sig);
1983
1984    ret=celt_decode_with_ec_float(st, data, len, out, frame_size, dec);
1985
1986    if (ret==0)
1987       for (j=0;j<C*N;j++)
1988          pcm[j] = FLOAT2INT16 (out[j]);
1989    
1990    RESTORE_STACK;
1991    return ret;
1992 }
1993 #endif
1994
1995 int celt_decode(CELTDecoder * restrict st, const unsigned char *data, int len, celt_int16 * restrict pcm, int frame_size)
1996 {
1997    return celt_decode_with_ec(st, data, len, pcm, frame_size, NULL);
1998 }
1999
2000 #ifndef DISABLE_FLOAT_API
2001 int celt_decode_float(CELTDecoder * restrict st, const unsigned char *data, int len, float * restrict pcm, int frame_size)
2002 {
2003    return celt_decode_with_ec_float(st, data, len, pcm, frame_size, NULL);
2004 }
2005 #endif /* DISABLE_FLOAT_API */
2006
2007 int celt_decoder_ctl(CELTDecoder * restrict st, int request, ...)
2008 {
2009    va_list ap;
2010
2011    va_start(ap, request);
2012    switch (request)
2013    {
2014       case CELT_GET_MODE_REQUEST:
2015       {
2016          const CELTMode ** value = va_arg(ap, const CELTMode**);
2017          if (value==0)
2018             goto bad_arg;
2019          *value=st->mode;
2020       }
2021       break;
2022       case CELT_SET_START_BAND_REQUEST:
2023       {
2024          celt_int32 value = va_arg(ap, celt_int32);
2025          if (value<0 || value>=st->mode->nbEBands)
2026             goto bad_arg;
2027          st->start = value;
2028       }
2029       break;
2030       case CELT_SET_END_BAND_REQUEST:
2031       {
2032          celt_int32 value = va_arg(ap, celt_int32);
2033          if (value<0 || value>=st->mode->nbEBands)
2034             goto bad_arg;
2035          st->end = value;
2036       }
2037       break;
2038       case CELT_RESET_STATE:
2039       {
2040          CELT_MEMSET((char*)&st->DECODER_RESET_START, 0,
2041                celt_decoder_get_size(st->mode, st->channels)-
2042                ((char*)&st->DECODER_RESET_START - (char*)st));
2043       }
2044       break;
2045       default:
2046          goto bad_request;
2047    }
2048    va_end(ap);
2049    return CELT_OK;
2050 bad_arg:
2051    va_end(ap);
2052    return CELT_BAD_ARG;
2053 bad_request:
2054       va_end(ap);
2055   return CELT_UNIMPLEMENTED;
2056 }
2057
2058 const char *celt_strerror(int error)
2059 {
2060    static const char *error_strings[8] = {
2061       "success",
2062       "invalid argument",
2063       "invalid mode",
2064       "internal error",
2065       "corrupted stream",
2066       "request not implemented",
2067       "invalid state",
2068       "memory allocation failed"
2069    };
2070    if (error > 0 || error < -7)
2071       return "unknown error";
2072    else 
2073       return error_strings[-error];
2074 }
2075