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