Fixes constrained VBR
[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
1168      nbCompressedBytes = IMIN(nbCompressedBytes,nbAvailableBytes+nbFilledBytes);
1169      nbAvailableBytes = nbCompressedBytes - nbFilledBytes;
1170
1171      target=nbAvailableBytes<<(BITRES+3);
1172
1173      if (st->vbr_count < 970)
1174      {
1175         st->vbr_count++;
1176         alpha = celt_rcp(SHL32(EXTEND32(st->vbr_count+20),16));
1177      } else
1178         alpha = QCONST16(.001f,15);
1179      /* How many bits have we used in excess of what we're allowed */
1180      if (st->constrained_vbr)
1181         st->vbr_reservoir += target - vbr_rate;
1182      /*printf ("%d\n", st->vbr_reservoir);*/
1183
1184      /* Compute the offset we need to apply in order to reach the target */
1185      st->vbr_drift += (celt_int32)MULT16_32_Q15(alpha,delta-st->vbr_offset-st->vbr_drift);
1186      st->vbr_offset = -st->vbr_drift;
1187      /*printf ("%d\n", st->vbr_drift);*/
1188
1189      if (st->constrained_vbr && st->vbr_reservoir < 0)
1190      {
1191         /* We're under the min value -- increase rate */
1192         int adjust = (-st->vbr_reservoir)/(8<<BITRES);
1193         nbAvailableBytes += adjust;
1194         st->vbr_reservoir = 0;
1195         /*printf ("+%d\n", adjust);*/
1196      }
1197
1198      /* This moves the raw bits to take into account the new compressed size */
1199      ec_byte_shrink(&buf, nbCompressedBytes);
1200    }
1201
1202    if (C==2)
1203    {
1204       int effectiveRate;
1205
1206       /* Always use MS for 2.5 ms frames until we can do a better analysis */
1207       if (LM!=0)
1208          dual_stereo = stereo_analysis(st->mode, X, st->mode->nbEBands, LM, C, N);
1209
1210       /* Account for coarse energy */
1211       effectiveRate = (8*effectiveBytes - 80)>>LM;
1212
1213       /* effectiveRate in kb/s */
1214       effectiveRate = 2*effectiveRate/5;
1215       if (effectiveRate<35)
1216          intensity = 8;
1217       else if (effectiveRate<50)
1218          intensity = 12;
1219       else if (effectiveRate<68)
1220          intensity = 16;
1221       else if (effectiveRate<84)
1222          intensity = 18;
1223       else if (effectiveRate<102)
1224          intensity = 19;
1225       else if (effectiveRate<130)
1226          intensity = 20;
1227       else
1228          intensity = 100;
1229       intensity = IMIN(st->end,IMAX(st->start, intensity));
1230    }
1231
1232    /* Bit allocation */
1233    ALLOC(fine_quant, st->mode->nbEBands, int);
1234    ALLOC(pulses, st->mode->nbEBands, int);
1235    ALLOC(fine_priority, st->mode->nbEBands, int);
1236
1237    /* bits =   packet size        -       where we are           - safety */
1238    bits = (nbCompressedBytes*8<<BITRES) - ec_enc_tell(enc, BITRES) - 1;
1239    codedBands = compute_allocation(st->mode, st->start, st->end, offsets,
1240          alloc_trim, &intensity, &dual_stereo, bits, pulses, fine_quant,
1241          fine_priority, C, LM, enc, 1, st->lastCodedBands);
1242    st->lastCodedBands = codedBands;
1243
1244    quant_fine_energy(st->mode, st->start, st->end, bandE, oldBandE, error, fine_quant, enc, C);
1245
1246 #ifdef MEASURE_NORM_MSE
1247    float X0[3000];
1248    float bandE0[60];
1249    c=0; do 
1250       for (i=0;i<N;i++)
1251          X0[i+c*N] = X[i+c*N];
1252    while (++c<C);
1253    for (i=0;i<C*st->mode->nbEBands;i++)
1254       bandE0[i] = bandE[i];
1255 #endif
1256
1257    /* Residual quantisation */
1258    quant_all_bands(1, st->mode, st->start, st->end, X, C==2 ? X+N : NULL,
1259          bandE, pulses, shortBlocks, st->spread_decision, dual_stereo, intensity, tf_res, resynth,
1260          nbCompressedBytes*8, enc, LM, codedBands);
1261
1262    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);
1263
1264 #ifdef RESYNTH
1265    /* Re-synthesis of the coded audio if required */
1266    if (resynth)
1267    {
1268       celt_sig *out_mem[2];
1269       celt_sig *overlap_mem[2];
1270
1271       log2Amp(st->mode, st->start, st->end, bandE, oldBandE, C);
1272
1273 #ifdef MEASURE_NORM_MSE
1274       measure_norm_mse(st->mode, X, X0, bandE, bandE0, M, N, C);
1275 #endif
1276
1277       /* Synthesis */
1278       denormalise_bands(st->mode, X, freq, bandE, effEnd, C, M);
1279
1280       CELT_MOVE(st->syn_mem[0], st->syn_mem[0]+N, MAX_PERIOD);
1281       if (C==2)
1282          CELT_MOVE(st->syn_mem[1], st->syn_mem[1]+N, MAX_PERIOD);
1283
1284       c=0; do
1285          for (i=0;i<M*st->mode->eBands[st->start];i++)
1286             freq[c*N+i] = 0;
1287       while (++c<C);
1288       c=0; do
1289          for (i=M*st->mode->eBands[st->end];i<N;i++)
1290             freq[c*N+i] = 0;
1291       while (++c<C);
1292
1293       out_mem[0] = st->syn_mem[0]+MAX_PERIOD;
1294       if (C==2)
1295          out_mem[1] = st->syn_mem[1]+MAX_PERIOD;
1296
1297       c=0; do
1298          overlap_mem[c] = _overlap_mem + c*st->overlap;
1299       while (++c<C);
1300
1301       compute_inv_mdcts(st->mode, shortBlocks, freq, out_mem, overlap_mem, C, LM);
1302
1303 #ifdef ENABLE_POSTFILTER
1304       c=0; do {
1305          st->prefilter_period=IMAX(st->prefilter_period, COMBFILTER_MINPERIOD);
1306          st->prefilter_period_old=IMAX(st->prefilter_period_old, COMBFILTER_MINPERIOD);
1307          if (LM!=0)
1308          {
1309             comb_filter(out_mem[c], out_mem[c], st->prefilter_period, st->prefilter_period, st->overlap, C,
1310                   st->prefilter_gain, st->prefilter_gain, NULL, 0);
1311             comb_filter(out_mem[c]+st->overlap, out_mem[c]+st->overlap, st->prefilter_period, pitch_index, N-st->overlap, C,
1312                   st->prefilter_gain, gain1, st->mode->window, st->mode->overlap);
1313          } else {
1314             comb_filter(out_mem[c], out_mem[c], st->prefilter_period_old, st->prefilter_period, N, C,
1315                   st->prefilter_gain_old, st->prefilter_gain, st->mode->window, st->mode->overlap);
1316          }
1317       } while (++c<C);
1318 #endif /* ENABLE_POSTFILTER */
1319
1320       deemphasis(out_mem, (celt_word16*)pcm, N, C, st->mode->preemph, st->preemph_memD);
1321       st->prefilter_period_old = st->prefilter_period;
1322       st->prefilter_gain_old = st->prefilter_gain;
1323    }
1324 #endif
1325
1326    st->prefilter_period = pitch_index;
1327    st->prefilter_gain = gain1;
1328
1329    /* If there's any room left (can only happen for very high rates),
1330       fill it with zeros */
1331    while (ec_enc_tell(enc,0) + 8 <= nbCompressedBytes*8)
1332       ec_enc_bits(enc, 0, 8);
1333    ec_enc_done(enc);
1334    
1335    RESTORE_STACK;
1336    if (ec_enc_get_error(enc))
1337       return CELT_CORRUPTED_DATA;
1338    else
1339       return nbCompressedBytes;
1340 }
1341
1342 #ifdef FIXED_POINT
1343 #ifndef DISABLE_FLOAT_API
1344 int celt_encode_with_ec_float(CELTEncoder * restrict st, const float * pcm, int frame_size, unsigned char *compressed, int nbCompressedBytes, ec_enc *enc)
1345 {
1346    int j, ret, C, N, LM, M;
1347    VARDECL(celt_int16, in);
1348    SAVE_STACK;
1349
1350    if (pcm==NULL)
1351       return CELT_BAD_ARG;
1352
1353    for (LM=0;LM<4;LM++)
1354       if (st->mode->shortMdctSize<<LM==frame_size)
1355          break;
1356    if (LM>=MAX_CONFIG_SIZES)
1357       return CELT_BAD_ARG;
1358    M=1<<LM;
1359
1360    C = CHANNELS(st->channels);
1361    N = M*st->mode->shortMdctSize;
1362    ALLOC(in, C*N, celt_int16);
1363
1364    for (j=0;j<C*N;j++)
1365      in[j] = FLOAT2INT16(pcm[j]);
1366
1367    ret=celt_encode_with_ec(st,in,frame_size,compressed,nbCompressedBytes, enc);
1368 #ifdef RESYNTH
1369    for (j=0;j<C*N;j++)
1370       ((float*)pcm)[j]=in[j]*(1.f/32768.f);
1371 #endif
1372    RESTORE_STACK;
1373    return ret;
1374
1375 }
1376 #endif /*DISABLE_FLOAT_API*/
1377 #else
1378 int celt_encode_with_ec(CELTEncoder * restrict st, const celt_int16 * pcm, int frame_size, unsigned char *compressed, int nbCompressedBytes, ec_enc *enc)
1379 {
1380    int j, ret, C, N, LM, M;
1381    VARDECL(celt_sig, in);
1382    SAVE_STACK;
1383
1384    if (pcm==NULL)
1385       return CELT_BAD_ARG;
1386
1387    for (LM=0;LM<4;LM++)
1388       if (st->mode->shortMdctSize<<LM==frame_size)
1389          break;
1390    if (LM>=MAX_CONFIG_SIZES)
1391       return CELT_BAD_ARG;
1392    M=1<<LM;
1393
1394    C=CHANNELS(st->channels);
1395    N=M*st->mode->shortMdctSize;
1396    ALLOC(in, C*N, celt_sig);
1397    for (j=0;j<C*N;j++) {
1398      in[j] = SCALEOUT(pcm[j]);
1399    }
1400
1401    ret = celt_encode_with_ec_float(st,in,frame_size,compressed,nbCompressedBytes, enc);
1402 #ifdef RESYNTH
1403    for (j=0;j<C*N;j++)
1404       ((celt_int16*)pcm)[j] = FLOAT2INT16(in[j]);
1405 #endif
1406    RESTORE_STACK;
1407    return ret;
1408 }
1409 #endif
1410
1411 int celt_encode(CELTEncoder * restrict st, const celt_int16 * pcm, int frame_size, unsigned char *compressed, int nbCompressedBytes)
1412 {
1413    return celt_encode_with_ec(st, pcm, frame_size, compressed, nbCompressedBytes, NULL);
1414 }
1415
1416 #ifndef DISABLE_FLOAT_API
1417 int celt_encode_float(CELTEncoder * restrict st, const float * pcm, int frame_size, unsigned char *compressed, int nbCompressedBytes)
1418 {
1419    return celt_encode_with_ec_float(st, pcm, frame_size, compressed, nbCompressedBytes, NULL);
1420 }
1421 #endif /* DISABLE_FLOAT_API */
1422
1423 int celt_encoder_ctl(CELTEncoder * restrict st, int request, ...)
1424 {
1425    va_list ap;
1426    
1427    va_start(ap, request);
1428    switch (request)
1429    {
1430       case CELT_GET_MODE_REQUEST:
1431       {
1432          const CELTMode ** value = va_arg(ap, const CELTMode**);
1433          if (value==0)
1434             goto bad_arg;
1435          *value=st->mode;
1436       }
1437       break;
1438       case CELT_SET_COMPLEXITY_REQUEST:
1439       {
1440          int value = va_arg(ap, celt_int32);
1441          if (value<0 || value>10)
1442             goto bad_arg;
1443          st->complexity = value;
1444       }
1445       break;
1446       case CELT_SET_START_BAND_REQUEST:
1447       {
1448          celt_int32 value = va_arg(ap, celt_int32);
1449          if (value<0 || value>=st->mode->nbEBands)
1450             goto bad_arg;
1451          st->start = value;
1452       }
1453       break;
1454       case CELT_SET_END_BAND_REQUEST:
1455       {
1456          celt_int32 value = va_arg(ap, celt_int32);
1457          if (value<1 || value>st->mode->nbEBands)
1458             goto bad_arg;
1459          st->end = value;
1460       }
1461       break;
1462       case CELT_SET_PREDICTION_REQUEST:
1463       {
1464          int value = va_arg(ap, celt_int32);
1465          if (value<0 || value>2)
1466             goto bad_arg;
1467          if (value==0)
1468          {
1469             st->force_intra   = 1;
1470          } else if (value==1) {
1471             st->force_intra   = 0;
1472          } else {
1473             st->force_intra   = 0;
1474          }   
1475       }
1476       break;
1477       case CELT_SET_VBR_CONSTRAINT_REQUEST:
1478       {
1479          celt_int32 value = va_arg(ap, celt_int32);
1480          st->constrained_vbr = value;
1481       }
1482       break;
1483       case CELT_SET_VBR_RATE_REQUEST:
1484       {
1485          celt_int32 value = va_arg(ap, celt_int32);
1486          int frame_rate;
1487          int N = st->mode->shortMdctSize;
1488          if (value<0)
1489             goto bad_arg;
1490          if (value>3072000)
1491             value = 3072000;
1492          frame_rate = ((st->mode->Fs<<3)+(N>>1))/N;
1493          st->vbr_rate_norm = ((value<<(BITRES+3))+(frame_rate>>1))/frame_rate;
1494       }
1495       break;
1496       case CELT_RESET_STATE:
1497       {
1498          CELT_MEMSET((char*)&st->ENCODER_RESET_START, 0,
1499                celt_encoder_get_size(st->mode, st->channels)-
1500                ((char*)&st->ENCODER_RESET_START - (char*)st));
1501          st->vbr_offset = 0;
1502          st->delayedIntra = 1;
1503          st->spread_decision = SPREAD_NORMAL;
1504          st->tonal_average = QCONST16(1.f,8);
1505       }
1506       break;
1507       default:
1508          goto bad_request;
1509    }
1510    va_end(ap);
1511    return CELT_OK;
1512 bad_arg:
1513    va_end(ap);
1514    return CELT_BAD_ARG;
1515 bad_request:
1516    va_end(ap);
1517    return CELT_UNIMPLEMENTED;
1518 }
1519
1520 /**********************************************************************/
1521 /*                                                                    */
1522 /*                             DECODER                                */
1523 /*                                                                    */
1524 /**********************************************************************/
1525 #define DECODE_BUFFER_SIZE 2048
1526
1527 /** Decoder state 
1528  @brief Decoder state
1529  */
1530 struct CELTDecoder {
1531    const CELTMode *mode;
1532    int overlap;
1533    int channels;
1534
1535    int start, end;
1536
1537    /* Everything beyond this point gets cleared on a reset */
1538 #define DECODER_RESET_START last_pitch_index
1539
1540    int last_pitch_index;
1541    int loss_count;
1542    int postfilter_period;
1543    int postfilter_period_old;
1544    celt_word16 postfilter_gain;
1545    celt_word16 postfilter_gain_old;
1546
1547    celt_sig preemph_memD[2];
1548    
1549    celt_sig _decode_mem[1]; /* Size = channels*(DECODE_BUFFER_SIZE+mode->overlap) */
1550    /* celt_word16 lpc[],  Size = channels*LPC_ORDER */
1551    /* celt_word16 oldEBands[], Size = channels*mode->nbEBands */
1552 };
1553
1554 int celt_decoder_get_size(const CELTMode *mode, int channels)
1555 {
1556    int size = sizeof(struct CELTDecoder)
1557             + (channels*(DECODE_BUFFER_SIZE+mode->overlap)-1)*sizeof(celt_sig)
1558             + channels*LPC_ORDER*sizeof(celt_word16)
1559             + channels*mode->nbEBands*sizeof(celt_word16);
1560    return size;
1561 }
1562
1563 CELTDecoder *celt_decoder_create(const CELTMode *mode, int channels, int *error)
1564 {
1565    return celt_decoder_init(
1566          (CELTDecoder *)celt_alloc(celt_decoder_get_size(mode, channels)),
1567          mode, channels, error);
1568 }
1569
1570 CELTDecoder *celt_decoder_init(CELTDecoder *st, const CELTMode *mode, int channels, int *error)
1571 {
1572    if (channels < 0 || channels > 2)
1573    {
1574       if (error)
1575          *error = CELT_BAD_ARG;
1576       return NULL;
1577    }
1578
1579    if (st==NULL)
1580    {
1581       if (error)
1582          *error = CELT_ALLOC_FAIL;
1583       return NULL;
1584    }
1585
1586    CELT_MEMSET((char*)st, 0, celt_decoder_get_size(mode, channels));
1587
1588    st->mode = mode;
1589    st->overlap = mode->overlap;
1590    st->channels = channels;
1591
1592    st->start = 0;
1593    st->end = st->mode->effEBands;
1594
1595    st->loss_count = 0;
1596
1597    if (error)
1598       *error = CELT_OK;
1599    return st;
1600 }
1601
1602 void celt_decoder_destroy(CELTDecoder *st)
1603 {
1604    celt_free(st);
1605 }
1606
1607 static void celt_decode_lost(CELTDecoder * restrict st, celt_word16 * restrict pcm, int N, int LM)
1608 {
1609    int c;
1610    int pitch_index;
1611    int overlap = st->mode->overlap;
1612    celt_word16 fade = Q15ONE;
1613    int i, len;
1614    const int C = CHANNELS(st->channels);
1615    int offset;
1616    celt_sig *out_mem[2];
1617    celt_sig *decode_mem[2];
1618    celt_sig *overlap_mem[2];
1619    celt_word16 *lpc;
1620    SAVE_STACK;
1621    
1622    c=0; do {
1623       decode_mem[c] = st->_decode_mem + c*(DECODE_BUFFER_SIZE+st->overlap);
1624       out_mem[c] = decode_mem[c]+DECODE_BUFFER_SIZE-MAX_PERIOD;
1625       overlap_mem[c] = decode_mem[c]+DECODE_BUFFER_SIZE;
1626    } while (++c<C);
1627    lpc = (celt_word16*)(st->_decode_mem+(DECODE_BUFFER_SIZE+st->overlap)*C);
1628
1629    len = N+st->mode->overlap;
1630    
1631    if (st->loss_count == 0)
1632    {
1633       celt_word16 pitch_buf[MAX_PERIOD>>1];
1634       celt_word32 tmp=0;
1635       celt_word32 mem0[2]={0,0};
1636       celt_word16 mem1[2]={0,0};
1637       int len2 = len;
1638       /* FIXME: This is a kludge */
1639       if (len2>MAX_PERIOD>>1)
1640          len2 = MAX_PERIOD>>1;
1641       pitch_downsample(out_mem, pitch_buf, MAX_PERIOD, MAX_PERIOD,
1642                        C, mem0, mem1);
1643       pitch_search(st->mode, pitch_buf+((MAX_PERIOD-len2)>>1), pitch_buf, len2,
1644                    MAX_PERIOD-len2-100, &pitch_index, &tmp, 1<<LM);
1645       pitch_index = MAX_PERIOD-len2-pitch_index;
1646       st->last_pitch_index = pitch_index;
1647    } else {
1648       pitch_index = st->last_pitch_index;
1649       if (st->loss_count < 5)
1650          fade = QCONST16(.8f,15);
1651       else
1652          fade = 0;
1653    }
1654
1655    c=0; do {
1656       /* FIXME: This is more memory than necessary */
1657       celt_word32 e[2*MAX_PERIOD];
1658       celt_word16 exc[2*MAX_PERIOD];
1659       celt_word32 ac[LPC_ORDER+1];
1660       celt_word16 decay = 1;
1661       celt_word32 S1=0;
1662       celt_word16 mem[LPC_ORDER]={0};
1663
1664       offset = MAX_PERIOD-pitch_index;
1665       for (i=0;i<MAX_PERIOD;i++)
1666          exc[i] = ROUND16(out_mem[c][i], SIG_SHIFT);
1667
1668       if (st->loss_count == 0)
1669       {
1670          _celt_autocorr(exc, ac, st->mode->window, st->mode->overlap,
1671                         LPC_ORDER, MAX_PERIOD);
1672
1673          /* Noise floor -40 dB */
1674 #ifdef FIXED_POINT
1675          ac[0] += SHR32(ac[0],13);
1676 #else
1677          ac[0] *= 1.0001f;
1678 #endif
1679          /* Lag windowing */
1680          for (i=1;i<=LPC_ORDER;i++)
1681          {
1682             /*ac[i] *= exp(-.5*(2*M_PI*.002*i)*(2*M_PI*.002*i));*/
1683 #ifdef FIXED_POINT
1684             ac[i] -= MULT16_32_Q15(2*i*i, ac[i]);
1685 #else
1686             ac[i] -= ac[i]*(.008f*i)*(.008f*i);
1687 #endif
1688          }
1689
1690          _celt_lpc(lpc+c*LPC_ORDER, ac, LPC_ORDER);
1691       }
1692       for (i=0;i<LPC_ORDER;i++)
1693          mem[i] = ROUND16(out_mem[c][MAX_PERIOD-1-i], SIG_SHIFT);
1694       fir(exc, lpc+c*LPC_ORDER, exc, MAX_PERIOD, LPC_ORDER, mem);
1695       /*for (i=0;i<MAX_PERIOD;i++)printf("%d ", exc[i]); printf("\n");*/
1696       /* Check if the waveform is decaying (and if so how fast) */
1697       {
1698          celt_word32 E1=1, E2=1;
1699          int period;
1700          if (pitch_index <= MAX_PERIOD/2)
1701             period = pitch_index;
1702          else
1703             period = MAX_PERIOD/2;
1704          for (i=0;i<period;i++)
1705          {
1706             E1 += SHR32(MULT16_16(exc[MAX_PERIOD-period+i],exc[MAX_PERIOD-period+i]),8);
1707             E2 += SHR32(MULT16_16(exc[MAX_PERIOD-2*period+i],exc[MAX_PERIOD-2*period+i]),8);
1708          }
1709          if (E1 > E2)
1710             E1 = E2;
1711          decay = celt_sqrt(frac_div32(SHR(E1,1),E2));
1712       }
1713
1714       /* Copy excitation, taking decay into account */
1715       for (i=0;i<len+st->mode->overlap;i++)
1716       {
1717          celt_word16 tmp;
1718          if (offset+i >= MAX_PERIOD)
1719          {
1720             offset -= pitch_index;
1721             decay = MULT16_16_Q15(decay, decay);
1722          }
1723          e[i] = SHL32(EXTEND32(MULT16_16_Q15(decay, exc[offset+i])), SIG_SHIFT);
1724          tmp = ROUND16(out_mem[c][offset+i],SIG_SHIFT);
1725          S1 += SHR32(MULT16_16(tmp,tmp),8);
1726       }
1727       for (i=0;i<LPC_ORDER;i++)
1728          mem[i] = ROUND16(out_mem[c][MAX_PERIOD-1-i], SIG_SHIFT);
1729       for (i=0;i<len+st->mode->overlap;i++)
1730          e[i] = MULT16_32_Q15(fade, e[i]);
1731       iir(e, lpc+c*LPC_ORDER, e, len+st->mode->overlap, LPC_ORDER, mem);
1732
1733       {
1734          celt_word32 S2=0;
1735          for (i=0;i<len+overlap;i++)
1736          {
1737             celt_word16 tmp = ROUND16(e[i],SIG_SHIFT);
1738             S2 += SHR32(MULT16_16(tmp,tmp),8);
1739          }
1740          /* This checks for an "explosion" in the synthesis */
1741 #ifdef FIXED_POINT
1742          if (!(S1 > SHR32(S2,2)))
1743 #else
1744          /* Float test is written this way to catch NaNs at the same time */
1745          if (!(S1 > 0.2f*S2))
1746 #endif
1747          {
1748             for (i=0;i<len+overlap;i++)
1749                e[i] = 0;
1750          } else if (S1 < S2)
1751          {
1752             celt_word16 ratio = celt_sqrt(frac_div32(SHR32(S1,1)+1,S2+1));
1753             for (i=0;i<len+overlap;i++)
1754                e[i] = MULT16_32_Q15(ratio, e[i]);
1755          }
1756       }
1757
1758 #ifdef ENABLE_POSTFILTER
1759       /* Apply post-filter to the MDCT overlap of the previous frame */
1760       comb_filter(out_mem[c]+MAX_PERIOD, out_mem[c]+MAX_PERIOD, st->postfilter_period, st->postfilter_period, st->overlap, C,
1761                   st->postfilter_gain, st->postfilter_gain, NULL, 0);
1762 #endif /* ENABLE_POSTFILTER */
1763
1764       for (i=0;i<MAX_PERIOD+st->mode->overlap-N;i++)
1765          out_mem[c][i] = out_mem[c][N+i];
1766
1767       /* Apply TDAC to the concealed audio so that it blends with the
1768          previous and next frames */
1769       for (i=0;i<overlap/2;i++)
1770       {
1771          celt_word32 tmp;
1772          tmp = MULT16_32_Q15(st->mode->window[i],           e[N+overlap-1-i]) +
1773                MULT16_32_Q15(st->mode->window[overlap-i-1], e[N+i          ]);
1774          out_mem[c][MAX_PERIOD+i] = MULT16_32_Q15(st->mode->window[overlap-i-1], tmp);
1775          out_mem[c][MAX_PERIOD+overlap-i-1] = MULT16_32_Q15(st->mode->window[i], tmp);
1776       }
1777       for (i=0;i<N;i++)
1778          out_mem[c][MAX_PERIOD-N+i] = e[i];
1779
1780 #ifdef ENABLE_POSTFILTER
1781       /* Apply pre-filter to the MDCT overlap for the next frame (post-filter will be applied then) */
1782       comb_filter(e, out_mem[c]+MAX_PERIOD, st->postfilter_period, st->postfilter_period, st->overlap, C,
1783                   -st->postfilter_gain, -st->postfilter_gain, NULL, 0);
1784 #endif /* ENABLE_POSTFILTER */
1785       for (i=0;i<overlap;i++)
1786          out_mem[c][MAX_PERIOD+i] = e[i];
1787    } while (++c<C);
1788
1789    {
1790       celt_word32 *out_syn[2];
1791       out_syn[0] = out_mem[0]+MAX_PERIOD-N;
1792       if (C==2)
1793          out_syn[1] = out_mem[1]+MAX_PERIOD-N;
1794       deemphasis(out_syn, pcm, N, C, st->mode->preemph, st->preemph_memD);
1795    }
1796    
1797    st->loss_count++;
1798
1799    RESTORE_STACK;
1800 }
1801
1802 #ifdef FIXED_POINT
1803 int celt_decode_with_ec(CELTDecoder * restrict st, const unsigned char *data, int len, celt_int16 * restrict pcm, int frame_size, ec_dec *dec)
1804 {
1805 #else
1806 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)
1807 {
1808 #endif
1809    int c, i, N;
1810    int spread_decision;
1811    int bits;
1812    ec_dec _dec;
1813    ec_byte_buffer buf;
1814    VARDECL(celt_sig, freq);
1815    VARDECL(celt_norm, X);
1816    VARDECL(celt_ener, bandE);
1817    VARDECL(int, fine_quant);
1818    VARDECL(int, pulses);
1819    VARDECL(int, offsets);
1820    VARDECL(int, fine_priority);
1821    VARDECL(int, tf_res);
1822    celt_sig *out_mem[2];
1823    celt_sig *decode_mem[2];
1824    celt_sig *overlap_mem[2];
1825    celt_sig *out_syn[2];
1826    celt_word16 *lpc;
1827    celt_word16 *oldBandE;
1828
1829    int shortBlocks;
1830    int isTransient;
1831    int intra_ener;
1832    const int C = CHANNELS(st->channels);
1833    int LM, M;
1834    int effEnd;
1835    int codedBands;
1836    int alloc_trim;
1837    int postfilter_pitch;
1838    celt_word16 postfilter_gain;
1839    int intensity=0;
1840    int dual_stereo=0;
1841    celt_int32 total_bits;
1842    celt_int32 tell;
1843    int dynalloc_logp;
1844    SAVE_STACK;
1845
1846    if (pcm==NULL)
1847       return CELT_BAD_ARG;
1848
1849    for (LM=0;LM<4;LM++)
1850       if (st->mode->shortMdctSize<<LM==frame_size)
1851          break;
1852    if (LM>=MAX_CONFIG_SIZES)
1853       return CELT_BAD_ARG;
1854    M=1<<LM;
1855
1856    c=0; do {
1857       decode_mem[c] = st->_decode_mem + c*(DECODE_BUFFER_SIZE+st->overlap);
1858       out_mem[c] = decode_mem[c]+DECODE_BUFFER_SIZE-MAX_PERIOD;
1859       overlap_mem[c] = decode_mem[c]+DECODE_BUFFER_SIZE;
1860    } while (++c<C);
1861    lpc = (celt_word16*)(st->_decode_mem+(DECODE_BUFFER_SIZE+st->overlap)*C);
1862    oldBandE = lpc+C*LPC_ORDER;
1863
1864    N = M*st->mode->shortMdctSize;
1865
1866    effEnd = st->end;
1867    if (effEnd > st->mode->effEBands)
1868       effEnd = st->mode->effEBands;
1869
1870    ALLOC(freq, C*N, celt_sig); /**< Interleaved signal MDCTs */
1871    ALLOC(X, C*N, celt_norm);   /**< Interleaved normalised MDCTs */
1872    ALLOC(bandE, st->mode->nbEBands*C, celt_ener);
1873    c=0; do
1874       for (i=0;i<M*st->mode->eBands[st->start];i++)
1875          X[c*N+i] = 0;
1876    while (++c<C);
1877    c=0; do   
1878       for (i=M*st->mode->eBands[effEnd];i<N;i++)
1879          X[c*N+i] = 0;
1880    while (++c<C);
1881
1882    if (data == NULL)
1883    {
1884       celt_decode_lost(st, pcm, N, LM);
1885       RESTORE_STACK;
1886       return CELT_OK;
1887    }
1888    if (len<0) {
1889      RESTORE_STACK;
1890      return CELT_BAD_ARG;
1891    }
1892    
1893    if (dec == NULL)
1894    {
1895       ec_byte_readinit(&buf,(unsigned char*)data,len);
1896       ec_dec_init(&_dec,&buf);
1897       dec = &_dec;
1898    }
1899
1900    total_bits = len*8;
1901    tell = ec_dec_tell(dec, 0);
1902
1903    postfilter_gain = 0;
1904    postfilter_pitch = 0;
1905    if (tell+15 <= total_bits)
1906    {
1907       if(ec_dec_bit_logp(dec, 1))
1908       {
1909 #ifdef ENABLE_POSTFILTER
1910          int qg, octave;
1911          octave = ec_dec_uint(dec, 6);
1912          postfilter_pitch = (16<<octave)+ec_dec_bits(dec, 4+octave);
1913          qg = ec_dec_bits(dec, 2);
1914          postfilter_gain = QCONST16(.125f,15)*(qg+2);
1915          tell = ec_dec_tell(dec, 0);
1916 #else /* ENABLE_POSTFILTER */
1917          RESTORE_STACK;
1918          return CELT_CORRUPTED_DATA;
1919 #endif /* ENABLE_POSTFILTER */
1920       }
1921       tell = ec_dec_tell(dec, 0);
1922    }
1923
1924    if (LM > 0 && tell+3 <= total_bits)
1925    {
1926       isTransient = ec_dec_bit_logp(dec, 3);
1927       tell = ec_dec_tell(dec, 0);
1928    }
1929    else
1930       isTransient = 0;
1931
1932    if (isTransient)
1933       shortBlocks = M;
1934    else
1935       shortBlocks = 0;
1936
1937    /* Decode the global flags (first symbols in the stream) */
1938    intra_ener = tell+3<=total_bits ? ec_dec_bit_logp(dec, 3) : 0;
1939    /* Get band energies */
1940    unquant_coarse_energy(st->mode, st->start, st->end, bandE, oldBandE,
1941          intra_ener, dec, C, LM);
1942
1943    ALLOC(tf_res, st->mode->nbEBands, int);
1944    tf_decode(st->start, st->end, C, isTransient, tf_res, LM, dec);
1945
1946    tell = ec_dec_tell(dec, 0);
1947    spread_decision = SPREAD_NORMAL;
1948    if (tell+4 <= total_bits)
1949       spread_decision = ec_dec_icdf(dec, spread_icdf, 5);
1950
1951    ALLOC(pulses, st->mode->nbEBands, int);
1952    ALLOC(offsets, st->mode->nbEBands, int);
1953    ALLOC(fine_priority, st->mode->nbEBands, int);
1954
1955    dynalloc_logp = 6;
1956    total_bits<<=BITRES;
1957    tell = ec_dec_tell(dec, BITRES);
1958    for (i=st->start;i<st->end;i++)
1959    {
1960       int width, quanta;
1961       int dynalloc_loop_logp;
1962       int boost;
1963       width = C*(st->mode->eBands[i+1]-st->mode->eBands[i])<<LM;
1964       /* quanta is 6 bits, but no more than 1 bit/sample
1965          and no less than 1/8 bit/sample */
1966       quanta = IMIN(width<<BITRES, IMAX(6<<BITRES, width));
1967       dynalloc_loop_logp = dynalloc_logp;
1968       boost = 0;
1969       while (tell+(dynalloc_loop_logp<<BITRES) < total_bits &&
1970             boost < (64<<LM)*(C<<BITRES))
1971       {
1972          int flag;
1973          flag = ec_dec_bit_logp(dec, dynalloc_loop_logp);
1974          tell = ec_dec_tell(dec, BITRES);
1975          if (!flag)
1976             break;
1977          boost += quanta;
1978          total_bits -= quanta;
1979          dynalloc_loop_logp = 1;
1980       }
1981       offsets[i] = boost;
1982       /* Making dynalloc more likely */
1983       if (boost>0)
1984          dynalloc_logp = IMAX(2, dynalloc_logp-1);
1985    }
1986
1987    ALLOC(fine_quant, st->mode->nbEBands, int);
1988    alloc_trim = tell+(6<<BITRES) <= total_bits ?
1989          ec_dec_icdf(dec, trim_icdf, 7) : 5;
1990
1991    bits = (len*8<<BITRES) - ec_dec_tell(dec, BITRES) - 1;
1992    codedBands = compute_allocation(st->mode, st->start, st->end, offsets,
1993          alloc_trim, &intensity, &dual_stereo, bits, pulses, fine_quant,
1994          fine_priority, C, LM, dec, 0, 0);
1995    
1996    unquant_fine_energy(st->mode, st->start, st->end, bandE, oldBandE, fine_quant, dec, C);
1997
1998    /* Decode fixed codebook */
1999    quant_all_bands(0, st->mode, st->start, st->end, X, C==2 ? X+N : NULL,
2000          NULL, pulses, shortBlocks, spread_decision, dual_stereo, intensity, tf_res, 1,
2001          len*8, dec, LM, codedBands);
2002
2003    unquant_energy_finalise(st->mode, st->start, st->end, bandE, oldBandE,
2004          fine_quant, fine_priority, len*8-ec_dec_tell(dec, 0), dec, C);
2005
2006    log2Amp(st->mode, st->start, st->end, bandE, oldBandE, C);
2007
2008    /* Synthesis */
2009    denormalise_bands(st->mode, X, freq, bandE, effEnd, C, M);
2010
2011    CELT_MOVE(decode_mem[0], decode_mem[0]+N, DECODE_BUFFER_SIZE-N);
2012    if (C==2)
2013       CELT_MOVE(decode_mem[1], decode_mem[1]+N, DECODE_BUFFER_SIZE-N);
2014
2015    c=0; do
2016       for (i=0;i<M*st->mode->eBands[st->start];i++)
2017          freq[c*N+i] = 0;
2018    while (++c<C);
2019    c=0; do
2020       for (i=M*st->mode->eBands[effEnd];i<N;i++)
2021          freq[c*N+i] = 0;
2022    while (++c<C);
2023
2024    out_syn[0] = out_mem[0]+MAX_PERIOD-N;
2025    if (C==2)
2026       out_syn[1] = out_mem[1]+MAX_PERIOD-N;
2027
2028    /* Compute inverse MDCTs */
2029    compute_inv_mdcts(st->mode, shortBlocks, freq, out_syn, overlap_mem, C, LM);
2030
2031 #ifdef ENABLE_POSTFILTER
2032    c=0; do {
2033       st->postfilter_period=IMAX(st->postfilter_period, COMBFILTER_MINPERIOD);
2034       st->postfilter_period_old=IMAX(st->postfilter_period_old, COMBFILTER_MINPERIOD);
2035       if (LM!=0)
2036       {
2037          comb_filter(out_syn[c], out_syn[c], st->postfilter_period, st->postfilter_period, st->overlap, C,
2038                st->postfilter_gain, st->postfilter_gain, NULL, 0);
2039          comb_filter(out_syn[c]+st->overlap, out_syn[c]+st->overlap, st->postfilter_period, postfilter_pitch, N-st->overlap, C,
2040                st->postfilter_gain, postfilter_gain, st->mode->window, st->mode->overlap);
2041       } else {
2042          comb_filter(out_syn[c], out_syn[c], st->postfilter_period_old, st->postfilter_period, N-st->overlap, C,
2043                st->postfilter_gain_old, st->postfilter_gain, st->mode->window, st->mode->overlap);
2044       }
2045    } while (++c<C);
2046    st->postfilter_period_old = st->postfilter_period;
2047    st->postfilter_gain_old = st->postfilter_gain;
2048    st->postfilter_period = postfilter_pitch;
2049    st->postfilter_gain = postfilter_gain;
2050 #endif /* ENABLE_POSTFILTER */
2051
2052    deemphasis(out_syn, pcm, N, C, st->mode->preemph, st->preemph_memD);
2053    st->loss_count = 0;
2054    RESTORE_STACK;
2055    if (ec_dec_tell(dec,0) > 8*len || ec_dec_get_error(dec))
2056       return CELT_CORRUPTED_DATA;
2057    else
2058       return CELT_OK;
2059 }
2060
2061 #ifdef FIXED_POINT
2062 #ifndef DISABLE_FLOAT_API
2063 int celt_decode_with_ec_float(CELTDecoder * restrict st, const unsigned char *data, int len, float * restrict pcm, int frame_size, ec_dec *dec)
2064 {
2065    int j, ret, C, N, LM, M;
2066    VARDECL(celt_int16, out);
2067    SAVE_STACK;
2068
2069    if (pcm==NULL)
2070       return CELT_BAD_ARG;
2071
2072    for (LM=0;LM<4;LM++)
2073       if (st->mode->shortMdctSize<<LM==frame_size)
2074          break;
2075    if (LM>=MAX_CONFIG_SIZES)
2076       return CELT_BAD_ARG;
2077    M=1<<LM;
2078
2079    C = CHANNELS(st->channels);
2080    N = M*st->mode->shortMdctSize;
2081    
2082    ALLOC(out, C*N, celt_int16);
2083    ret=celt_decode_with_ec(st, data, len, out, frame_size, dec);
2084    if (ret==0)
2085       for (j=0;j<C*N;j++)
2086          pcm[j]=out[j]*(1.f/32768.f);
2087      
2088    RESTORE_STACK;
2089    return ret;
2090 }
2091 #endif /*DISABLE_FLOAT_API*/
2092 #else
2093 int celt_decode_with_ec(CELTDecoder * restrict st, const unsigned char *data, int len, celt_int16 * restrict pcm, int frame_size, ec_dec *dec)
2094 {
2095    int j, ret, C, N, LM, M;
2096    VARDECL(celt_sig, out);
2097    SAVE_STACK;
2098
2099    if (pcm==NULL)
2100       return CELT_BAD_ARG;
2101
2102    for (LM=0;LM<4;LM++)
2103       if (st->mode->shortMdctSize<<LM==frame_size)
2104          break;
2105    if (LM>=MAX_CONFIG_SIZES)
2106       return CELT_BAD_ARG;
2107    M=1<<LM;
2108
2109    C = CHANNELS(st->channels);
2110    N = M*st->mode->shortMdctSize;
2111    ALLOC(out, C*N, celt_sig);
2112
2113    ret=celt_decode_with_ec_float(st, data, len, out, frame_size, dec);
2114
2115    if (ret==0)
2116       for (j=0;j<C*N;j++)
2117          pcm[j] = FLOAT2INT16 (out[j]);
2118    
2119    RESTORE_STACK;
2120    return ret;
2121 }
2122 #endif
2123
2124 int celt_decode(CELTDecoder * restrict st, const unsigned char *data, int len, celt_int16 * restrict pcm, int frame_size)
2125 {
2126    return celt_decode_with_ec(st, data, len, pcm, frame_size, NULL);
2127 }
2128
2129 #ifndef DISABLE_FLOAT_API
2130 int celt_decode_float(CELTDecoder * restrict st, const unsigned char *data, int len, float * restrict pcm, int frame_size)
2131 {
2132    return celt_decode_with_ec_float(st, data, len, pcm, frame_size, NULL);
2133 }
2134 #endif /* DISABLE_FLOAT_API */
2135
2136 int celt_decoder_ctl(CELTDecoder * restrict st, int request, ...)
2137 {
2138    va_list ap;
2139
2140    va_start(ap, request);
2141    switch (request)
2142    {
2143       case CELT_GET_MODE_REQUEST:
2144       {
2145          const CELTMode ** value = va_arg(ap, const CELTMode**);
2146          if (value==0)
2147             goto bad_arg;
2148          *value=st->mode;
2149       }
2150       break;
2151       case CELT_SET_START_BAND_REQUEST:
2152       {
2153          celt_int32 value = va_arg(ap, celt_int32);
2154          if (value<0 || value>=st->mode->nbEBands)
2155             goto bad_arg;
2156          st->start = value;
2157       }
2158       break;
2159       case CELT_SET_END_BAND_REQUEST:
2160       {
2161          celt_int32 value = va_arg(ap, celt_int32);
2162          if (value<0 || value>=st->mode->nbEBands)
2163             goto bad_arg;
2164          st->end = value;
2165       }
2166       break;
2167       case CELT_RESET_STATE:
2168       {
2169          CELT_MEMSET((char*)&st->DECODER_RESET_START, 0,
2170                celt_decoder_get_size(st->mode, st->channels)-
2171                ((char*)&st->DECODER_RESET_START - (char*)st));
2172       }
2173       break;
2174       default:
2175          goto bad_request;
2176    }
2177    va_end(ap);
2178    return CELT_OK;
2179 bad_arg:
2180    va_end(ap);
2181    return CELT_BAD_ARG;
2182 bad_request:
2183       va_end(ap);
2184   return CELT_UNIMPLEMENTED;
2185 }
2186
2187 const char *celt_strerror(int error)
2188 {
2189    static const char *error_strings[8] = {
2190       "success",
2191       "invalid argument",
2192       "invalid mode",
2193       "internal error",
2194       "corrupted stream",
2195       "request not implemented",
2196       "invalid state",
2197       "memory allocation failed"
2198    };
2199    if (error > 0 || error < -7)
2200       return "unknown error";
2201    else 
2202       return error_strings[-error];
2203 }
2204