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