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