Relicensing under the simplified (2-clause) BSD license
[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;
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          st->prefilter_period=IMAX(st->prefilter_period, COMBFILTER_MINPERIOD);
1152          CELT_COPY(in+c*(N+st->overlap), st->in_mem+c*(st->overlap), st->overlap);
1153 #ifdef ENABLE_POSTFILTER
1154          comb_filter(in+c*(N+st->overlap)+st->overlap, pre[c]+COMBFILTER_MAXPERIOD,
1155                st->prefilter_period, pitch_index, N, -st->prefilter_gain, -gain1,
1156                st->prefilter_tapset, prefilter_tapset, st->mode->window, st->mode->overlap);
1157 #endif /* ENABLE_POSTFILTER */
1158          CELT_COPY(st->in_mem+c*(st->overlap), in+c*(N+st->overlap)+N, st->overlap);
1159
1160 #ifdef ENABLE_POSTFILTER
1161          if (N>COMBFILTER_MAXPERIOD)
1162          {
1163             CELT_MOVE(prefilter_mem+c*COMBFILTER_MAXPERIOD, pre[c]+N, COMBFILTER_MAXPERIOD);
1164          } else {
1165             CELT_MOVE(prefilter_mem+c*COMBFILTER_MAXPERIOD, prefilter_mem+c*COMBFILTER_MAXPERIOD+N, COMBFILTER_MAXPERIOD-N);
1166             CELT_MOVE(prefilter_mem+c*COMBFILTER_MAXPERIOD+COMBFILTER_MAXPERIOD-N, pre[c]+COMBFILTER_MAXPERIOD, N);
1167          }
1168 #endif /* ENABLE_POSTFILTER */
1169       } while (++c<CC);
1170
1171       RESTORE_STACK;
1172    }
1173
1174 #ifdef RESYNTH
1175    resynth = 1;
1176 #else
1177    resynth = 0;
1178 #endif
1179
1180    isTransient = 0;
1181    shortBlocks = 0;
1182    if (LM>0 && ec_tell(enc)+3<=total_bits)
1183    {
1184       if (st->complexity > 1)
1185       {
1186          isTransient = transient_analysis(in, N+st->overlap, CC,
1187                   st->overlap);
1188          if (isTransient)
1189             shortBlocks = M;
1190       }
1191       ec_enc_bit_logp(enc, isTransient, 3);
1192    }
1193
1194    ALLOC(freq, CC*N, celt_sig); /**< Interleaved signal MDCTs */
1195    ALLOC(bandE,st->mode->nbEBands*CC, celt_ener);
1196    ALLOC(bandLogE,st->mode->nbEBands*CC, celt_word16);
1197    /* Compute MDCTs */
1198    compute_mdcts(st->mode, shortBlocks, in, freq, CC, LM);
1199
1200    if (CC==2&&C==1)
1201    {
1202       for (i=0;i<N;i++)
1203          freq[i] = ADD32(HALF32(freq[i]), HALF32(freq[N+i]));
1204    }
1205    if (st->upsample != 1)
1206    {
1207       c=0; do
1208       {
1209          int bound = N/st->upsample;
1210          for (i=0;i<bound;i++)
1211             freq[c*N+i] *= st->upsample;
1212          for (;i<N;i++)
1213             freq[c*N+i] = 0;
1214       } while (++c<C);
1215    }
1216    ALLOC(X, C*N, celt_norm);         /**< Interleaved normalised MDCTs */
1217
1218    compute_band_energies(st->mode, freq, bandE, effEnd, C, M);
1219
1220    amp2Log2(st->mode, effEnd, st->end, bandE, bandLogE, C);
1221
1222    /* Band normalisation */
1223    normalise_bands(st->mode, freq, X, bandE, effEnd, C, M);
1224
1225    ALLOC(tf_res, st->mode->nbEBands, int);
1226    /* Needs to be before coarse energy quantization because otherwise the energy gets modified */
1227    tf_select = tf_analysis(st->mode, bandLogE, oldBandE, effEnd, C, isTransient, tf_res, effectiveBytes, X, N, LM, &tf_sum);
1228    for (i=effEnd;i<st->end;i++)
1229       tf_res[i] = tf_res[effEnd-1];
1230
1231    ALLOC(error, C*st->mode->nbEBands, celt_word16);
1232    quant_coarse_energy(st->mode, st->start, st->end, effEnd, bandLogE,
1233          oldBandE, total_bits, error, enc,
1234          C, LM, nbAvailableBytes, st->force_intra,
1235          &st->delayedIntra, st->complexity >= 4);
1236
1237    tf_encode(st->start, st->end, isTransient, tf_res, LM, tf_select, enc);
1238
1239    st->spread_decision = SPREAD_NORMAL;
1240    if (ec_tell(enc)+4<=total_bits)
1241    {
1242       if (shortBlocks || st->complexity < 3 || nbAvailableBytes < 10*C)
1243       {
1244          if (st->complexity == 0)
1245             st->spread_decision = SPREAD_NONE;
1246       } else {
1247          st->spread_decision = spreading_decision(st->mode, X,
1248                &st->tonal_average, st->spread_decision, &st->hf_average,
1249                &st->tapset_decision, pf_on&&!shortBlocks, effEnd, C, M);
1250       }
1251       ec_enc_icdf(enc, st->spread_decision, spread_icdf, 5);
1252    }
1253
1254    ALLOC(cap, st->mode->nbEBands, int);
1255    ALLOC(offsets, st->mode->nbEBands, int);
1256
1257    init_caps(st->mode,cap,LM,C);
1258    for (i=0;i<st->mode->nbEBands;i++)
1259       offsets[i] = 0;
1260    /* Dynamic allocation code */
1261    /* Make sure that dynamic allocation can't make us bust the budget */
1262    if (effectiveBytes > 50 && LM>=1)
1263    {
1264       int t1, t2;
1265       if (LM <= 1)
1266       {
1267          t1 = 3;
1268          t2 = 5;
1269       } else {
1270          t1 = 2;
1271          t2 = 4;
1272       }
1273       for (i=st->start+1;i<st->end-1;i++)
1274       {
1275          celt_word32 d2;
1276          d2 = 2*bandLogE[i]-bandLogE[i-1]-bandLogE[i+1];
1277          if (C==2)
1278             d2 = HALF32(d2 + 2*bandLogE[i+st->mode->nbEBands]-
1279                   bandLogE[i-1+st->mode->nbEBands]-bandLogE[i+1+st->mode->nbEBands]);
1280          if (d2 > SHL16(t1,DB_SHIFT))
1281             offsets[i] += 1;
1282          if (d2 > SHL16(t2,DB_SHIFT))
1283             offsets[i] += 1;
1284       }
1285    }
1286    dynalloc_logp = 6;
1287    total_bits<<=BITRES;
1288    total_boost = 0;
1289    tell = ec_tell_frac(enc);
1290    for (i=st->start;i<st->end;i++)
1291    {
1292       int width, quanta;
1293       int dynalloc_loop_logp;
1294       int boost;
1295       int j;
1296       width = C*(st->mode->eBands[i+1]-st->mode->eBands[i])<<LM;
1297       /* quanta is 6 bits, but no more than 1 bit/sample
1298          and no less than 1/8 bit/sample */
1299       quanta = IMIN(width<<BITRES, IMAX(6<<BITRES, width));
1300       dynalloc_loop_logp = dynalloc_logp;
1301       boost = 0;
1302       for (j = 0; tell+(dynalloc_loop_logp<<BITRES) < total_bits-total_boost
1303             && boost < cap[i]; j++)
1304       {
1305          int flag;
1306          flag = j<offsets[i];
1307          ec_enc_bit_logp(enc, flag, dynalloc_loop_logp);
1308          tell = ec_tell_frac(enc);
1309          if (!flag)
1310             break;
1311          boost += quanta;
1312          total_boost += quanta;
1313          dynalloc_loop_logp = 1;
1314       }
1315       /* Making dynalloc more likely */
1316       if (j)
1317          dynalloc_logp = IMAX(2, dynalloc_logp-1);
1318       offsets[i] = boost;
1319    }
1320    alloc_trim = 5;
1321    if (tell+(6<<BITRES) <= total_bits - total_boost)
1322    {
1323       alloc_trim = alloc_trim_analysis(st->mode, X, bandLogE,
1324             st->end, LM, C, N);
1325       ec_enc_icdf(enc, alloc_trim, trim_icdf, 7);
1326       tell = ec_tell_frac(enc);
1327    }
1328
1329    /* Variable bitrate */
1330    if (vbr_rate>0)
1331    {
1332      celt_word16 alpha;
1333      celt_int32 delta;
1334      /* The target rate in 8th bits per frame */
1335      celt_int32 target;
1336      celt_int32 min_allowed;
1337
1338      target = vbr_rate + st->vbr_offset - ((40*C+20)<<BITRES);
1339
1340      /* Shortblocks get a large boost in bitrate, but since they
1341         are uncommon long blocks are not greatly affected */
1342      if (shortBlocks || tf_sum < -2*(st->end-st->start))
1343         target = 7*target/4;
1344      else if (tf_sum < -(st->end-st->start))
1345         target = 3*target/2;
1346      else if (M > 1)
1347         target-=(target+14)/28;
1348
1349      /* The current offset is removed from the target and the space used
1350         so far is added*/
1351      target=target+tell;
1352
1353      /* In VBR mode the frame size must not be reduced so much that it would
1354          result in the encoder running out of bits.
1355         The margin of 2 bytes ensures that none of the bust-prevention logic
1356          in the decoder will have triggered so far. */
1357      min_allowed = (tell+total_boost+(1<<BITRES+3)-1>>(BITRES+3)) + 2 - nbFilledBytes;
1358
1359      nbAvailableBytes = target+(1<<(BITRES+2))>>(BITRES+3);
1360      nbAvailableBytes = IMAX(min_allowed,nbAvailableBytes);
1361      nbAvailableBytes = IMIN(nbCompressedBytes,nbAvailableBytes+nbFilledBytes) - nbFilledBytes;
1362
1363      if(silence)
1364      {
1365        nbAvailableBytes = 2;
1366        target = 2*8<<BITRES;
1367      }
1368
1369      /* By how much did we "miss" the target on that frame */
1370      delta = target - vbr_rate;
1371
1372      target=nbAvailableBytes<<(BITRES+3);
1373
1374      if (st->vbr_count < 970)
1375      {
1376         st->vbr_count++;
1377         alpha = celt_rcp(SHL32(EXTEND32(st->vbr_count+20),16));
1378      } else
1379         alpha = QCONST16(.001f,15);
1380      /* How many bits have we used in excess of what we're allowed */
1381      if (st->constrained_vbr)
1382         st->vbr_reservoir += target - vbr_rate;
1383      /*printf ("%d\n", st->vbr_reservoir);*/
1384
1385      /* Compute the offset we need to apply in order to reach the target */
1386      st->vbr_drift += (celt_int32)MULT16_32_Q15(alpha,delta-st->vbr_offset-st->vbr_drift);
1387      st->vbr_offset = -st->vbr_drift;
1388      /*printf ("%d\n", st->vbr_drift);*/
1389
1390      if (st->constrained_vbr && st->vbr_reservoir < 0)
1391      {
1392         /* We're under the min value -- increase rate */
1393         int adjust = (-st->vbr_reservoir)/(8<<BITRES);
1394         /* Unless we're just coding silence */
1395         nbAvailableBytes += silence?0:adjust;
1396         st->vbr_reservoir = 0;
1397         /*printf ("+%d\n", adjust);*/
1398      }
1399      nbCompressedBytes = IMIN(nbCompressedBytes,nbAvailableBytes+nbFilledBytes);
1400      /* This moves the raw bits to take into account the new compressed size */
1401      ec_enc_shrink(enc, nbCompressedBytes);
1402    }
1403    if (C==2)
1404    {
1405       int effectiveRate;
1406
1407       /* Always use MS for 2.5 ms frames until we can do a better analysis */
1408       if (LM!=0)
1409          dual_stereo = stereo_analysis(st->mode, X, LM, N);
1410
1411       /* Account for coarse energy */
1412       effectiveRate = (8*effectiveBytes - 80)>>LM;
1413
1414       /* effectiveRate in kb/s */
1415       effectiveRate = 2*effectiveRate/5;
1416       if (effectiveRate<35)
1417          intensity = 8;
1418       else if (effectiveRate<50)
1419          intensity = 12;
1420       else if (effectiveRate<68)
1421          intensity = 16;
1422       else if (effectiveRate<84)
1423          intensity = 18;
1424       else if (effectiveRate<102)
1425          intensity = 19;
1426       else if (effectiveRate<130)
1427          intensity = 20;
1428       else
1429          intensity = 100;
1430       intensity = IMIN(st->end,IMAX(st->start, intensity));
1431    }
1432
1433    /* Bit allocation */
1434    ALLOC(fine_quant, st->mode->nbEBands, int);
1435    ALLOC(pulses, st->mode->nbEBands, int);
1436    ALLOC(fine_priority, st->mode->nbEBands, int);
1437
1438    /* bits =           packet size                    - where we are - safety*/
1439    bits = ((celt_int32)nbCompressedBytes*8<<BITRES) - ec_tell_frac(enc) - 1;
1440    anti_collapse_rsv = isTransient&&LM>=2&&bits>=(LM+2<<BITRES) ? (1<<BITRES) : 0;
1441    bits -= anti_collapse_rsv;
1442    codedBands = compute_allocation(st->mode, st->start, st->end, offsets, cap,
1443          alloc_trim, &intensity, &dual_stereo, bits, &balance, pulses,
1444          fine_quant, fine_priority, C, LM, enc, 1, st->lastCodedBands);
1445    st->lastCodedBands = codedBands;
1446
1447    quant_fine_energy(st->mode, st->start, st->end, oldBandE, error, fine_quant, enc, C);
1448
1449 #ifdef MEASURE_NORM_MSE
1450    float X0[3000];
1451    float bandE0[60];
1452    c=0; do 
1453       for (i=0;i<N;i++)
1454          X0[i+c*N] = X[i+c*N];
1455    while (++c<C);
1456    for (i=0;i<C*st->mode->nbEBands;i++)
1457       bandE0[i] = bandE[i];
1458 #endif
1459
1460    /* Residual quantisation */
1461    ALLOC(collapse_masks, C*st->mode->nbEBands, unsigned char);
1462    quant_all_bands(1, st->mode, st->start, st->end, X, C==2 ? X+N : NULL, collapse_masks,
1463          bandE, pulses, shortBlocks, st->spread_decision, dual_stereo, intensity, tf_res, resynth,
1464          nbCompressedBytes*(8<<BITRES)-anti_collapse_rsv, balance, enc, LM, codedBands, &st->rng);
1465
1466    if (anti_collapse_rsv > 0)
1467    {
1468       anti_collapse_on = st->consec_transient<2;
1469       ec_enc_bits(enc, anti_collapse_on, 1);
1470    }
1471    quant_energy_finalise(st->mode, st->start, st->end, oldBandE, error, fine_quant, fine_priority, nbCompressedBytes*8-ec_tell(enc), enc, C);
1472
1473    if (silence)
1474    {
1475       for (i=0;i<C*st->mode->nbEBands;i++)
1476          oldBandE[i] = -QCONST16(28.f,DB_SHIFT);
1477    }
1478
1479 #ifdef RESYNTH
1480    /* Re-synthesis of the coded audio if required */
1481    if (resynth)
1482    {
1483       celt_sig *out_mem[2];
1484       celt_sig *overlap_mem[2];
1485
1486       log2Amp(st->mode, st->start, st->end, bandE, oldBandE, C);
1487       if (silence)
1488       {
1489          for (i=0;i<C*st->mode->nbEBands;i++)
1490             bandE[i] = 0;
1491       }
1492
1493 #ifdef MEASURE_NORM_MSE
1494       measure_norm_mse(st->mode, X, X0, bandE, bandE0, M, N, C);
1495 #endif
1496       if (anti_collapse_on)
1497       {
1498          anti_collapse(st->mode, X, collapse_masks, LM, C, CC, N,
1499                st->start, st->end, oldBandE, oldLogE, oldLogE2, pulses, st->rng);
1500       }
1501
1502       /* Synthesis */
1503       denormalise_bands(st->mode, X, freq, bandE, effEnd, C, M);
1504
1505       CELT_MOVE(st->syn_mem[0], st->syn_mem[0]+N, MAX_PERIOD);
1506       if (CC==2)
1507          CELT_MOVE(st->syn_mem[1], st->syn_mem[1]+N, MAX_PERIOD);
1508
1509       c=0; do
1510          for (i=0;i<M*st->mode->eBands[st->start];i++)
1511             freq[c*N+i] = 0;
1512       while (++c<C);
1513       c=0; do
1514          for (i=M*st->mode->eBands[st->end];i<N;i++)
1515             freq[c*N+i] = 0;
1516       while (++c<C);
1517
1518       if (CC==2&&C==1)
1519       {
1520          for (i=0;i<N;i++)
1521             freq[N+i] = freq[i];
1522       }
1523
1524       out_mem[0] = st->syn_mem[0]+MAX_PERIOD;
1525       if (CC==2)
1526          out_mem[1] = st->syn_mem[1]+MAX_PERIOD;
1527
1528       c=0; do
1529          overlap_mem[c] = _overlap_mem + c*st->overlap;
1530       while (++c<CC);
1531
1532       compute_inv_mdcts(st->mode, shortBlocks, freq, out_mem, overlap_mem, CC, LM);
1533
1534 #ifdef ENABLE_POSTFILTER
1535       c=0; do {
1536          st->prefilter_period=IMAX(st->prefilter_period, COMBFILTER_MINPERIOD);
1537          st->prefilter_period_old=IMAX(st->prefilter_period_old, COMBFILTER_MINPERIOD);
1538          if (LM!=0)
1539          {
1540             comb_filter(out_mem[c], out_mem[c], st->prefilter_period, st->prefilter_period, st->overlap,
1541                   st->prefilter_gain, st->prefilter_gain, st->prefilter_tapset, st->prefilter_tapset,
1542                   NULL, 0);
1543             comb_filter(out_mem[c]+st->overlap, out_mem[c]+st->overlap, st->prefilter_period, pitch_index, N-st->overlap,
1544                   st->prefilter_gain, gain1, st->prefilter_tapset, prefilter_tapset,
1545                   st->mode->window, st->mode->overlap);
1546          } else {
1547             comb_filter(out_mem[c], out_mem[c], st->prefilter_period_old, st->prefilter_period, N,
1548                   st->prefilter_gain_old, st->prefilter_gain, st->prefilter_tapset_old, st->prefilter_tapset,
1549                   st->mode->window, st->mode->overlap);
1550          }
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
1565    if (CC==2&&C==1) {
1566       for (i=0;i<st->mode->nbEBands;i++)
1567          oldBandE[st->mode->nbEBands+i]=oldBandE[i];
1568    }
1569
1570    /* In case start or end were to change */
1571    c=0; do
1572    {
1573       for (i=0;i<st->start;i++)
1574          oldBandE[c*st->mode->nbEBands+i]=0;
1575       for (i=st->end;i<st->mode->nbEBands;i++)
1576          oldBandE[c*st->mode->nbEBands+i]=0;
1577    } while (++c<CC);
1578    if (!isTransient)
1579    {
1580       for (i=0;i<CC*st->mode->nbEBands;i++)
1581          oldLogE2[i] = oldLogE[i];
1582       for (i=0;i<CC*st->mode->nbEBands;i++)
1583          oldLogE[i] = oldBandE[i];
1584    } else {
1585       for (i=0;i<CC*st->mode->nbEBands;i++)
1586          oldLogE[i] = MIN16(oldLogE[i], oldBandE[i]);
1587    }
1588    if (isTransient)
1589       st->consec_transient++;
1590    else
1591       st->consec_transient=0;
1592    st->rng = enc->rng;
1593
1594    /* If there's any room left (can only happen for very high rates),
1595       it's already filled with zeros */
1596    ec_enc_done(enc);
1597    
1598    RESTORE_STACK;
1599    if (ec_get_error(enc))
1600       return CELT_CORRUPTED_DATA;
1601    else
1602       return nbCompressedBytes;
1603 }
1604
1605 #ifdef FIXED_POINT
1606 #ifndef DISABLE_FLOAT_API
1607 CELT_STATIC
1608 int celt_encode_with_ec_float(CELTEncoder * restrict st, const float * pcm, int frame_size, unsigned char *compressed, int nbCompressedBytes, ec_enc *enc)
1609 {
1610    int j, ret, C, N;
1611    VARDECL(celt_int16, in);
1612    ALLOC_STACK;
1613    SAVE_STACK;
1614
1615    if (pcm==NULL)
1616       return CELT_BAD_ARG;
1617
1618    C = CHANNELS(st->channels);
1619    N = frame_size;
1620    ALLOC(in, C*N, celt_int16);
1621
1622    for (j=0;j<C*N;j++)
1623      in[j] = FLOAT2INT16(pcm[j]);
1624
1625    ret=celt_encode_with_ec(st,in,frame_size,compressed,nbCompressedBytes, enc);
1626 #ifdef RESYNTH
1627    for (j=0;j<C*N;j++)
1628       ((float*)pcm)[j]=in[j]*(1.f/32768.f);
1629 #endif
1630    RESTORE_STACK;
1631    return ret;
1632
1633 }
1634 #endif /*DISABLE_FLOAT_API*/
1635 #else
1636 CELT_STATIC
1637 int celt_encode_with_ec(CELTEncoder * restrict st, const celt_int16 * pcm, int frame_size, unsigned char *compressed, int nbCompressedBytes, ec_enc *enc)
1638 {
1639    int j, ret, C, N;
1640    VARDECL(celt_sig, in);
1641    ALLOC_STACK;
1642    SAVE_STACK;
1643
1644    if (pcm==NULL)
1645       return CELT_BAD_ARG;
1646
1647    C=CHANNELS(st->channels);
1648    N=frame_size;
1649    ALLOC(in, C*N, celt_sig);
1650    for (j=0;j<C*N;j++) {
1651      in[j] = SCALEOUT(pcm[j]);
1652    }
1653
1654    ret = celt_encode_with_ec_float(st,in,frame_size,compressed,nbCompressedBytes, enc);
1655 #ifdef RESYNTH
1656    for (j=0;j<C*N;j++)
1657       ((celt_int16*)pcm)[j] = FLOAT2INT16(in[j]);
1658 #endif
1659    RESTORE_STACK;
1660    return ret;
1661 }
1662 #endif
1663
1664 int celt_encode(CELTEncoder * restrict st, const celt_int16 * pcm, int frame_size, unsigned char *compressed, int nbCompressedBytes)
1665 {
1666    return celt_encode_with_ec(st, pcm, frame_size, compressed, nbCompressedBytes, NULL);
1667 }
1668
1669 #ifndef DISABLE_FLOAT_API
1670 int celt_encode_float(CELTEncoder * restrict st, const float * pcm, int frame_size, unsigned char *compressed, int nbCompressedBytes)
1671 {
1672    return celt_encode_with_ec_float(st, pcm, frame_size, compressed, nbCompressedBytes, NULL);
1673 }
1674 #endif /* DISABLE_FLOAT_API */
1675
1676 int celt_encoder_ctl(CELTEncoder * restrict st, int request, ...)
1677 {
1678    va_list ap;
1679    
1680    va_start(ap, request);
1681    switch (request)
1682    {
1683       case CELT_GET_MODE_REQUEST:
1684       {
1685          const CELTMode ** value = va_arg(ap, const CELTMode**);
1686          if (value==0)
1687             goto bad_arg;
1688          *value=st->mode;
1689       }
1690       break;
1691       case CELT_SET_COMPLEXITY_REQUEST:
1692       {
1693          int value = va_arg(ap, celt_int32);
1694          if (value<0 || value>10)
1695             goto bad_arg;
1696          st->complexity = value;
1697       }
1698       break;
1699       case CELT_SET_START_BAND_REQUEST:
1700       {
1701          celt_int32 value = va_arg(ap, celt_int32);
1702          if (value<0 || value>=st->mode->nbEBands)
1703             goto bad_arg;
1704          st->start = value;
1705       }
1706       break;
1707       case CELT_SET_END_BAND_REQUEST:
1708       {
1709          celt_int32 value = va_arg(ap, celt_int32);
1710          if (value<1 || value>st->mode->nbEBands)
1711             goto bad_arg;
1712          st->end = value;
1713       }
1714       break;
1715       case CELT_SET_PREDICTION_REQUEST:
1716       {
1717          int value = va_arg(ap, celt_int32);
1718          if (value<0 || value>2)
1719             goto bad_arg;
1720          st->disable_pf = value<=1;
1721          st->force_intra = value==0;
1722       }
1723       break;
1724       case CELT_SET_VBR_CONSTRAINT_REQUEST:
1725       {
1726          celt_int32 value = va_arg(ap, celt_int32);
1727          st->constrained_vbr = value;
1728       }
1729       break;
1730       case CELT_SET_VBR_REQUEST:
1731       {
1732          celt_int32 value = va_arg(ap, celt_int32);
1733          st->vbr = value;
1734       }
1735       break;
1736       case CELT_SET_BITRATE_REQUEST:
1737       {
1738          celt_int32 value = va_arg(ap, celt_int32);
1739          if (value<=500)
1740             goto bad_arg;
1741          value = IMIN(value, 260000*st->channels);
1742          st->bitrate = value;
1743       }
1744       break;
1745       case CELT_SET_CHANNELS_REQUEST:
1746       {
1747          celt_int32 value = va_arg(ap, celt_int32);
1748          if (value<1 || value>2)
1749             goto bad_arg;
1750          st->stream_channels = value;
1751       }
1752       break;
1753       case CELT_RESET_STATE:
1754       {
1755          CELT_MEMSET((char*)&st->ENCODER_RESET_START, 0,
1756                celt_encoder_get_size_custom(st->mode, st->channels)-
1757                ((char*)&st->ENCODER_RESET_START - (char*)st));
1758          st->vbr_offset = 0;
1759          st->delayedIntra = 1;
1760          st->spread_decision = SPREAD_NORMAL;
1761          st->tonal_average = QCONST16(1.f,8);
1762       }
1763       break;
1764       case CELT_SET_INPUT_CLIPPING_REQUEST:
1765       {
1766          celt_int32 value = va_arg(ap, celt_int32);
1767          st->clip = value;
1768       }
1769       break;
1770       default:
1771          goto bad_request;
1772    }
1773    va_end(ap);
1774    return CELT_OK;
1775 bad_arg:
1776    va_end(ap);
1777    return CELT_BAD_ARG;
1778 bad_request:
1779    va_end(ap);
1780    return CELT_UNIMPLEMENTED;
1781 }
1782
1783 /**********************************************************************/
1784 /*                                                                    */
1785 /*                             DECODER                                */
1786 /*                                                                    */
1787 /**********************************************************************/
1788 #define DECODE_BUFFER_SIZE 2048
1789
1790 /** Decoder state 
1791  @brief Decoder state
1792  */
1793 struct CELTDecoder {
1794    const CELTMode *mode;
1795    int overlap;
1796    int channels;
1797    int stream_channels;
1798
1799    int downsample;
1800    int start, end;
1801
1802    /* Everything beyond this point gets cleared on a reset */
1803 #define DECODER_RESET_START rng
1804
1805    ec_uint32 rng;
1806    int last_pitch_index;
1807    int loss_count;
1808    int postfilter_period;
1809    int postfilter_period_old;
1810    celt_word16 postfilter_gain;
1811    celt_word16 postfilter_gain_old;
1812    int postfilter_tapset;
1813    int postfilter_tapset_old;
1814
1815    celt_sig preemph_memD[2];
1816    
1817    celt_sig _decode_mem[1]; /* Size = channels*(DECODE_BUFFER_SIZE+mode->overlap) */
1818    /* celt_word16 lpc[],  Size = channels*LPC_ORDER */
1819    /* celt_word16 oldEBands[], Size = channels*mode->nbEBands */
1820    /* celt_word16 oldLogE[], Size = channels*mode->nbEBands */
1821    /* celt_word16 oldLogE2[], Size = channels*mode->nbEBands */
1822    /* celt_word16 backgroundLogE[], Size = channels*mode->nbEBands */
1823 };
1824
1825 int celt_decoder_get_size(int channels)
1826 {
1827    const CELTMode *mode = celt_mode_create(48000, 960, NULL);
1828    return celt_decoder_get_size_custom(mode, channels);
1829 }
1830
1831 int celt_decoder_get_size_custom(const CELTMode *mode, int channels)
1832 {
1833    int size = sizeof(struct CELTDecoder)
1834             + (channels*(DECODE_BUFFER_SIZE+mode->overlap)-1)*sizeof(celt_sig)
1835             + channels*LPC_ORDER*sizeof(celt_word16)
1836             + 4*channels*mode->nbEBands*sizeof(celt_word16);
1837    return size;
1838 }
1839
1840 CELTDecoder *celt_decoder_create(int sampling_rate, int channels, int *error)
1841 {
1842    CELTDecoder *st;
1843    st = (CELTDecoder *)celt_alloc(celt_decoder_get_size(channels));
1844    if (st!=NULL && celt_decoder_init(st, sampling_rate, channels, error)==NULL)
1845    {
1846       celt_decoder_destroy(st);
1847       st = NULL;
1848    }
1849    return st;
1850 }
1851
1852 CELTDecoder *celt_decoder_create_custom(const CELTMode *mode, int channels, int *error)
1853 {
1854    CELTDecoder *st = (CELTDecoder *)celt_alloc(celt_decoder_get_size_custom(mode, channels));
1855    if (st!=NULL && celt_decoder_init_custom(st, mode, channels, error)==NULL)
1856    {
1857       celt_decoder_destroy(st);
1858       st = NULL;
1859    }
1860    return st;
1861 }
1862
1863 CELTDecoder *celt_decoder_init(CELTDecoder *st, int sampling_rate, int channels, int *error)
1864 {
1865    celt_decoder_init_custom(st, celt_mode_create(48000, 960, NULL), channels, error);
1866    st->downsample = resampling_factor(sampling_rate);
1867    if (st->downsample==0)
1868    {
1869       if (error)
1870          *error = CELT_BAD_ARG;
1871       return NULL;
1872    }
1873    return st;
1874 }
1875
1876 CELTDecoder *celt_decoder_init_custom(CELTDecoder *st, const CELTMode *mode, int channels, int *error)
1877 {
1878    if (channels < 0 || channels > 2)
1879    {
1880       if (error)
1881          *error = CELT_BAD_ARG;
1882       return NULL;
1883    }
1884
1885    if (st==NULL)
1886    {
1887       if (error)
1888          *error = CELT_ALLOC_FAIL;
1889       return NULL;
1890    }
1891
1892    CELT_MEMSET((char*)st, 0, celt_decoder_get_size_custom(mode, channels));
1893
1894    st->mode = mode;
1895    st->overlap = mode->overlap;
1896    st->stream_channels = st->channels = channels;
1897
1898    st->downsample = 1;
1899    st->start = 0;
1900    st->end = st->mode->effEBands;
1901
1902    st->loss_count = 0;
1903
1904    if (error)
1905       *error = CELT_OK;
1906    return st;
1907 }
1908
1909 void celt_decoder_destroy(CELTDecoder *st)
1910 {
1911    celt_free(st);
1912 }
1913
1914 static void celt_decode_lost(CELTDecoder * restrict st, celt_word16 * restrict pcm, int N, int LM)
1915 {
1916    int c;
1917    int pitch_index;
1918    int overlap = st->mode->overlap;
1919    celt_word16 fade = Q15ONE;
1920    int i, len;
1921    const int C = CHANNELS(st->channels);
1922    int offset;
1923    celt_sig *out_mem[2];
1924    celt_sig *decode_mem[2];
1925    celt_sig *overlap_mem[2];
1926    celt_word16 *lpc;
1927    celt_word32 *out_syn[2];
1928    celt_word16 *oldBandE, *oldLogE2, *backgroundLogE;
1929    int plc=1;
1930    SAVE_STACK;
1931    
1932    c=0; do {
1933       decode_mem[c] = st->_decode_mem + c*(DECODE_BUFFER_SIZE+st->overlap);
1934       out_mem[c] = decode_mem[c]+DECODE_BUFFER_SIZE-MAX_PERIOD;
1935       overlap_mem[c] = decode_mem[c]+DECODE_BUFFER_SIZE;
1936    } while (++c<C);
1937    lpc = (celt_word16*)(st->_decode_mem+(DECODE_BUFFER_SIZE+st->overlap)*C);
1938    oldBandE = lpc+C*LPC_ORDER;
1939    oldLogE2 = oldBandE + C*st->mode->nbEBands;
1940    backgroundLogE = oldLogE2  + C*st->mode->nbEBands;
1941
1942    out_syn[0] = out_mem[0]+MAX_PERIOD-N;
1943    if (C==2)
1944       out_syn[1] = out_mem[1]+MAX_PERIOD-N;
1945
1946    len = N+st->mode->overlap;
1947    
1948    if (st->loss_count >= 5)
1949    {
1950       VARDECL(celt_sig, freq);
1951       VARDECL(celt_norm, X);
1952       VARDECL(celt_ener, bandE);
1953       celt_uint32 seed;
1954
1955       ALLOC(freq, C*N, celt_sig); /**< Interleaved signal MDCTs */
1956       ALLOC(X, C*N, celt_norm);   /**< Interleaved normalised MDCTs */
1957       ALLOC(bandE, st->mode->nbEBands*C, celt_ener);
1958
1959       log2Amp(st->mode, st->start, st->end, bandE, backgroundLogE, C);
1960
1961       seed = st->rng;
1962       for (i=0;i<C*N;i++)
1963       {
1964             seed = lcg_rand(seed);
1965             X[i] = (celt_int32)(seed)>>20;
1966       }
1967       st->rng = seed;
1968       for (c=0;c<C;c++)
1969          for (i=0;i<st->mode->nbEBands;i++)
1970             renormalise_vector(X+N*c+(st->mode->eBands[i]<<LM), (st->mode->eBands[i+1]-st->mode->eBands[i])<<LM, Q15ONE);
1971
1972       denormalise_bands(st->mode, X, freq, bandE, st->mode->nbEBands, C, 1<<LM);
1973
1974       compute_inv_mdcts(st->mode, 0, freq, out_syn, overlap_mem, C, LM);
1975       plc = 0;
1976    } else if (st->loss_count == 0)
1977    {
1978       celt_word16 pitch_buf[MAX_PERIOD>>1];
1979       int len2 = len;
1980       /* FIXME: This is a kludge */
1981       if (len2>MAX_PERIOD>>1)
1982          len2 = MAX_PERIOD>>1;
1983       pitch_downsample(out_mem, pitch_buf, MAX_PERIOD, C);
1984       pitch_search(pitch_buf+((MAX_PERIOD-len2)>>1), pitch_buf, len2,
1985                    MAX_PERIOD-len2-100, &pitch_index);
1986       pitch_index = MAX_PERIOD-len2-pitch_index;
1987       st->last_pitch_index = pitch_index;
1988    } else {
1989       pitch_index = st->last_pitch_index;
1990       fade = QCONST16(.8f,15);
1991    }
1992
1993    if (plc)
1994    {
1995       c=0; do {
1996          /* FIXME: This is more memory than necessary */
1997          celt_word32 e[2*MAX_PERIOD];
1998          celt_word16 exc[2*MAX_PERIOD];
1999          celt_word32 ac[LPC_ORDER+1];
2000          celt_word16 decay = 1;
2001          celt_word32 S1=0;
2002          celt_word16 mem[LPC_ORDER]={0};
2003
2004          offset = MAX_PERIOD-pitch_index;
2005          for (i=0;i<MAX_PERIOD;i++)
2006             exc[i] = ROUND16(out_mem[c][i], SIG_SHIFT);
2007
2008          if (st->loss_count == 0)
2009          {
2010             _celt_autocorr(exc, ac, st->mode->window, st->mode->overlap,
2011                   LPC_ORDER, MAX_PERIOD);
2012
2013             /* Noise floor -40 dB */
2014 #ifdef FIXED_POINT
2015             ac[0] += SHR32(ac[0],13);
2016 #else
2017             ac[0] *= 1.0001f;
2018 #endif
2019             /* Lag windowing */
2020             for (i=1;i<=LPC_ORDER;i++)
2021             {
2022                /*ac[i] *= exp(-.5*(2*M_PI*.002*i)*(2*M_PI*.002*i));*/
2023 #ifdef FIXED_POINT
2024                ac[i] -= MULT16_32_Q15(2*i*i, ac[i]);
2025 #else
2026                ac[i] -= ac[i]*(.008f*i)*(.008f*i);
2027 #endif
2028             }
2029
2030             _celt_lpc(lpc+c*LPC_ORDER, ac, LPC_ORDER);
2031          }
2032          for (i=0;i<LPC_ORDER;i++)
2033             mem[i] = ROUND16(out_mem[c][MAX_PERIOD-1-i], SIG_SHIFT);
2034          fir(exc, lpc+c*LPC_ORDER, exc, MAX_PERIOD, LPC_ORDER, mem);
2035          /*for (i=0;i<MAX_PERIOD;i++)printf("%d ", exc[i]); printf("\n");*/
2036          /* Check if the waveform is decaying (and if so how fast) */
2037          {
2038             celt_word32 E1=1, E2=1;
2039             int period;
2040             if (pitch_index <= MAX_PERIOD/2)
2041                period = pitch_index;
2042             else
2043                period = MAX_PERIOD/2;
2044             for (i=0;i<period;i++)
2045             {
2046                E1 += SHR32(MULT16_16(exc[MAX_PERIOD-period+i],exc[MAX_PERIOD-period+i]),8);
2047                E2 += SHR32(MULT16_16(exc[MAX_PERIOD-2*period+i],exc[MAX_PERIOD-2*period+i]),8);
2048             }
2049             if (E1 > E2)
2050                E1 = E2;
2051             decay = celt_sqrt(frac_div32(SHR(E1,1),E2));
2052          }
2053
2054          /* Copy excitation, taking decay into account */
2055          for (i=0;i<len+st->mode->overlap;i++)
2056          {
2057             celt_word16 tmp;
2058             if (offset+i >= MAX_PERIOD)
2059             {
2060                offset -= pitch_index;
2061                decay = MULT16_16_Q15(decay, decay);
2062             }
2063             e[i] = SHL32(EXTEND32(MULT16_16_Q15(decay, exc[offset+i])), SIG_SHIFT);
2064             tmp = ROUND16(out_mem[c][offset+i],SIG_SHIFT);
2065             S1 += SHR32(MULT16_16(tmp,tmp),8);
2066          }
2067          for (i=0;i<LPC_ORDER;i++)
2068             mem[i] = ROUND16(out_mem[c][MAX_PERIOD-1-i], SIG_SHIFT);
2069          for (i=0;i<len+st->mode->overlap;i++)
2070             e[i] = MULT16_32_Q15(fade, e[i]);
2071          iir(e, lpc+c*LPC_ORDER, e, len+st->mode->overlap, LPC_ORDER, mem);
2072
2073          {
2074             celt_word32 S2=0;
2075             for (i=0;i<len+overlap;i++)
2076             {
2077                celt_word16 tmp = ROUND16(e[i],SIG_SHIFT);
2078                S2 += SHR32(MULT16_16(tmp,tmp),8);
2079             }
2080             /* This checks for an "explosion" in the synthesis */
2081 #ifdef FIXED_POINT
2082             if (!(S1 > SHR32(S2,2)))
2083 #else
2084                /* Float test is written this way to catch NaNs at the same time */
2085                if (!(S1 > 0.2f*S2))
2086 #endif
2087                {
2088                   for (i=0;i<len+overlap;i++)
2089                      e[i] = 0;
2090                } else if (S1 < S2)
2091                {
2092                   celt_word16 ratio = celt_sqrt(frac_div32(SHR32(S1,1)+1,S2+1));
2093                   for (i=0;i<len+overlap;i++)
2094                      e[i] = MULT16_32_Q15(ratio, e[i]);
2095                }
2096          }
2097
2098 #ifdef ENABLE_POSTFILTER
2099          /* Apply post-filter to the MDCT overlap of the previous frame */
2100          comb_filter(out_mem[c]+MAX_PERIOD, out_mem[c]+MAX_PERIOD, st->postfilter_period, st->postfilter_period, st->overlap,
2101                st->postfilter_gain, st->postfilter_gain, st->postfilter_tapset, st->postfilter_tapset,
2102                NULL, 0);
2103 #endif /* ENABLE_POSTFILTER */
2104
2105          for (i=0;i<MAX_PERIOD+st->mode->overlap-N;i++)
2106             out_mem[c][i] = out_mem[c][N+i];
2107
2108          /* Apply TDAC to the concealed audio so that it blends with the
2109          previous and next frames */
2110          for (i=0;i<overlap/2;i++)
2111          {
2112             celt_word32 tmp;
2113             tmp = MULT16_32_Q15(st->mode->window[i],           e[N+overlap-1-i]) +
2114                   MULT16_32_Q15(st->mode->window[overlap-i-1], e[N+i          ]);
2115             out_mem[c][MAX_PERIOD+i] = MULT16_32_Q15(st->mode->window[overlap-i-1], tmp);
2116             out_mem[c][MAX_PERIOD+overlap-i-1] = MULT16_32_Q15(st->mode->window[i], tmp);
2117          }
2118          for (i=0;i<N;i++)
2119             out_mem[c][MAX_PERIOD-N+i] = e[i];
2120
2121 #ifdef ENABLE_POSTFILTER
2122          /* Apply pre-filter to the MDCT overlap for the next frame (post-filter will be applied then) */
2123          comb_filter(e, out_mem[c]+MAX_PERIOD, st->postfilter_period, st->postfilter_period, st->overlap,
2124                -st->postfilter_gain, -st->postfilter_gain, st->postfilter_tapset, st->postfilter_tapset,
2125                NULL, 0);
2126 #endif /* ENABLE_POSTFILTER */
2127          for (i=0;i<overlap;i++)
2128             out_mem[c][MAX_PERIOD+i] = e[i];
2129       } while (++c<C);
2130    }
2131
2132    deemphasis(out_syn, pcm, N, C, st->downsample, st->mode->preemph, st->preemph_memD);
2133    
2134    st->loss_count++;
2135
2136    RESTORE_STACK;
2137 }
2138
2139 #ifdef FIXED_POINT
2140 CELT_STATIC
2141 int celt_decode_with_ec(CELTDecoder * restrict st, const unsigned char *data, int len, celt_int16 * restrict pcm, int frame_size, ec_dec *dec)
2142 {
2143 #else
2144 CELT_STATIC
2145 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)
2146 {
2147 #endif
2148    int c, i, N;
2149    int spread_decision;
2150    celt_int32 bits;
2151    ec_dec _dec;
2152    VARDECL(celt_sig, freq);
2153    VARDECL(celt_norm, X);
2154    VARDECL(celt_ener, bandE);
2155    VARDECL(int, fine_quant);
2156    VARDECL(int, pulses);
2157    VARDECL(int, cap);
2158    VARDECL(int, offsets);
2159    VARDECL(int, fine_priority);
2160    VARDECL(int, tf_res);
2161    VARDECL(unsigned char, collapse_masks);
2162    celt_sig *out_mem[2];
2163    celt_sig *decode_mem[2];
2164    celt_sig *overlap_mem[2];
2165    celt_sig *out_syn[2];
2166    celt_word16 *lpc;
2167    celt_word16 *oldBandE, *oldLogE, *oldLogE2, *backgroundLogE;
2168
2169    int shortBlocks;
2170    int isTransient;
2171    int intra_ener;
2172    const int CC = CHANNELS(st->channels);
2173    int LM, M;
2174    int effEnd;
2175    int codedBands;
2176    int alloc_trim;
2177    int postfilter_pitch;
2178    celt_word16 postfilter_gain;
2179    int intensity=0;
2180    int dual_stereo=0;
2181    celt_int32 total_bits;
2182    celt_int32 balance;
2183    celt_int32 tell;
2184    int dynalloc_logp;
2185    int postfilter_tapset;
2186    int anti_collapse_rsv;
2187    int anti_collapse_on=0;
2188    int silence;
2189    const int C = CHANNELS(st->stream_channels);
2190
2191    SAVE_STACK;
2192
2193    if (len<0 || len>1275 || pcm==NULL)
2194       return CELT_BAD_ARG;
2195
2196    frame_size *= st->downsample;
2197    for (LM=0;LM<=st->mode->maxLM;LM++)
2198       if (st->mode->shortMdctSize<<LM==frame_size)
2199          break;
2200    if (LM>st->mode->maxLM)
2201       return CELT_BAD_ARG;
2202    M=1<<LM;
2203
2204    c=0; do {
2205       decode_mem[c] = st->_decode_mem + c*(DECODE_BUFFER_SIZE+st->overlap);
2206       out_mem[c] = decode_mem[c]+DECODE_BUFFER_SIZE-MAX_PERIOD;
2207       overlap_mem[c] = decode_mem[c]+DECODE_BUFFER_SIZE;
2208    } while (++c<CC);
2209    lpc = (celt_word16*)(st->_decode_mem+(DECODE_BUFFER_SIZE+st->overlap)*CC);
2210    oldBandE = lpc+CC*LPC_ORDER;
2211    oldLogE = oldBandE + CC*st->mode->nbEBands;
2212    oldLogE2 = oldLogE + CC*st->mode->nbEBands;
2213    backgroundLogE = oldLogE2  + CC*st->mode->nbEBands;
2214
2215    N = M*st->mode->shortMdctSize;
2216
2217    effEnd = st->end;
2218    if (effEnd > st->mode->effEBands)
2219       effEnd = st->mode->effEBands;
2220
2221    ALLOC(freq, CC*N, celt_sig); /**< Interleaved signal MDCTs */
2222    ALLOC(X, CC*N, celt_norm);   /**< Interleaved normalised MDCTs */
2223    ALLOC(bandE, st->mode->nbEBands*CC, celt_ener);
2224    c=0; do
2225       for (i=0;i<M*st->mode->eBands[st->start];i++)
2226          X[c*N+i] = 0;
2227    while (++c<CC);
2228    c=0; do   
2229       for (i=M*st->mode->eBands[effEnd];i<N;i++)
2230          X[c*N+i] = 0;
2231    while (++c<CC);
2232
2233    if (data == NULL || len<=1)
2234    {
2235       celt_decode_lost(st, pcm, N, LM);
2236       RESTORE_STACK;
2237       return CELT_OK;
2238    }
2239    if (len<0) {
2240      RESTORE_STACK;
2241      return CELT_BAD_ARG;
2242    }
2243    
2244    if (dec == NULL)
2245    {
2246       ec_dec_init(&_dec,(unsigned char*)data,len);
2247       dec = &_dec;
2248    }
2249
2250    if (C>CC)
2251    {
2252       RESTORE_STACK;
2253       return CELT_CORRUPTED_DATA;
2254    } else if (C<CC)
2255    {
2256       for (i=0;i<st->mode->nbEBands;i++)
2257          oldBandE[i]=MAX16(oldBandE[i],oldBandE[st->mode->nbEBands+i]);
2258    }
2259
2260    total_bits = len*8;
2261    tell = ec_tell(dec);
2262
2263    if (tell==1)
2264       silence = ec_dec_bit_logp(dec, 15);
2265    else
2266       silence = 0;
2267    if (silence)
2268    {
2269       /* Pretend we've read all the remaining bits */
2270       tell = len*8;
2271       dec->nbits_total+=tell-ec_tell(dec);
2272    }
2273
2274    postfilter_gain = 0;
2275    postfilter_pitch = 0;
2276    postfilter_tapset = 0;
2277    if (st->start==0 && tell+16 <= total_bits)
2278    {
2279       if(ec_dec_bit_logp(dec, 1))
2280       {
2281 #ifdef ENABLE_POSTFILTER
2282          int qg, octave;
2283          octave = ec_dec_uint(dec, 6);
2284          postfilter_pitch = (16<<octave)+ec_dec_bits(dec, 4+octave)-1;
2285          qg = ec_dec_bits(dec, 3);
2286          if (ec_tell(dec)+2<=total_bits)
2287             postfilter_tapset = ec_dec_icdf(dec, tapset_icdf, 2);
2288          postfilter_gain = QCONST16(.09375f,15)*(qg+1);
2289 #else /* ENABLE_POSTFILTER */
2290          RESTORE_STACK;
2291          return CELT_CORRUPTED_DATA;
2292 #endif /* ENABLE_POSTFILTER */
2293       }
2294       tell = ec_tell(dec);
2295    }
2296
2297    if (LM > 0 && tell+3 <= total_bits)
2298    {
2299       isTransient = ec_dec_bit_logp(dec, 3);
2300       tell = ec_tell(dec);
2301    }
2302    else
2303       isTransient = 0;
2304
2305    if (isTransient)
2306       shortBlocks = M;
2307    else
2308       shortBlocks = 0;
2309
2310    /* Decode the global flags (first symbols in the stream) */
2311    intra_ener = tell+3<=total_bits ? ec_dec_bit_logp(dec, 3) : 0;
2312    /* Get band energies */
2313    unquant_coarse_energy(st->mode, st->start, st->end, oldBandE,
2314          intra_ener, dec, C, LM);
2315
2316    ALLOC(tf_res, st->mode->nbEBands, int);
2317    tf_decode(st->start, st->end, isTransient, tf_res, LM, dec);
2318
2319    tell = ec_tell(dec);
2320    spread_decision = SPREAD_NORMAL;
2321    if (tell+4 <= total_bits)
2322       spread_decision = ec_dec_icdf(dec, spread_icdf, 5);
2323
2324    ALLOC(pulses, st->mode->nbEBands, int);
2325    ALLOC(cap, st->mode->nbEBands, int);
2326    ALLOC(offsets, st->mode->nbEBands, int);
2327    ALLOC(fine_priority, st->mode->nbEBands, int);
2328
2329    init_caps(st->mode,cap,LM,C);
2330
2331    dynalloc_logp = 6;
2332    total_bits<<=BITRES;
2333    tell = ec_tell_frac(dec);
2334    for (i=st->start;i<st->end;i++)
2335    {
2336       int width, quanta;
2337       int dynalloc_loop_logp;
2338       int boost;
2339       width = C*(st->mode->eBands[i+1]-st->mode->eBands[i])<<LM;
2340       /* quanta is 6 bits, but no more than 1 bit/sample
2341          and no less than 1/8 bit/sample */
2342       quanta = IMIN(width<<BITRES, IMAX(6<<BITRES, width));
2343       dynalloc_loop_logp = dynalloc_logp;
2344       boost = 0;
2345       while (tell+(dynalloc_loop_logp<<BITRES) < total_bits && boost < cap[i])
2346       {
2347          int flag;
2348          flag = ec_dec_bit_logp(dec, dynalloc_loop_logp);
2349          tell = ec_tell_frac(dec);
2350          if (!flag)
2351             break;
2352          boost += quanta;
2353          total_bits -= quanta;
2354          dynalloc_loop_logp = 1;
2355       }
2356       offsets[i] = boost;
2357       /* Making dynalloc more likely */
2358       if (boost>0)
2359          dynalloc_logp = IMAX(2, dynalloc_logp-1);
2360    }
2361
2362    ALLOC(fine_quant, st->mode->nbEBands, int);
2363    alloc_trim = tell+(6<<BITRES) <= total_bits ?
2364          ec_dec_icdf(dec, trim_icdf, 7) : 5;
2365
2366    bits = ((celt_int32)len*8<<BITRES) - ec_tell_frac(dec) - 1;
2367    anti_collapse_rsv = isTransient&&LM>=2&&bits>=(LM+2<<BITRES) ? (1<<BITRES) : 0;
2368    bits -= anti_collapse_rsv;
2369    codedBands = compute_allocation(st->mode, st->start, st->end, offsets, cap,
2370          alloc_trim, &intensity, &dual_stereo, bits, &balance, pulses,
2371          fine_quant, fine_priority, C, LM, dec, 0, 0);
2372    
2373    unquant_fine_energy(st->mode, st->start, st->end, oldBandE, fine_quant, dec, C);
2374
2375    /* Decode fixed codebook */
2376    ALLOC(collapse_masks, C*st->mode->nbEBands, unsigned char);
2377    quant_all_bands(0, st->mode, st->start, st->end, X, C==2 ? X+N : NULL, collapse_masks,
2378          NULL, pulses, shortBlocks, spread_decision, dual_stereo, intensity, tf_res, 1,
2379          len*(8<<BITRES)-anti_collapse_rsv, balance, dec, LM, codedBands, &st->rng);
2380
2381    if (anti_collapse_rsv > 0)
2382    {
2383       anti_collapse_on = ec_dec_bits(dec, 1);
2384    }
2385
2386    unquant_energy_finalise(st->mode, st->start, st->end, oldBandE,
2387          fine_quant, fine_priority, len*8-ec_tell(dec), dec, C);
2388
2389    if (anti_collapse_on)
2390       anti_collapse(st->mode, X, collapse_masks, LM, C, CC, N,
2391             st->start, st->end, oldBandE, oldLogE, oldLogE2, pulses, st->rng);
2392
2393    log2Amp(st->mode, st->start, st->end, bandE, oldBandE, C);
2394
2395    if (silence)
2396    {
2397       for (i=0;i<C*st->mode->nbEBands;i++)
2398       {
2399          bandE[i] = 0;
2400          oldBandE[i] = -QCONST16(28.f,DB_SHIFT);
2401       }
2402    }
2403    /* Synthesis */
2404    denormalise_bands(st->mode, X, freq, bandE, effEnd, C, M);
2405
2406    CELT_MOVE(decode_mem[0], decode_mem[0]+N, DECODE_BUFFER_SIZE-N);
2407    if (CC==2)
2408       CELT_MOVE(decode_mem[1], decode_mem[1]+N, DECODE_BUFFER_SIZE-N);
2409
2410    c=0; do
2411       for (i=0;i<M*st->mode->eBands[st->start];i++)
2412          freq[c*N+i] = 0;
2413    while (++c<C);
2414    c=0; do {
2415       int bound = M*st->mode->eBands[effEnd];
2416       if (st->downsample!=1)
2417          bound = IMIN(bound, N/st->downsample);
2418       for (i=bound;i<N;i++)
2419          freq[c*N+i] = 0;
2420    } while (++c<C);
2421
2422    out_syn[0] = out_mem[0]+MAX_PERIOD-N;
2423    if (CC==2)
2424       out_syn[1] = out_mem[1]+MAX_PERIOD-N;
2425
2426    if (CC==2&&C==1)
2427    {
2428       for (i=0;i<N;i++)
2429          freq[N+i] = freq[i];
2430    }
2431
2432    /* Compute inverse MDCTs */
2433    compute_inv_mdcts(st->mode, shortBlocks, freq, out_syn, overlap_mem, CC, LM);
2434
2435 #ifdef ENABLE_POSTFILTER
2436    c=0; do {
2437       st->postfilter_period=IMAX(st->postfilter_period, COMBFILTER_MINPERIOD);
2438       st->postfilter_period_old=IMAX(st->postfilter_period_old, COMBFILTER_MINPERIOD);
2439       if (LM!=0)
2440       {
2441          comb_filter(out_syn[c], out_syn[c], st->postfilter_period, st->postfilter_period, st->overlap,
2442                st->postfilter_gain, st->postfilter_gain, st->postfilter_tapset, st->postfilter_tapset,
2443                NULL, 0);
2444          comb_filter(out_syn[c]+st->overlap, out_syn[c]+st->overlap, st->postfilter_period, postfilter_pitch, N-st->overlap,
2445                st->postfilter_gain, postfilter_gain, st->postfilter_tapset, postfilter_tapset,
2446                st->mode->window, st->mode->overlap);
2447       } else {
2448          comb_filter(out_syn[c], out_syn[c], st->postfilter_period_old, st->postfilter_period, N-st->overlap,
2449                st->postfilter_gain_old, st->postfilter_gain, st->postfilter_tapset_old, st->postfilter_tapset,
2450                st->mode->window, st->mode->overlap);
2451       }
2452    } while (++c<CC);
2453    st->postfilter_period_old = st->postfilter_period;
2454    st->postfilter_gain_old = st->postfilter_gain;
2455    st->postfilter_tapset_old = st->postfilter_tapset;
2456    st->postfilter_period = postfilter_pitch;
2457    st->postfilter_gain = postfilter_gain;
2458    st->postfilter_tapset = postfilter_tapset;
2459 #endif /* ENABLE_POSTFILTER */
2460
2461    if (CC==2&&C==1) {
2462       for (i=0;i<st->mode->nbEBands;i++)
2463          oldBandE[st->mode->nbEBands+i]=oldBandE[i];
2464    }
2465
2466    /* In case start or end were to change */
2467    c=0; do
2468    {
2469       for (i=0;i<st->start;i++)
2470          oldBandE[c*st->mode->nbEBands+i]=0;
2471       for (i=st->end;i<st->mode->nbEBands;i++)
2472          oldBandE[c*st->mode->nbEBands+i]=0;
2473    } while (++c<CC);
2474    if (!isTransient)
2475    {
2476       for (i=0;i<CC*st->mode->nbEBands;i++)
2477          oldLogE2[i] = oldLogE[i];
2478       for (i=0;i<CC*st->mode->nbEBands;i++)
2479          oldLogE[i] = oldBandE[i];
2480       for (i=0;i<CC*st->mode->nbEBands;i++)
2481          backgroundLogE[i] = MIN16(backgroundLogE[i] + M*QCONST16(0.001f,DB_SHIFT), oldBandE[i]);
2482    } else {
2483       for (i=0;i<CC*st->mode->nbEBands;i++)
2484          oldLogE[i] = MIN16(oldLogE[i], oldBandE[i]);
2485    }
2486    st->rng = dec->rng;
2487
2488    deemphasis(out_syn, pcm, N, CC, st->downsample, st->mode->preemph, st->preemph_memD);
2489    st->loss_count = 0;
2490    RESTORE_STACK;
2491    if (ec_tell(dec) > 8*len || ec_get_error(dec))
2492       return CELT_CORRUPTED_DATA;
2493    else
2494       return CELT_OK;
2495 }
2496
2497 #ifdef FIXED_POINT
2498 #ifndef DISABLE_FLOAT_API
2499 CELT_STATIC
2500 int celt_decode_with_ec_float(CELTDecoder * restrict st, const unsigned char *data, int len, float * restrict pcm, int frame_size, ec_dec *dec)
2501 {
2502    int j, ret, C, N;
2503    VARDECL(celt_int16, out);
2504    ALLOC_STACK;
2505    SAVE_STACK;
2506
2507    if (pcm==NULL)
2508       return CELT_BAD_ARG;
2509
2510    C = CHANNELS(st->channels);
2511    N = frame_size;
2512    
2513    ALLOC(out, C*N, celt_int16);
2514    ret=celt_decode_with_ec(st, data, len, out, frame_size, dec);
2515    if (ret==0)
2516       for (j=0;j<C*N;j++)
2517          pcm[j]=out[j]*(1.f/32768.f);
2518      
2519    RESTORE_STACK;
2520    return ret;
2521 }
2522 #endif /*DISABLE_FLOAT_API*/
2523 #else
2524 CELT_STATIC
2525 int celt_decode_with_ec(CELTDecoder * restrict st, const unsigned char *data, int len, celt_int16 * restrict pcm, int frame_size, ec_dec *dec)
2526 {
2527    int j, ret, C, N;
2528    VARDECL(celt_sig, out);
2529    ALLOC_STACK;
2530    SAVE_STACK;
2531
2532    if (pcm==NULL)
2533       return CELT_BAD_ARG;
2534
2535    C = CHANNELS(st->channels);
2536    N = frame_size;
2537    ALLOC(out, C*N, celt_sig);
2538
2539    ret=celt_decode_with_ec_float(st, data, len, out, frame_size, dec);
2540
2541    if (ret==0)
2542       for (j=0;j<C*N;j++)
2543          pcm[j] = FLOAT2INT16 (out[j]);
2544    
2545    RESTORE_STACK;
2546    return ret;
2547 }
2548 #endif
2549
2550 int celt_decode(CELTDecoder * restrict st, const unsigned char *data, int len, celt_int16 * restrict pcm, int frame_size)
2551 {
2552    return celt_decode_with_ec(st, data, len, pcm, frame_size, NULL);
2553 }
2554
2555 #ifndef DISABLE_FLOAT_API
2556 int celt_decode_float(CELTDecoder * restrict st, const unsigned char *data, int len, float * restrict pcm, int frame_size)
2557 {
2558    return celt_decode_with_ec_float(st, data, len, pcm, frame_size, NULL);
2559 }
2560 #endif /* DISABLE_FLOAT_API */
2561
2562 int celt_decoder_ctl(CELTDecoder * restrict st, int request, ...)
2563 {
2564    va_list ap;
2565
2566    va_start(ap, request);
2567    switch (request)
2568    {
2569       case CELT_GET_MODE_REQUEST:
2570       {
2571          const CELTMode ** value = va_arg(ap, const CELTMode**);
2572          if (value==0)
2573             goto bad_arg;
2574          *value=st->mode;
2575       }
2576       break;
2577       case CELT_SET_START_BAND_REQUEST:
2578       {
2579          celt_int32 value = va_arg(ap, celt_int32);
2580          if (value<0 || value>=st->mode->nbEBands)
2581             goto bad_arg;
2582          st->start = value;
2583       }
2584       break;
2585       case CELT_SET_END_BAND_REQUEST:
2586       {
2587          celt_int32 value = va_arg(ap, celt_int32);
2588          if (value<1 || value>st->mode->nbEBands)
2589             goto bad_arg;
2590          st->end = value;
2591       }
2592       break;
2593       case CELT_SET_CHANNELS_REQUEST:
2594       {
2595          celt_int32 value = va_arg(ap, celt_int32);
2596          if (value<1 || value>2)
2597             goto bad_arg;
2598          st->stream_channels = value;
2599       }
2600       break;
2601       case CELT_RESET_STATE:
2602       {
2603          CELT_MEMSET((char*)&st->DECODER_RESET_START, 0,
2604                celt_decoder_get_size_custom(st->mode, st->channels)-
2605                ((char*)&st->DECODER_RESET_START - (char*)st));
2606       }
2607       break;
2608       default:
2609          goto bad_request;
2610    }
2611    va_end(ap);
2612    return CELT_OK;
2613 bad_arg:
2614    va_end(ap);
2615    return CELT_BAD_ARG;
2616 bad_request:
2617       va_end(ap);
2618   return CELT_UNIMPLEMENTED;
2619 }
2620
2621 const char *celt_strerror(int error)
2622 {
2623    static const char *error_strings[8] = {
2624       "success",
2625       "invalid argument",
2626       "invalid mode",
2627       "internal error",
2628       "corrupted stream",
2629       "request not implemented",
2630       "invalid state",
2631       "memory allocation failed"
2632    };
2633    if (error > 0 || error < -7)
2634       return "unknown error";
2635    else 
2636       return error_strings[-error];
2637 }
2638