Fix effectiveBytes computation for VBR
[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    ec_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       float max_abs=0;
321       for (j=0;j<block;j++)
322          max_abs = MAX32(max_abs, 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    ec_uint32 budget;
691    ec_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    ec_uint32 budget;
731    ec_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    ec_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 (i=0;i<C*N;i++)
1971       {
1972             seed = lcg_rand(seed);
1973             X[i] = (celt_int32)(seed)>>20;
1974       }
1975       st->rng = seed;
1976       for (c=0;c<C;c++)
1977          for (i=0;i<st->mode->nbEBands;i++)
1978             renormalise_vector(X+N*c+(st->mode->eBands[i]<<LM), (st->mode->eBands[i+1]-st->mode->eBands[i])<<LM, Q15ONE);
1979
1980       denormalise_bands(st->mode, X, freq, bandE, st->mode->nbEBands, C, 1<<LM);
1981
1982       compute_inv_mdcts(st->mode, 0, freq, out_syn, overlap_mem, C, LM);
1983       plc = 0;
1984    } else if (st->loss_count == 0)
1985    {
1986       celt_word16 pitch_buf[MAX_PERIOD>>1];
1987       int len2 = len;
1988       /* FIXME: This is a kludge */
1989       if (len2>MAX_PERIOD>>1)
1990          len2 = MAX_PERIOD>>1;
1991       pitch_downsample(out_mem, pitch_buf, MAX_PERIOD, C);
1992       pitch_search(pitch_buf+((MAX_PERIOD-len2)>>1), pitch_buf, len2,
1993                    MAX_PERIOD-len2-100, &pitch_index);
1994       pitch_index = MAX_PERIOD-len2-pitch_index;
1995       st->last_pitch_index = pitch_index;
1996    } else {
1997       pitch_index = st->last_pitch_index;
1998       fade = QCONST16(.8f,15);
1999    }
2000
2001    if (plc)
2002    {
2003       c=0; do {
2004          /* FIXME: This is more memory than necessary */
2005          celt_word32 e[2*MAX_PERIOD];
2006          celt_word16 exc[2*MAX_PERIOD];
2007          celt_word32 ac[LPC_ORDER+1];
2008          celt_word16 decay = 1;
2009          celt_word32 S1=0;
2010          celt_word16 mem[LPC_ORDER]={0};
2011
2012          offset = MAX_PERIOD-pitch_index;
2013          for (i=0;i<MAX_PERIOD;i++)
2014             exc[i] = ROUND16(out_mem[c][i], SIG_SHIFT);
2015
2016          if (st->loss_count == 0)
2017          {
2018             _celt_autocorr(exc, ac, st->mode->window, st->mode->overlap,
2019                   LPC_ORDER, MAX_PERIOD);
2020
2021             /* Noise floor -40 dB */
2022 #ifdef FIXED_POINT
2023             ac[0] += SHR32(ac[0],13);
2024 #else
2025             ac[0] *= 1.0001f;
2026 #endif
2027             /* Lag windowing */
2028             for (i=1;i<=LPC_ORDER;i++)
2029             {
2030                /*ac[i] *= exp(-.5*(2*M_PI*.002*i)*(2*M_PI*.002*i));*/
2031 #ifdef FIXED_POINT
2032                ac[i] -= MULT16_32_Q15(2*i*i, ac[i]);
2033 #else
2034                ac[i] -= ac[i]*(.008f*i)*(.008f*i);
2035 #endif
2036             }
2037
2038             _celt_lpc(lpc+c*LPC_ORDER, ac, LPC_ORDER);
2039          }
2040          for (i=0;i<LPC_ORDER;i++)
2041             mem[i] = ROUND16(out_mem[c][MAX_PERIOD-1-i], SIG_SHIFT);
2042          fir(exc, lpc+c*LPC_ORDER, exc, MAX_PERIOD, LPC_ORDER, mem);
2043          /*for (i=0;i<MAX_PERIOD;i++)printf("%d ", exc[i]); printf("\n");*/
2044          /* Check if the waveform is decaying (and if so how fast) */
2045          {
2046             celt_word32 E1=1, E2=1;
2047             int period;
2048             if (pitch_index <= MAX_PERIOD/2)
2049                period = pitch_index;
2050             else
2051                period = MAX_PERIOD/2;
2052             for (i=0;i<period;i++)
2053             {
2054                E1 += SHR32(MULT16_16(exc[MAX_PERIOD-period+i],exc[MAX_PERIOD-period+i]),8);
2055                E2 += SHR32(MULT16_16(exc[MAX_PERIOD-2*period+i],exc[MAX_PERIOD-2*period+i]),8);
2056             }
2057             if (E1 > E2)
2058                E1 = E2;
2059             decay = celt_sqrt(frac_div32(SHR(E1,1),E2));
2060          }
2061
2062          /* Copy excitation, taking decay into account */
2063          for (i=0;i<len+st->mode->overlap;i++)
2064          {
2065             celt_word16 tmp;
2066             if (offset+i >= MAX_PERIOD)
2067             {
2068                offset -= pitch_index;
2069                decay = MULT16_16_Q15(decay, decay);
2070             }
2071             e[i] = SHL32(EXTEND32(MULT16_16_Q15(decay, exc[offset+i])), SIG_SHIFT);
2072             tmp = ROUND16(out_mem[c][offset+i],SIG_SHIFT);
2073             S1 += SHR32(MULT16_16(tmp,tmp),8);
2074          }
2075          for (i=0;i<LPC_ORDER;i++)
2076             mem[i] = ROUND16(out_mem[c][MAX_PERIOD-1-i], SIG_SHIFT);
2077          for (i=0;i<len+st->mode->overlap;i++)
2078             e[i] = MULT16_32_Q15(fade, e[i]);
2079          iir(e, lpc+c*LPC_ORDER, e, len+st->mode->overlap, LPC_ORDER, mem);
2080
2081          {
2082             celt_word32 S2=0;
2083             for (i=0;i<len+overlap;i++)
2084             {
2085                celt_word16 tmp = ROUND16(e[i],SIG_SHIFT);
2086                S2 += SHR32(MULT16_16(tmp,tmp),8);
2087             }
2088             /* This checks for an "explosion" in the synthesis */
2089 #ifdef FIXED_POINT
2090             if (!(S1 > SHR32(S2,2)))
2091 #else
2092                /* Float test is written this way to catch NaNs at the same time */
2093                if (!(S1 > 0.2f*S2))
2094 #endif
2095                {
2096                   for (i=0;i<len+overlap;i++)
2097                      e[i] = 0;
2098                } else if (S1 < S2)
2099                {
2100                   celt_word16 ratio = celt_sqrt(frac_div32(SHR32(S1,1)+1,S2+1));
2101                   for (i=0;i<len+overlap;i++)
2102                      e[i] = MULT16_32_Q15(ratio, e[i]);
2103                }
2104          }
2105
2106 #ifdef ENABLE_POSTFILTER
2107          /* Apply post-filter to the MDCT overlap of the previous frame */
2108          comb_filter(out_mem[c]+MAX_PERIOD, out_mem[c]+MAX_PERIOD, st->postfilter_period, st->postfilter_period, st->overlap,
2109                st->postfilter_gain, st->postfilter_gain, st->postfilter_tapset, st->postfilter_tapset,
2110                NULL, 0);
2111 #endif /* ENABLE_POSTFILTER */
2112
2113          for (i=0;i<MAX_PERIOD+st->mode->overlap-N;i++)
2114             out_mem[c][i] = out_mem[c][N+i];
2115
2116          /* Apply TDAC to the concealed audio so that it blends with the
2117          previous and next frames */
2118          for (i=0;i<overlap/2;i++)
2119          {
2120             celt_word32 tmp;
2121             tmp = MULT16_32_Q15(st->mode->window[i],           e[N+overlap-1-i]) +
2122                   MULT16_32_Q15(st->mode->window[overlap-i-1], e[N+i          ]);
2123             out_mem[c][MAX_PERIOD+i] = MULT16_32_Q15(st->mode->window[overlap-i-1], tmp);
2124             out_mem[c][MAX_PERIOD+overlap-i-1] = MULT16_32_Q15(st->mode->window[i], tmp);
2125          }
2126          for (i=0;i<N;i++)
2127             out_mem[c][MAX_PERIOD-N+i] = e[i];
2128
2129 #ifdef ENABLE_POSTFILTER
2130          /* Apply pre-filter to the MDCT overlap for the next frame (post-filter will be applied then) */
2131          comb_filter(e, out_mem[c]+MAX_PERIOD, st->postfilter_period, st->postfilter_period, st->overlap,
2132                -st->postfilter_gain, -st->postfilter_gain, st->postfilter_tapset, st->postfilter_tapset,
2133                NULL, 0);
2134 #endif /* ENABLE_POSTFILTER */
2135          for (i=0;i<overlap;i++)
2136             out_mem[c][MAX_PERIOD+i] = e[i];
2137       } while (++c<C);
2138    }
2139
2140    deemphasis(out_syn, pcm, N, C, st->downsample, st->mode->preemph, st->preemph_memD);
2141    
2142    st->loss_count++;
2143
2144    RESTORE_STACK;
2145 }
2146
2147 #ifdef FIXED_POINT
2148 CELT_STATIC
2149 int celt_decode_with_ec(CELTDecoder * restrict st, const unsigned char *data, int len, celt_int16 * restrict pcm, int frame_size, ec_dec *dec)
2150 {
2151 #else
2152 CELT_STATIC
2153 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)
2154 {
2155 #endif
2156    int c, i, N;
2157    int spread_decision;
2158    celt_int32 bits;
2159    ec_dec _dec;
2160    VARDECL(celt_sig, freq);
2161    VARDECL(celt_norm, X);
2162    VARDECL(celt_ener, bandE);
2163    VARDECL(int, fine_quant);
2164    VARDECL(int, pulses);
2165    VARDECL(int, cap);
2166    VARDECL(int, offsets);
2167    VARDECL(int, fine_priority);
2168    VARDECL(int, tf_res);
2169    VARDECL(unsigned char, collapse_masks);
2170    celt_sig *out_mem[2];
2171    celt_sig *decode_mem[2];
2172    celt_sig *overlap_mem[2];
2173    celt_sig *out_syn[2];
2174    celt_word16 *lpc;
2175    celt_word16 *oldBandE, *oldLogE, *oldLogE2, *backgroundLogE;
2176
2177    int shortBlocks;
2178    int isTransient;
2179    int intra_ener;
2180    const int CC = CHANNELS(st->channels);
2181    int LM, M;
2182    int effEnd;
2183    int codedBands;
2184    int alloc_trim;
2185    int postfilter_pitch;
2186    celt_word16 postfilter_gain;
2187    int intensity=0;
2188    int dual_stereo=0;
2189    celt_int32 total_bits;
2190    celt_int32 balance;
2191    celt_int32 tell;
2192    int dynalloc_logp;
2193    int postfilter_tapset;
2194    int anti_collapse_rsv;
2195    int anti_collapse_on=0;
2196    int silence;
2197    const int C = CHANNELS(st->stream_channels);
2198
2199    SAVE_STACK;
2200
2201    if (len<0 || len>1275 || pcm==NULL)
2202       return CELT_BAD_ARG;
2203
2204    frame_size *= st->downsample;
2205    for (LM=0;LM<=st->mode->maxLM;LM++)
2206       if (st->mode->shortMdctSize<<LM==frame_size)
2207          break;
2208    if (LM>st->mode->maxLM)
2209       return CELT_BAD_ARG;
2210    M=1<<LM;
2211
2212    c=0; do {
2213       decode_mem[c] = st->_decode_mem + c*(DECODE_BUFFER_SIZE+st->overlap);
2214       out_mem[c] = decode_mem[c]+DECODE_BUFFER_SIZE-MAX_PERIOD;
2215       overlap_mem[c] = decode_mem[c]+DECODE_BUFFER_SIZE;
2216    } while (++c<CC);
2217    lpc = (celt_word16*)(st->_decode_mem+(DECODE_BUFFER_SIZE+st->overlap)*CC);
2218    oldBandE = lpc+CC*LPC_ORDER;
2219    oldLogE = oldBandE + CC*st->mode->nbEBands;
2220    oldLogE2 = oldLogE + CC*st->mode->nbEBands;
2221    backgroundLogE = oldLogE2  + CC*st->mode->nbEBands;
2222
2223    N = M*st->mode->shortMdctSize;
2224
2225    effEnd = st->end;
2226    if (effEnd > st->mode->effEBands)
2227       effEnd = st->mode->effEBands;
2228
2229    ALLOC(freq, CC*N, celt_sig); /**< Interleaved signal MDCTs */
2230    ALLOC(X, CC*N, celt_norm);   /**< Interleaved normalised MDCTs */
2231    ALLOC(bandE, st->mode->nbEBands*CC, celt_ener);
2232    c=0; do
2233       for (i=0;i<M*st->mode->eBands[st->start];i++)
2234          X[c*N+i] = 0;
2235    while (++c<CC);
2236    c=0; do   
2237       for (i=M*st->mode->eBands[effEnd];i<N;i++)
2238          X[c*N+i] = 0;
2239    while (++c<CC);
2240
2241    if (data == NULL || len<=1)
2242    {
2243       celt_decode_lost(st, pcm, N, LM);
2244       RESTORE_STACK;
2245       return CELT_OK;
2246    }
2247    if (len<0) {
2248      RESTORE_STACK;
2249      return CELT_BAD_ARG;
2250    }
2251    
2252    if (dec == NULL)
2253    {
2254       ec_dec_init(&_dec,(unsigned char*)data,len);
2255       dec = &_dec;
2256    }
2257
2258    if (C>CC)
2259    {
2260       RESTORE_STACK;
2261       return CELT_CORRUPTED_DATA;
2262    } else if (C<CC)
2263    {
2264       for (i=0;i<st->mode->nbEBands;i++)
2265          oldBandE[i]=MAX16(oldBandE[i],oldBandE[st->mode->nbEBands+i]);
2266    }
2267
2268    total_bits = len*8;
2269    tell = ec_tell(dec);
2270
2271    if (tell==1)
2272       silence = ec_dec_bit_logp(dec, 15);
2273    else
2274       silence = 0;
2275    if (silence)
2276    {
2277       /* Pretend we've read all the remaining bits */
2278       tell = len*8;
2279       dec->nbits_total+=tell-ec_tell(dec);
2280    }
2281
2282    postfilter_gain = 0;
2283    postfilter_pitch = 0;
2284    postfilter_tapset = 0;
2285    if (st->start==0 && tell+16 <= total_bits)
2286    {
2287       if(ec_dec_bit_logp(dec, 1))
2288       {
2289 #ifdef ENABLE_POSTFILTER
2290          int qg, octave;
2291          octave = ec_dec_uint(dec, 6);
2292          postfilter_pitch = (16<<octave)+ec_dec_bits(dec, 4+octave)-1;
2293          qg = ec_dec_bits(dec, 3);
2294          if (ec_tell(dec)+2<=total_bits)
2295             postfilter_tapset = ec_dec_icdf(dec, tapset_icdf, 2);
2296          postfilter_gain = QCONST16(.09375f,15)*(qg+1);
2297 #else /* ENABLE_POSTFILTER */
2298          RESTORE_STACK;
2299          return CELT_CORRUPTED_DATA;
2300 #endif /* ENABLE_POSTFILTER */
2301       }
2302       tell = ec_tell(dec);
2303    }
2304
2305    if (LM > 0 && tell+3 <= total_bits)
2306    {
2307       isTransient = ec_dec_bit_logp(dec, 3);
2308       tell = ec_tell(dec);
2309    }
2310    else
2311       isTransient = 0;
2312
2313    if (isTransient)
2314       shortBlocks = M;
2315    else
2316       shortBlocks = 0;
2317
2318    /* Decode the global flags (first symbols in the stream) */
2319    intra_ener = tell+3<=total_bits ? ec_dec_bit_logp(dec, 3) : 0;
2320    /* Get band energies */
2321    unquant_coarse_energy(st->mode, st->start, st->end, oldBandE,
2322          intra_ener, dec, C, LM);
2323
2324    ALLOC(tf_res, st->mode->nbEBands, int);
2325    tf_decode(st->start, st->end, isTransient, tf_res, LM, dec);
2326
2327    tell = ec_tell(dec);
2328    spread_decision = SPREAD_NORMAL;
2329    if (tell+4 <= total_bits)
2330       spread_decision = ec_dec_icdf(dec, spread_icdf, 5);
2331
2332    ALLOC(pulses, st->mode->nbEBands, int);
2333    ALLOC(cap, st->mode->nbEBands, int);
2334    ALLOC(offsets, st->mode->nbEBands, int);
2335    ALLOC(fine_priority, st->mode->nbEBands, int);
2336
2337    init_caps(st->mode,cap,LM,C);
2338
2339    dynalloc_logp = 6;
2340    total_bits<<=BITRES;
2341    tell = ec_tell_frac(dec);
2342    for (i=st->start;i<st->end;i++)
2343    {
2344       int width, quanta;
2345       int dynalloc_loop_logp;
2346       int boost;
2347       width = C*(st->mode->eBands[i+1]-st->mode->eBands[i])<<LM;
2348       /* quanta is 6 bits, but no more than 1 bit/sample
2349          and no less than 1/8 bit/sample */
2350       quanta = IMIN(width<<BITRES, IMAX(6<<BITRES, width));
2351       dynalloc_loop_logp = dynalloc_logp;
2352       boost = 0;
2353       while (tell+(dynalloc_loop_logp<<BITRES) < total_bits && boost < cap[i])
2354       {
2355          int flag;
2356          flag = ec_dec_bit_logp(dec, dynalloc_loop_logp);
2357          tell = ec_tell_frac(dec);
2358          if (!flag)
2359             break;
2360          boost += quanta;
2361          total_bits -= quanta;
2362          dynalloc_loop_logp = 1;
2363       }
2364       offsets[i] = boost;
2365       /* Making dynalloc more likely */
2366       if (boost>0)
2367          dynalloc_logp = IMAX(2, dynalloc_logp-1);
2368    }
2369
2370    ALLOC(fine_quant, st->mode->nbEBands, int);
2371    alloc_trim = tell+(6<<BITRES) <= total_bits ?
2372          ec_dec_icdf(dec, trim_icdf, 7) : 5;
2373
2374    bits = ((celt_int32)len*8<<BITRES) - ec_tell_frac(dec) - 1;
2375    anti_collapse_rsv = isTransient&&LM>=2&&bits>=(LM+2<<BITRES) ? (1<<BITRES) : 0;
2376    bits -= anti_collapse_rsv;
2377    codedBands = compute_allocation(st->mode, st->start, st->end, offsets, cap,
2378          alloc_trim, &intensity, &dual_stereo, bits, &balance, pulses,
2379          fine_quant, fine_priority, C, LM, dec, 0, 0);
2380    
2381    unquant_fine_energy(st->mode, st->start, st->end, oldBandE, fine_quant, dec, C);
2382
2383    /* Decode fixed codebook */
2384    ALLOC(collapse_masks, C*st->mode->nbEBands, unsigned char);
2385    quant_all_bands(0, st->mode, st->start, st->end, X, C==2 ? X+N : NULL, collapse_masks,
2386          NULL, pulses, shortBlocks, spread_decision, dual_stereo, intensity, tf_res, 1,
2387          len*(8<<BITRES)-anti_collapse_rsv, balance, dec, LM, codedBands, &st->rng);
2388
2389    if (anti_collapse_rsv > 0)
2390    {
2391       anti_collapse_on = ec_dec_bits(dec, 1);
2392    }
2393
2394    unquant_energy_finalise(st->mode, st->start, st->end, oldBandE,
2395          fine_quant, fine_priority, len*8-ec_tell(dec), dec, C);
2396
2397    if (anti_collapse_on)
2398       anti_collapse(st->mode, X, collapse_masks, LM, C, CC, N,
2399             st->start, st->end, oldBandE, oldLogE, oldLogE2, pulses, st->rng);
2400
2401    log2Amp(st->mode, st->start, st->end, bandE, oldBandE, C);
2402
2403    if (silence)
2404    {
2405       for (i=0;i<C*st->mode->nbEBands;i++)
2406       {
2407          bandE[i] = 0;
2408          oldBandE[i] = -QCONST16(28.f,DB_SHIFT);
2409       }
2410    }
2411    /* Synthesis */
2412    denormalise_bands(st->mode, X, freq, bandE, effEnd, C, M);
2413
2414    CELT_MOVE(decode_mem[0], decode_mem[0]+N, DECODE_BUFFER_SIZE-N);
2415    if (CC==2)
2416       CELT_MOVE(decode_mem[1], decode_mem[1]+N, DECODE_BUFFER_SIZE-N);
2417
2418    c=0; do
2419       for (i=0;i<M*st->mode->eBands[st->start];i++)
2420          freq[c*N+i] = 0;
2421    while (++c<C);
2422    c=0; do {
2423       int bound = M*st->mode->eBands[effEnd];
2424       if (st->downsample!=1)
2425          bound = IMIN(bound, N/st->downsample);
2426       for (i=bound;i<N;i++)
2427          freq[c*N+i] = 0;
2428    } while (++c<C);
2429
2430    out_syn[0] = out_mem[0]+MAX_PERIOD-N;
2431    if (CC==2)
2432       out_syn[1] = out_mem[1]+MAX_PERIOD-N;
2433
2434    if (CC==2&&C==1)
2435    {
2436       for (i=0;i<N;i++)
2437          freq[N+i] = freq[i];
2438    }
2439
2440    /* Compute inverse MDCTs */
2441    compute_inv_mdcts(st->mode, shortBlocks, freq, out_syn, overlap_mem, CC, LM);
2442
2443 #ifdef ENABLE_POSTFILTER
2444    c=0; do {
2445       st->postfilter_period=IMAX(st->postfilter_period, COMBFILTER_MINPERIOD);
2446       st->postfilter_period_old=IMAX(st->postfilter_period_old, COMBFILTER_MINPERIOD);
2447       comb_filter(out_syn[c], out_syn[c], st->postfilter_period_old, st->postfilter_period, st->mode->shortMdctSize,
2448             st->postfilter_gain_old, st->postfilter_gain, st->postfilter_tapset_old, st->postfilter_tapset,
2449             st->mode->window, st->overlap);
2450       if (LM!=0)
2451          comb_filter(out_syn[c]+st->mode->shortMdctSize, out_syn[c]+st->mode->shortMdctSize, st->postfilter_period, postfilter_pitch, N-st->mode->shortMdctSize,
2452                st->postfilter_gain, postfilter_gain, st->postfilter_tapset, postfilter_tapset,
2453                st->mode->window, st->mode->overlap);
2454
2455    } while (++c<CC);
2456    st->postfilter_period_old = st->postfilter_period;
2457    st->postfilter_gain_old = st->postfilter_gain;
2458    st->postfilter_tapset_old = st->postfilter_tapset;
2459    st->postfilter_period = postfilter_pitch;
2460    st->postfilter_gain = postfilter_gain;
2461    st->postfilter_tapset = postfilter_tapset;
2462    if (LM!=0)
2463    {
2464       st->postfilter_period_old = st->postfilter_period;
2465       st->postfilter_gain_old = st->postfilter_gain;
2466       st->postfilter_tapset_old = st->postfilter_tapset;
2467    }
2468 #endif /* ENABLE_POSTFILTER */
2469
2470    if (CC==2&&C==1) {
2471       for (i=0;i<st->mode->nbEBands;i++)
2472          oldBandE[st->mode->nbEBands+i]=oldBandE[i];
2473    }
2474
2475    /* In case start or end were to change */
2476    c=0; do
2477    {
2478       for (i=0;i<st->start;i++)
2479          oldBandE[c*st->mode->nbEBands+i]=0;
2480       for (i=st->end;i<st->mode->nbEBands;i++)
2481          oldBandE[c*st->mode->nbEBands+i]=0;
2482    } while (++c<CC);
2483    if (!isTransient)
2484    {
2485       for (i=0;i<CC*st->mode->nbEBands;i++)
2486          oldLogE2[i] = oldLogE[i];
2487       for (i=0;i<CC*st->mode->nbEBands;i++)
2488          oldLogE[i] = oldBandE[i];
2489       for (i=0;i<CC*st->mode->nbEBands;i++)
2490          backgroundLogE[i] = MIN16(backgroundLogE[i] + M*QCONST16(0.001f,DB_SHIFT), oldBandE[i]);
2491    } else {
2492       for (i=0;i<CC*st->mode->nbEBands;i++)
2493          oldLogE[i] = MIN16(oldLogE[i], oldBandE[i]);
2494    }
2495    st->rng = dec->rng;
2496
2497    deemphasis(out_syn, pcm, N, CC, st->downsample, st->mode->preemph, st->preemph_memD);
2498    st->loss_count = 0;
2499    RESTORE_STACK;
2500    if (ec_tell(dec) > 8*len || ec_get_error(dec))
2501       return CELT_CORRUPTED_DATA;
2502    else
2503       return CELT_OK;
2504 }
2505
2506 #ifdef FIXED_POINT
2507 #ifndef DISABLE_FLOAT_API
2508 CELT_STATIC
2509 int celt_decode_with_ec_float(CELTDecoder * restrict st, const unsigned char *data, int len, float * restrict pcm, int frame_size, ec_dec *dec)
2510 {
2511    int j, ret, C, N;
2512    VARDECL(celt_int16, out);
2513    ALLOC_STACK;
2514    SAVE_STACK;
2515
2516    if (pcm==NULL)
2517       return CELT_BAD_ARG;
2518
2519    C = CHANNELS(st->channels);
2520    N = frame_size;
2521    
2522    ALLOC(out, C*N, celt_int16);
2523    ret=celt_decode_with_ec(st, data, len, out, frame_size, dec);
2524    if (ret==0)
2525       for (j=0;j<C*N;j++)
2526          pcm[j]=out[j]*(1.f/32768.f);
2527      
2528    RESTORE_STACK;
2529    return ret;
2530 }
2531 #endif /*DISABLE_FLOAT_API*/
2532 #else
2533 CELT_STATIC
2534 int celt_decode_with_ec(CELTDecoder * restrict st, const unsigned char *data, int len, celt_int16 * restrict pcm, int frame_size, ec_dec *dec)
2535 {
2536    int j, ret, C, N;
2537    VARDECL(celt_sig, out);
2538    ALLOC_STACK;
2539    SAVE_STACK;
2540
2541    if (pcm==NULL)
2542       return CELT_BAD_ARG;
2543
2544    C = CHANNELS(st->channels);
2545    N = frame_size;
2546    ALLOC(out, C*N, celt_sig);
2547
2548    ret=celt_decode_with_ec_float(st, data, len, out, frame_size, dec);
2549
2550    if (ret==0)
2551       for (j=0;j<C*N;j++)
2552          pcm[j] = FLOAT2INT16 (out[j]);
2553    
2554    RESTORE_STACK;
2555    return ret;
2556 }
2557 #endif
2558
2559 int celt_decode(CELTDecoder * restrict st, const unsigned char *data, int len, celt_int16 * restrict pcm, int frame_size)
2560 {
2561    return celt_decode_with_ec(st, data, len, pcm, frame_size, NULL);
2562 }
2563
2564 #ifndef DISABLE_FLOAT_API
2565 int celt_decode_float(CELTDecoder * restrict st, const unsigned char *data, int len, float * restrict pcm, int frame_size)
2566 {
2567    return celt_decode_with_ec_float(st, data, len, pcm, frame_size, NULL);
2568 }
2569 #endif /* DISABLE_FLOAT_API */
2570
2571 int celt_decoder_ctl(CELTDecoder * restrict st, int request, ...)
2572 {
2573    va_list ap;
2574
2575    va_start(ap, request);
2576    switch (request)
2577    {
2578       case CELT_GET_MODE_REQUEST:
2579       {
2580          const CELTMode ** value = va_arg(ap, const CELTMode**);
2581          if (value==0)
2582             goto bad_arg;
2583          *value=st->mode;
2584       }
2585       break;
2586       case CELT_SET_START_BAND_REQUEST:
2587       {
2588          celt_int32 value = va_arg(ap, celt_int32);
2589          if (value<0 || value>=st->mode->nbEBands)
2590             goto bad_arg;
2591          st->start = value;
2592       }
2593       break;
2594       case CELT_SET_END_BAND_REQUEST:
2595       {
2596          celt_int32 value = va_arg(ap, celt_int32);
2597          if (value<1 || value>st->mode->nbEBands)
2598             goto bad_arg;
2599          st->end = value;
2600       }
2601       break;
2602       case CELT_SET_CHANNELS_REQUEST:
2603       {
2604          celt_int32 value = va_arg(ap, celt_int32);
2605          if (value<1 || value>2)
2606             goto bad_arg;
2607          st->stream_channels = value;
2608       }
2609       break;
2610       case CELT_RESET_STATE:
2611       {
2612          CELT_MEMSET((char*)&st->DECODER_RESET_START, 0,
2613                celt_decoder_get_size_custom(st->mode, st->channels)-
2614                ((char*)&st->DECODER_RESET_START - (char*)st));
2615       }
2616       break;
2617       default:
2618          goto bad_request;
2619    }
2620    va_end(ap);
2621    return CELT_OK;
2622 bad_arg:
2623    va_end(ap);
2624    return CELT_BAD_ARG;
2625 bad_request:
2626       va_end(ap);
2627   return CELT_UNIMPLEMENTED;
2628 }
2629
2630 const char *celt_strerror(int error)
2631 {
2632    static const char *error_strings[8] = {
2633       "success",
2634       "invalid argument",
2635       "invalid mode",
2636       "internal error",
2637       "corrupted stream",
2638       "request not implemented",
2639       "invalid state",
2640       "memory allocation failed"
2641    };
2642    if (error > 0 || error < -7)
2643       return "unknown error";
2644    else 
2645       return error_strings[-error];
2646 }
2647