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