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