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