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