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