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