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