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