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