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