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