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