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