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