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