99ce9199608424a30ebcf5957d5aecb4ffc52dcb
[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    /* Clamp the band energy used for prediction */
1382    for (i=0;i<C*st->mode->nbEBands;i++)
1383       if (oldBandE[i] < -QCONST16(28.f,DB_SHIFT))
1384          oldBandE[i] = -QCONST16(28.f,DB_SHIFT);
1385
1386    /* In case start or end were to change */
1387    c=0; do
1388    {
1389       for (i=0;i<st->start;i++)
1390          oldBandE[c*st->mode->nbEBands+i]=0;
1391       for (i=st->end;i<st->mode->nbEBands;i++)
1392          oldBandE[c*st->mode->nbEBands+i]=0;
1393    } while (++c<C);
1394    for (i=0;i<C*st->mode->nbEBands;i++)
1395       oldLogE2[i] = oldLogE[i];
1396    if (isTransient)
1397       st->consec_transient++;
1398    else
1399       st->consec_transient=0;
1400    st->rng = enc->rng;
1401
1402    /* If there's any room left (can only happen for very high rates),
1403       fill it with zeros */
1404    while (ec_enc_tell(enc,0) + 8 <= nbCompressedBytes*8)
1405       ec_enc_bits(enc, 0, 8);
1406    ec_enc_done(enc);
1407    
1408    RESTORE_STACK;
1409    if (ec_enc_get_error(enc))
1410       return CELT_CORRUPTED_DATA;
1411    else
1412       return nbCompressedBytes;
1413 }
1414
1415 #ifdef FIXED_POINT
1416 #ifndef DISABLE_FLOAT_API
1417 int celt_encode_with_ec_float(CELTEncoder * restrict st, const float * pcm, int frame_size, unsigned char *compressed, int nbCompressedBytes, ec_enc *enc)
1418 {
1419    int j, ret, C, N, LM, M;
1420    VARDECL(celt_int16, in);
1421    SAVE_STACK;
1422
1423    if (pcm==NULL)
1424       return CELT_BAD_ARG;
1425
1426    for (LM=0;LM<4;LM++)
1427       if (st->mode->shortMdctSize<<LM==frame_size)
1428          break;
1429    if (LM>=MAX_CONFIG_SIZES)
1430       return CELT_BAD_ARG;
1431    M=1<<LM;
1432
1433    C = CHANNELS(st->channels);
1434    N = M*st->mode->shortMdctSize;
1435    ALLOC(in, C*N, celt_int16);
1436
1437    for (j=0;j<C*N;j++)
1438      in[j] = FLOAT2INT16(pcm[j]);
1439
1440    ret=celt_encode_with_ec(st,in,frame_size,compressed,nbCompressedBytes, enc);
1441 #ifdef RESYNTH
1442    for (j=0;j<C*N;j++)
1443       ((float*)pcm)[j]=in[j]*(1.f/32768.f);
1444 #endif
1445    RESTORE_STACK;
1446    return ret;
1447
1448 }
1449 #endif /*DISABLE_FLOAT_API*/
1450 #else
1451 int celt_encode_with_ec(CELTEncoder * restrict st, const celt_int16 * pcm, int frame_size, unsigned char *compressed, int nbCompressedBytes, ec_enc *enc)
1452 {
1453    int j, ret, C, N, LM, M;
1454    VARDECL(celt_sig, in);
1455    SAVE_STACK;
1456
1457    if (pcm==NULL)
1458       return CELT_BAD_ARG;
1459
1460    for (LM=0;LM<4;LM++)
1461       if (st->mode->shortMdctSize<<LM==frame_size)
1462          break;
1463    if (LM>=MAX_CONFIG_SIZES)
1464       return CELT_BAD_ARG;
1465    M=1<<LM;
1466
1467    C=CHANNELS(st->channels);
1468    N=M*st->mode->shortMdctSize;
1469    ALLOC(in, C*N, celt_sig);
1470    for (j=0;j<C*N;j++) {
1471      in[j] = SCALEOUT(pcm[j]);
1472    }
1473
1474    ret = celt_encode_with_ec_float(st,in,frame_size,compressed,nbCompressedBytes, enc);
1475 #ifdef RESYNTH
1476    for (j=0;j<C*N;j++)
1477       ((celt_int16*)pcm)[j] = FLOAT2INT16(in[j]);
1478 #endif
1479    RESTORE_STACK;
1480    return ret;
1481 }
1482 #endif
1483
1484 int celt_encode(CELTEncoder * restrict st, const celt_int16 * pcm, int frame_size, unsigned char *compressed, int nbCompressedBytes)
1485 {
1486    return celt_encode_with_ec(st, pcm, frame_size, compressed, nbCompressedBytes, NULL);
1487 }
1488
1489 #ifndef DISABLE_FLOAT_API
1490 int celt_encode_float(CELTEncoder * restrict st, const float * pcm, int frame_size, unsigned char *compressed, int nbCompressedBytes)
1491 {
1492    return celt_encode_with_ec_float(st, pcm, frame_size, compressed, nbCompressedBytes, NULL);
1493 }
1494 #endif /* DISABLE_FLOAT_API */
1495
1496 int celt_encoder_ctl(CELTEncoder * restrict st, int request, ...)
1497 {
1498    va_list ap;
1499    
1500    va_start(ap, request);
1501    switch (request)
1502    {
1503       case CELT_GET_MODE_REQUEST:
1504       {
1505          const CELTMode ** value = va_arg(ap, const CELTMode**);
1506          if (value==0)
1507             goto bad_arg;
1508          *value=st->mode;
1509       }
1510       break;
1511       case CELT_SET_COMPLEXITY_REQUEST:
1512       {
1513          int value = va_arg(ap, celt_int32);
1514          if (value<0 || value>10)
1515             goto bad_arg;
1516          st->complexity = value;
1517       }
1518       break;
1519       case CELT_SET_START_BAND_REQUEST:
1520       {
1521          celt_int32 value = va_arg(ap, celt_int32);
1522          if (value<0 || value>=st->mode->nbEBands)
1523             goto bad_arg;
1524          st->start = value;
1525       }
1526       break;
1527       case CELT_SET_END_BAND_REQUEST:
1528       {
1529          celt_int32 value = va_arg(ap, celt_int32);
1530          if (value<1 || value>st->mode->nbEBands)
1531             goto bad_arg;
1532          st->end = value;
1533       }
1534       break;
1535       case CELT_SET_PREDICTION_REQUEST:
1536       {
1537          int value = va_arg(ap, celt_int32);
1538          if (value<0 || value>2)
1539             goto bad_arg;
1540          if (value==0)
1541          {
1542             st->force_intra   = 1;
1543          } else if (value==1) {
1544             st->force_intra   = 0;
1545          } else {
1546             st->force_intra   = 0;
1547          }   
1548       }
1549       break;
1550       case CELT_SET_VBR_CONSTRAINT_REQUEST:
1551       {
1552          celt_int32 value = va_arg(ap, celt_int32);
1553          st->constrained_vbr = value;
1554       }
1555       break;
1556       case CELT_SET_VBR_RATE_REQUEST:
1557       {
1558          celt_int32 value = va_arg(ap, celt_int32);
1559          int frame_rate;
1560          int N = st->mode->shortMdctSize;
1561          if (value<0)
1562             goto bad_arg;
1563          if (value>3072000)
1564             value = 3072000;
1565          frame_rate = ((st->mode->Fs<<3)+(N>>1))/N;
1566          st->vbr_rate_norm = ((value<<(BITRES+3))+(frame_rate>>1))/frame_rate;
1567       }
1568       break;
1569       case CELT_RESET_STATE:
1570       {
1571          CELT_MEMSET((char*)&st->ENCODER_RESET_START, 0,
1572                celt_encoder_get_size(st->mode, st->channels)-
1573                ((char*)&st->ENCODER_RESET_START - (char*)st));
1574          st->vbr_offset = 0;
1575          st->delayedIntra = 1;
1576          st->spread_decision = SPREAD_NORMAL;
1577          st->tonal_average = QCONST16(1.f,8);
1578       }
1579       break;
1580       default:
1581          goto bad_request;
1582    }
1583    va_end(ap);
1584    return CELT_OK;
1585 bad_arg:
1586    va_end(ap);
1587    return CELT_BAD_ARG;
1588 bad_request:
1589    va_end(ap);
1590    return CELT_UNIMPLEMENTED;
1591 }
1592
1593 /**********************************************************************/
1594 /*                                                                    */
1595 /*                             DECODER                                */
1596 /*                                                                    */
1597 /**********************************************************************/
1598 #define DECODE_BUFFER_SIZE 2048
1599
1600 /** Decoder state 
1601  @brief Decoder state
1602  */
1603 struct CELTDecoder {
1604    const CELTMode *mode;
1605    int overlap;
1606    int channels;
1607
1608    int start, end;
1609
1610    /* Everything beyond this point gets cleared on a reset */
1611 #define DECODER_RESET_START rng
1612
1613    ec_uint32 rng;
1614    int last_pitch_index;
1615    int loss_count;
1616    int postfilter_period;
1617    int postfilter_period_old;
1618    celt_word16 postfilter_gain;
1619    celt_word16 postfilter_gain_old;
1620    int postfilter_tapset;
1621    int postfilter_tapset_old;
1622
1623    celt_sig preemph_memD[2];
1624    
1625    celt_sig _decode_mem[1]; /* Size = channels*(DECODE_BUFFER_SIZE+mode->overlap) */
1626    /* celt_word16 lpc[],  Size = channels*LPC_ORDER */
1627    /* celt_word16 oldEBands[], Size = channels*mode->nbEBands */
1628    /* celt_word16 oldLogE2[], Size = channels*mode->nbEBands */
1629    /* celt_word16 backgroundLogE[], Size = channels*mode->nbEBands */
1630 };
1631
1632 int celt_decoder_get_size(const CELTMode *mode, int channels)
1633 {
1634    int size = sizeof(struct CELTDecoder)
1635             + (channels*(DECODE_BUFFER_SIZE+mode->overlap)-1)*sizeof(celt_sig)
1636             + channels*LPC_ORDER*sizeof(celt_word16)
1637             + 3*channels*mode->nbEBands*sizeof(celt_word16);
1638    return size;
1639 }
1640
1641 CELTDecoder *celt_decoder_create(const CELTMode *mode, int channels, int *error)
1642 {
1643    return celt_decoder_init(
1644          (CELTDecoder *)celt_alloc(celt_decoder_get_size(mode, channels)),
1645          mode, channels, error);
1646 }
1647
1648 CELTDecoder *celt_decoder_init(CELTDecoder *st, const CELTMode *mode, int channels, int *error)
1649 {
1650    if (channels < 0 || channels > 2)
1651    {
1652       if (error)
1653          *error = CELT_BAD_ARG;
1654       return NULL;
1655    }
1656
1657    if (st==NULL)
1658    {
1659       if (error)
1660          *error = CELT_ALLOC_FAIL;
1661       return NULL;
1662    }
1663
1664    CELT_MEMSET((char*)st, 0, celt_decoder_get_size(mode, channels));
1665
1666    st->mode = mode;
1667    st->overlap = mode->overlap;
1668    st->channels = channels;
1669
1670    st->start = 0;
1671    st->end = st->mode->effEBands;
1672
1673    st->loss_count = 0;
1674
1675    if (error)
1676       *error = CELT_OK;
1677    return st;
1678 }
1679
1680 void celt_decoder_destroy(CELTDecoder *st)
1681 {
1682    celt_free(st);
1683 }
1684
1685 static void celt_decode_lost(CELTDecoder * restrict st, celt_word16 * restrict pcm, int N, int LM)
1686 {
1687    int c;
1688    int pitch_index;
1689    int overlap = st->mode->overlap;
1690    celt_word16 fade = Q15ONE;
1691    int i, len;
1692    const int C = CHANNELS(st->channels);
1693    int offset;
1694    celt_sig *out_mem[2];
1695    celt_sig *decode_mem[2];
1696    celt_sig *overlap_mem[2];
1697    celt_word16 *lpc;
1698    SAVE_STACK;
1699    
1700    c=0; do {
1701       decode_mem[c] = st->_decode_mem + c*(DECODE_BUFFER_SIZE+st->overlap);
1702       out_mem[c] = decode_mem[c]+DECODE_BUFFER_SIZE-MAX_PERIOD;
1703       overlap_mem[c] = decode_mem[c]+DECODE_BUFFER_SIZE;
1704    } while (++c<C);
1705    lpc = (celt_word16*)(st->_decode_mem+(DECODE_BUFFER_SIZE+st->overlap)*C);
1706
1707    len = N+st->mode->overlap;
1708    
1709    if (st->loss_count == 0)
1710    {
1711       celt_word16 pitch_buf[MAX_PERIOD>>1];
1712       celt_word32 tmp=0;
1713       celt_word32 mem0[2]={0,0};
1714       celt_word16 mem1[2]={0,0};
1715       int len2 = len;
1716       /* FIXME: This is a kludge */
1717       if (len2>MAX_PERIOD>>1)
1718          len2 = MAX_PERIOD>>1;
1719       pitch_downsample(out_mem, pitch_buf, MAX_PERIOD, MAX_PERIOD,
1720                        C, mem0, mem1);
1721       pitch_search(st->mode, pitch_buf+((MAX_PERIOD-len2)>>1), pitch_buf, len2,
1722                    MAX_PERIOD-len2-100, &pitch_index, &tmp, 1<<LM);
1723       pitch_index = MAX_PERIOD-len2-pitch_index;
1724       st->last_pitch_index = pitch_index;
1725    } else {
1726       pitch_index = st->last_pitch_index;
1727       if (st->loss_count < 5)
1728          fade = QCONST16(.8f,15);
1729       else
1730          fade = 0;
1731    }
1732
1733    c=0; do {
1734       /* FIXME: This is more memory than necessary */
1735       celt_word32 e[2*MAX_PERIOD];
1736       celt_word16 exc[2*MAX_PERIOD];
1737       celt_word32 ac[LPC_ORDER+1];
1738       celt_word16 decay = 1;
1739       celt_word32 S1=0;
1740       celt_word16 mem[LPC_ORDER]={0};
1741
1742       offset = MAX_PERIOD-pitch_index;
1743       for (i=0;i<MAX_PERIOD;i++)
1744          exc[i] = ROUND16(out_mem[c][i], SIG_SHIFT);
1745
1746       if (st->loss_count == 0)
1747       {
1748          _celt_autocorr(exc, ac, st->mode->window, st->mode->overlap,
1749                         LPC_ORDER, MAX_PERIOD);
1750
1751          /* Noise floor -40 dB */
1752 #ifdef FIXED_POINT
1753          ac[0] += SHR32(ac[0],13);
1754 #else
1755          ac[0] *= 1.0001f;
1756 #endif
1757          /* Lag windowing */
1758          for (i=1;i<=LPC_ORDER;i++)
1759          {
1760             /*ac[i] *= exp(-.5*(2*M_PI*.002*i)*(2*M_PI*.002*i));*/
1761 #ifdef FIXED_POINT
1762             ac[i] -= MULT16_32_Q15(2*i*i, ac[i]);
1763 #else
1764             ac[i] -= ac[i]*(.008f*i)*(.008f*i);
1765 #endif
1766          }
1767
1768          _celt_lpc(lpc+c*LPC_ORDER, ac, LPC_ORDER);
1769       }
1770       for (i=0;i<LPC_ORDER;i++)
1771          mem[i] = ROUND16(out_mem[c][MAX_PERIOD-1-i], SIG_SHIFT);
1772       fir(exc, lpc+c*LPC_ORDER, exc, MAX_PERIOD, LPC_ORDER, mem);
1773       /*for (i=0;i<MAX_PERIOD;i++)printf("%d ", exc[i]); printf("\n");*/
1774       /* Check if the waveform is decaying (and if so how fast) */
1775       {
1776          celt_word32 E1=1, E2=1;
1777          int period;
1778          if (pitch_index <= MAX_PERIOD/2)
1779             period = pitch_index;
1780          else
1781             period = MAX_PERIOD/2;
1782          for (i=0;i<period;i++)
1783          {
1784             E1 += SHR32(MULT16_16(exc[MAX_PERIOD-period+i],exc[MAX_PERIOD-period+i]),8);
1785             E2 += SHR32(MULT16_16(exc[MAX_PERIOD-2*period+i],exc[MAX_PERIOD-2*period+i]),8);
1786          }
1787          if (E1 > E2)
1788             E1 = E2;
1789          decay = celt_sqrt(frac_div32(SHR(E1,1),E2));
1790       }
1791
1792       /* Copy excitation, taking decay into account */
1793       for (i=0;i<len+st->mode->overlap;i++)
1794       {
1795          celt_word16 tmp;
1796          if (offset+i >= MAX_PERIOD)
1797          {
1798             offset -= pitch_index;
1799             decay = MULT16_16_Q15(decay, decay);
1800          }
1801          e[i] = SHL32(EXTEND32(MULT16_16_Q15(decay, exc[offset+i])), SIG_SHIFT);
1802          tmp = ROUND16(out_mem[c][offset+i],SIG_SHIFT);
1803          S1 += SHR32(MULT16_16(tmp,tmp),8);
1804       }
1805       for (i=0;i<LPC_ORDER;i++)
1806          mem[i] = ROUND16(out_mem[c][MAX_PERIOD-1-i], SIG_SHIFT);
1807       for (i=0;i<len+st->mode->overlap;i++)
1808          e[i] = MULT16_32_Q15(fade, e[i]);
1809       iir(e, lpc+c*LPC_ORDER, e, len+st->mode->overlap, LPC_ORDER, mem);
1810
1811       {
1812          celt_word32 S2=0;
1813          for (i=0;i<len+overlap;i++)
1814          {
1815             celt_word16 tmp = ROUND16(e[i],SIG_SHIFT);
1816             S2 += SHR32(MULT16_16(tmp,tmp),8);
1817          }
1818          /* This checks for an "explosion" in the synthesis */
1819 #ifdef FIXED_POINT
1820          if (!(S1 > SHR32(S2,2)))
1821 #else
1822          /* Float test is written this way to catch NaNs at the same time */
1823          if (!(S1 > 0.2f*S2))
1824 #endif
1825          {
1826             for (i=0;i<len+overlap;i++)
1827                e[i] = 0;
1828          } else if (S1 < S2)
1829          {
1830             celt_word16 ratio = celt_sqrt(frac_div32(SHR32(S1,1)+1,S2+1));
1831             for (i=0;i<len+overlap;i++)
1832                e[i] = MULT16_32_Q15(ratio, e[i]);
1833          }
1834       }
1835
1836 #ifdef ENABLE_POSTFILTER
1837       /* Apply post-filter to the MDCT overlap of the previous frame */
1838       comb_filter(out_mem[c]+MAX_PERIOD, out_mem[c]+MAX_PERIOD, st->postfilter_period, st->postfilter_period, st->overlap, C,
1839                   st->postfilter_gain, st->postfilter_gain, st->postfilter_tapset, st->postfilter_tapset,
1840                   NULL, 0);
1841 #endif /* ENABLE_POSTFILTER */
1842
1843       for (i=0;i<MAX_PERIOD+st->mode->overlap-N;i++)
1844          out_mem[c][i] = out_mem[c][N+i];
1845
1846       /* Apply TDAC to the concealed audio so that it blends with the
1847          previous and next frames */
1848       for (i=0;i<overlap/2;i++)
1849       {
1850          celt_word32 tmp;
1851          tmp = MULT16_32_Q15(st->mode->window[i],           e[N+overlap-1-i]) +
1852                MULT16_32_Q15(st->mode->window[overlap-i-1], e[N+i          ]);
1853          out_mem[c][MAX_PERIOD+i] = MULT16_32_Q15(st->mode->window[overlap-i-1], tmp);
1854          out_mem[c][MAX_PERIOD+overlap-i-1] = MULT16_32_Q15(st->mode->window[i], tmp);
1855       }
1856       for (i=0;i<N;i++)
1857          out_mem[c][MAX_PERIOD-N+i] = e[i];
1858
1859 #ifdef ENABLE_POSTFILTER
1860       /* Apply pre-filter to the MDCT overlap for the next frame (post-filter will be applied then) */
1861       comb_filter(e, out_mem[c]+MAX_PERIOD, st->postfilter_period, st->postfilter_period, st->overlap, C,
1862                   -st->postfilter_gain, -st->postfilter_gain, st->postfilter_tapset, st->postfilter_tapset,
1863                   NULL, 0);
1864 #endif /* ENABLE_POSTFILTER */
1865       for (i=0;i<overlap;i++)
1866          out_mem[c][MAX_PERIOD+i] = e[i];
1867    } while (++c<C);
1868
1869    {
1870       celt_word32 *out_syn[2];
1871       out_syn[0] = out_mem[0]+MAX_PERIOD-N;
1872       if (C==2)
1873          out_syn[1] = out_mem[1]+MAX_PERIOD-N;
1874       deemphasis(out_syn, pcm, N, C, st->mode->preemph, st->preemph_memD);
1875    }
1876    
1877    st->loss_count++;
1878
1879    RESTORE_STACK;
1880 }
1881
1882 #ifdef FIXED_POINT
1883 int celt_decode_with_ec(CELTDecoder * restrict st, const unsigned char *data, int len, celt_int16 * restrict pcm, int frame_size, ec_dec *dec)
1884 {
1885 #else
1886 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)
1887 {
1888 #endif
1889    int c, i, N;
1890    int spread_decision;
1891    int bits;
1892    ec_dec _dec;
1893    ec_byte_buffer buf;
1894    VARDECL(celt_sig, freq);
1895    VARDECL(celt_norm, X);
1896    VARDECL(celt_ener, bandE);
1897    VARDECL(celt_word16, oldLogE);
1898    VARDECL(int, fine_quant);
1899    VARDECL(int, pulses);
1900    VARDECL(int, offsets);
1901    VARDECL(int, fine_priority);
1902    VARDECL(int, tf_res);
1903    VARDECL(unsigned char, collapse_masks);
1904    celt_sig *out_mem[2];
1905    celt_sig *decode_mem[2];
1906    celt_sig *overlap_mem[2];
1907    celt_sig *out_syn[2];
1908    celt_word16 *lpc;
1909    celt_word16 *oldBandE, *oldLogE2, *backgroundLogE;
1910
1911    int shortBlocks;
1912    int isTransient;
1913    int intra_ener;
1914    const int C = CHANNELS(st->channels);
1915    int LM, M;
1916    int effEnd;
1917    int codedBands;
1918    int alloc_trim;
1919    int postfilter_pitch;
1920    celt_word16 postfilter_gain;
1921    int intensity=0;
1922    int dual_stereo=0;
1923    celt_int32 total_bits;
1924    celt_int32 tell;
1925    int dynalloc_logp;
1926    int postfilter_tapset;
1927    int anti_collapse_rsv;
1928    int anti_collapse_on=0;
1929    SAVE_STACK;
1930
1931    if (pcm==NULL)
1932       return CELT_BAD_ARG;
1933
1934    for (LM=0;LM<4;LM++)
1935       if (st->mode->shortMdctSize<<LM==frame_size)
1936          break;
1937    if (LM>=MAX_CONFIG_SIZES)
1938       return CELT_BAD_ARG;
1939    M=1<<LM;
1940
1941    c=0; do {
1942       decode_mem[c] = st->_decode_mem + c*(DECODE_BUFFER_SIZE+st->overlap);
1943       out_mem[c] = decode_mem[c]+DECODE_BUFFER_SIZE-MAX_PERIOD;
1944       overlap_mem[c] = decode_mem[c]+DECODE_BUFFER_SIZE;
1945    } while (++c<C);
1946    lpc = (celt_word16*)(st->_decode_mem+(DECODE_BUFFER_SIZE+st->overlap)*C);
1947    oldBandE = lpc+C*LPC_ORDER;
1948    oldLogE2 = oldBandE + C*st->mode->nbEBands;
1949    backgroundLogE = oldLogE2  + C*st->mode->nbEBands;
1950
1951    N = M*st->mode->shortMdctSize;
1952
1953    effEnd = st->end;
1954    if (effEnd > st->mode->effEBands)
1955       effEnd = st->mode->effEBands;
1956
1957    ALLOC(freq, C*N, celt_sig); /**< Interleaved signal MDCTs */
1958    ALLOC(X, C*N, celt_norm);   /**< Interleaved normalised MDCTs */
1959    ALLOC(bandE, st->mode->nbEBands*C, celt_ener);
1960    c=0; do
1961       for (i=0;i<M*st->mode->eBands[st->start];i++)
1962          X[c*N+i] = 0;
1963    while (++c<C);
1964    c=0; do   
1965       for (i=M*st->mode->eBands[effEnd];i<N;i++)
1966          X[c*N+i] = 0;
1967    while (++c<C);
1968
1969    if (data == NULL)
1970    {
1971       celt_decode_lost(st, pcm, N, LM);
1972       RESTORE_STACK;
1973       return CELT_OK;
1974    }
1975    if (len<0) {
1976      RESTORE_STACK;
1977      return CELT_BAD_ARG;
1978    }
1979    
1980    if (dec == NULL)
1981    {
1982       ec_byte_readinit(&buf,(unsigned char*)data,len);
1983       ec_dec_init(&_dec,&buf);
1984       dec = &_dec;
1985    }
1986
1987    total_bits = len*8;
1988    tell = ec_dec_tell(dec, 0);
1989
1990    postfilter_gain = 0;
1991    postfilter_pitch = 0;
1992    postfilter_tapset = 0;
1993    if (tell+17 <= total_bits)
1994    {
1995       if(ec_dec_bit_logp(dec, 1))
1996       {
1997 #ifdef ENABLE_POSTFILTER
1998          int qg, octave;
1999          octave = ec_dec_uint(dec, 6);
2000          postfilter_pitch = (16<<octave)+ec_dec_bits(dec, 4+octave)-1;
2001          qg = ec_dec_bits(dec, 2);
2002          postfilter_tapset = ec_dec_icdf(dec, tapset_icdf, 2);
2003          postfilter_gain = QCONST16(.125f,15)*(qg+2);
2004 #else /* ENABLE_POSTFILTER */
2005          RESTORE_STACK;
2006          return CELT_CORRUPTED_DATA;
2007 #endif /* ENABLE_POSTFILTER */
2008       }
2009       tell = ec_dec_tell(dec, 0);
2010    }
2011
2012    if (LM > 0 && tell+3 <= total_bits)
2013    {
2014       isTransient = ec_dec_bit_logp(dec, 3);
2015       tell = ec_dec_tell(dec, 0);
2016    }
2017    else
2018       isTransient = 0;
2019
2020    if (isTransient)
2021       shortBlocks = M;
2022    else
2023       shortBlocks = 0;
2024
2025    ALLOC(oldLogE, C*st->mode->nbEBands, celt_word16);
2026    for (i=0;i<C*st->mode->nbEBands;i++)
2027       oldLogE[i] = oldBandE[i];
2028
2029    /* Decode the global flags (first symbols in the stream) */
2030    intra_ener = tell+3<=total_bits ? ec_dec_bit_logp(dec, 3) : 0;
2031    /* Get band energies */
2032    unquant_coarse_energy(st->mode, st->start, st->end, bandE, oldBandE,
2033          intra_ener, dec, C, LM);
2034
2035    ALLOC(tf_res, st->mode->nbEBands, int);
2036    tf_decode(st->start, st->end, C, isTransient, tf_res, LM, dec);
2037
2038    tell = ec_dec_tell(dec, 0);
2039    spread_decision = SPREAD_NORMAL;
2040    if (tell+4 <= total_bits)
2041       spread_decision = ec_dec_icdf(dec, spread_icdf, 5);
2042
2043    ALLOC(pulses, st->mode->nbEBands, int);
2044    ALLOC(offsets, st->mode->nbEBands, int);
2045    ALLOC(fine_priority, st->mode->nbEBands, int);
2046
2047    dynalloc_logp = 6;
2048    total_bits<<=BITRES;
2049    tell = ec_dec_tell(dec, BITRES);
2050    for (i=st->start;i<st->end;i++)
2051    {
2052       int width, quanta;
2053       int dynalloc_loop_logp;
2054       int boost;
2055       width = C*(st->mode->eBands[i+1]-st->mode->eBands[i])<<LM;
2056       /* quanta is 6 bits, but no more than 1 bit/sample
2057          and no less than 1/8 bit/sample */
2058       quanta = IMIN(width<<BITRES, IMAX(6<<BITRES, width));
2059       dynalloc_loop_logp = dynalloc_logp;
2060       boost = 0;
2061       while (tell+(dynalloc_loop_logp<<BITRES) < total_bits &&
2062             boost < (64<<LM)*(C<<BITRES))
2063       {
2064          int flag;
2065          flag = ec_dec_bit_logp(dec, dynalloc_loop_logp);
2066          tell = ec_dec_tell(dec, BITRES);
2067          if (!flag)
2068             break;
2069          boost += quanta;
2070          total_bits -= quanta;
2071          dynalloc_loop_logp = 1;
2072       }
2073       offsets[i] = boost;
2074       /* Making dynalloc more likely */
2075       if (boost>0)
2076          dynalloc_logp = IMAX(2, dynalloc_logp-1);
2077    }
2078
2079    ALLOC(fine_quant, st->mode->nbEBands, int);
2080    alloc_trim = tell+(6<<BITRES) <= total_bits ?
2081          ec_dec_icdf(dec, trim_icdf, 7) : 5;
2082
2083    bits = (len*8<<BITRES) - ec_dec_tell(dec, BITRES) - 1;
2084    anti_collapse_rsv = isTransient&&LM>=2&&bits>=(LM+2<<BITRES) ? (1<<BITRES) : 0;
2085    bits -= anti_collapse_rsv;
2086    codedBands = compute_allocation(st->mode, st->start, st->end, offsets,
2087          alloc_trim, &intensity, &dual_stereo, bits, pulses, fine_quant,
2088          fine_priority, C, LM, dec, 0, 0);
2089    
2090    unquant_fine_energy(st->mode, st->start, st->end, bandE, oldBandE, fine_quant, dec, C);
2091
2092    /* Decode fixed codebook */
2093    ALLOC(collapse_masks, st->mode->nbEBands, unsigned char);
2094    quant_all_bands(0, st->mode, st->start, st->end, X, C==2 ? X+N : NULL, collapse_masks,
2095          NULL, pulses, shortBlocks, spread_decision, dual_stereo, intensity, tf_res, 1,
2096          len*8, dec, LM, codedBands, &st->rng);
2097
2098    if (anti_collapse_rsv > 0)
2099    {
2100       anti_collapse_on = ec_dec_bits(dec, 1);
2101    }
2102
2103    unquant_energy_finalise(st->mode, st->start, st->end, bandE, oldBandE,
2104          fine_quant, fine_priority, len*8-ec_dec_tell(dec, 0), dec, C);
2105
2106    if (anti_collapse_on)
2107       anti_collapse(st->mode, X, collapse_masks, LM, C, N,
2108             st->start, st->end, oldBandE, oldLogE, oldLogE2, pulses, st->rng);
2109
2110    log2Amp(st->mode, st->start, st->end, bandE, oldBandE, C);
2111
2112    /* Synthesis */
2113    denormalise_bands(st->mode, X, freq, bandE, effEnd, C, M);
2114
2115    CELT_MOVE(decode_mem[0], decode_mem[0]+N, DECODE_BUFFER_SIZE-N);
2116    if (C==2)
2117       CELT_MOVE(decode_mem[1], decode_mem[1]+N, DECODE_BUFFER_SIZE-N);
2118
2119    c=0; do
2120       for (i=0;i<M*st->mode->eBands[st->start];i++)
2121          freq[c*N+i] = 0;
2122    while (++c<C);
2123    c=0; do
2124       for (i=M*st->mode->eBands[effEnd];i<N;i++)
2125          freq[c*N+i] = 0;
2126    while (++c<C);
2127
2128    out_syn[0] = out_mem[0]+MAX_PERIOD-N;
2129    if (C==2)
2130       out_syn[1] = out_mem[1]+MAX_PERIOD-N;
2131
2132    /* Compute inverse MDCTs */
2133    compute_inv_mdcts(st->mode, shortBlocks, freq, out_syn, overlap_mem, C, LM);
2134
2135 #ifdef ENABLE_POSTFILTER
2136    c=0; do {
2137       st->postfilter_period=IMAX(st->postfilter_period, COMBFILTER_MINPERIOD);
2138       st->postfilter_period_old=IMAX(st->postfilter_period_old, COMBFILTER_MINPERIOD);
2139       if (LM!=0)
2140       {
2141          comb_filter(out_syn[c], out_syn[c], st->postfilter_period, st->postfilter_period, st->overlap, C,
2142                st->postfilter_gain, st->postfilter_gain, st->postfilter_tapset, st->postfilter_tapset,
2143                NULL, 0);
2144          comb_filter(out_syn[c]+st->overlap, out_syn[c]+st->overlap, st->postfilter_period, postfilter_pitch, N-st->overlap, C,
2145                st->postfilter_gain, postfilter_gain, st->postfilter_tapset, postfilter_tapset,
2146                st->mode->window, st->mode->overlap);
2147       } else {
2148          comb_filter(out_syn[c], out_syn[c], st->postfilter_period_old, st->postfilter_period, N-st->overlap, C,
2149                st->postfilter_gain_old, st->postfilter_gain, st->postfilter_tapset_old, st->postfilter_tapset,
2150                st->mode->window, st->mode->overlap);
2151       }
2152    } while (++c<C);
2153    st->postfilter_period_old = st->postfilter_period;
2154    st->postfilter_gain_old = st->postfilter_gain;
2155    st->postfilter_tapset_old = st->postfilter_tapset;
2156    st->postfilter_period = postfilter_pitch;
2157    st->postfilter_gain = postfilter_gain;
2158    st->postfilter_tapset = postfilter_tapset;
2159 #endif /* ENABLE_POSTFILTER */
2160
2161    /* Clamp the band energy used for prediction */
2162    for (i=0;i<C*st->mode->nbEBands;i++)
2163       if (oldBandE[i] < -QCONST16(28.f,DB_SHIFT))
2164          oldBandE[i] = -QCONST16(28.f,DB_SHIFT);
2165
2166    /* In case start or end were to change */
2167    c=0; do
2168    {
2169       for (i=0;i<st->start;i++)
2170          oldBandE[c*st->mode->nbEBands+i]=0;
2171       for (i=st->end;i<st->mode->nbEBands;i++)
2172          oldBandE[c*st->mode->nbEBands+i]=0;
2173    } while (++c<C);
2174    for (i=0;i<C*st->mode->nbEBands;i++)
2175       oldLogE2[i] = oldLogE[i];
2176    for (i=0;i<C*st->mode->nbEBands;i++)
2177       backgroundLogE[i] = MIN16(backgroundLogE[i] + M*QCONST16(0.001f,DB_SHIFT), oldLogE[i]);
2178
2179    st->rng = dec->rng;
2180
2181    deemphasis(out_syn, pcm, N, C, st->mode->preemph, st->preemph_memD);
2182    st->loss_count = 0;
2183    RESTORE_STACK;
2184    if (ec_dec_tell(dec,0) > 8*len || ec_dec_get_error(dec))
2185       return CELT_CORRUPTED_DATA;
2186    else
2187       return CELT_OK;
2188 }
2189
2190 #ifdef FIXED_POINT
2191 #ifndef DISABLE_FLOAT_API
2192 int celt_decode_with_ec_float(CELTDecoder * restrict st, const unsigned char *data, int len, float * restrict pcm, int frame_size, ec_dec *dec)
2193 {
2194    int j, ret, C, N, LM, M;
2195    VARDECL(celt_int16, out);
2196    SAVE_STACK;
2197
2198    if (pcm==NULL)
2199       return CELT_BAD_ARG;
2200
2201    for (LM=0;LM<4;LM++)
2202       if (st->mode->shortMdctSize<<LM==frame_size)
2203          break;
2204    if (LM>=MAX_CONFIG_SIZES)
2205       return CELT_BAD_ARG;
2206    M=1<<LM;
2207
2208    C = CHANNELS(st->channels);
2209    N = M*st->mode->shortMdctSize;
2210    
2211    ALLOC(out, C*N, celt_int16);
2212    ret=celt_decode_with_ec(st, data, len, out, frame_size, dec);
2213    if (ret==0)
2214       for (j=0;j<C*N;j++)
2215          pcm[j]=out[j]*(1.f/32768.f);
2216      
2217    RESTORE_STACK;
2218    return ret;
2219 }
2220 #endif /*DISABLE_FLOAT_API*/
2221 #else
2222 int celt_decode_with_ec(CELTDecoder * restrict st, const unsigned char *data, int len, celt_int16 * restrict pcm, int frame_size, ec_dec *dec)
2223 {
2224    int j, ret, C, N, LM, M;
2225    VARDECL(celt_sig, out);
2226    SAVE_STACK;
2227
2228    if (pcm==NULL)
2229       return CELT_BAD_ARG;
2230
2231    for (LM=0;LM<4;LM++)
2232       if (st->mode->shortMdctSize<<LM==frame_size)
2233          break;
2234    if (LM>=MAX_CONFIG_SIZES)
2235       return CELT_BAD_ARG;
2236    M=1<<LM;
2237
2238    C = CHANNELS(st->channels);
2239    N = M*st->mode->shortMdctSize;
2240    ALLOC(out, C*N, celt_sig);
2241
2242    ret=celt_decode_with_ec_float(st, data, len, out, frame_size, dec);
2243
2244    if (ret==0)
2245       for (j=0;j<C*N;j++)
2246          pcm[j] = FLOAT2INT16 (out[j]);
2247    
2248    RESTORE_STACK;
2249    return ret;
2250 }
2251 #endif
2252
2253 int celt_decode(CELTDecoder * restrict st, const unsigned char *data, int len, celt_int16 * restrict pcm, int frame_size)
2254 {
2255    return celt_decode_with_ec(st, data, len, pcm, frame_size, NULL);
2256 }
2257
2258 #ifndef DISABLE_FLOAT_API
2259 int celt_decode_float(CELTDecoder * restrict st, const unsigned char *data, int len, float * restrict pcm, int frame_size)
2260 {
2261    return celt_decode_with_ec_float(st, data, len, pcm, frame_size, NULL);
2262 }
2263 #endif /* DISABLE_FLOAT_API */
2264
2265 int celt_decoder_ctl(CELTDecoder * restrict st, int request, ...)
2266 {
2267    va_list ap;
2268
2269    va_start(ap, request);
2270    switch (request)
2271    {
2272       case CELT_GET_MODE_REQUEST:
2273       {
2274          const CELTMode ** value = va_arg(ap, const CELTMode**);
2275          if (value==0)
2276             goto bad_arg;
2277          *value=st->mode;
2278       }
2279       break;
2280       case CELT_SET_START_BAND_REQUEST:
2281       {
2282          celt_int32 value = va_arg(ap, celt_int32);
2283          if (value<0 || value>=st->mode->nbEBands)
2284             goto bad_arg;
2285          st->start = value;
2286       }
2287       break;
2288       case CELT_SET_END_BAND_REQUEST:
2289       {
2290          celt_int32 value = va_arg(ap, celt_int32);
2291          if (value<0 || value>=st->mode->nbEBands)
2292             goto bad_arg;
2293          st->end = value;
2294       }
2295       break;
2296       case CELT_RESET_STATE:
2297       {
2298          CELT_MEMSET((char*)&st->DECODER_RESET_START, 0,
2299                celt_decoder_get_size(st->mode, st->channels)-
2300                ((char*)&st->DECODER_RESET_START - (char*)st));
2301       }
2302       break;
2303       default:
2304          goto bad_request;
2305    }
2306    va_end(ap);
2307    return CELT_OK;
2308 bad_arg:
2309    va_end(ap);
2310    return CELT_BAD_ARG;
2311 bad_request:
2312       va_end(ap);
2313   return CELT_UNIMPLEMENTED;
2314 }
2315
2316 const char *celt_strerror(int error)
2317 {
2318    static const char *error_strings[8] = {
2319       "success",
2320       "invalid argument",
2321       "invalid mode",
2322       "internal error",
2323       "corrupted stream",
2324       "request not implemented",
2325       "invalid state",
2326       "memory allocation failed"
2327    };
2328    if (error > 0 || error < -7)
2329       return "unknown error";
2330    else 
2331       return error_strings[-error];
2332 }
2333