469360400e7dcf6608f803a00bfa0a2c41eae5ca
[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 = 0;
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 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 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    int effectiveBytes;
725    SAVE_STACK;
726
727    if (nbCompressedBytes<0 || pcm==NULL)
728      return CELT_BAD_ARG;
729
730    for (LM=0;LM<4;LM++)
731       if (st->mode->shortMdctSize<<LM==frame_size)
732          break;
733    if (LM>=MAX_CONFIG_SIZES)
734       return CELT_BAD_ARG;
735    M=1<<LM;
736
737    prefilter_mem = st->in_mem+C*(st->overlap);
738    _overlap_mem = prefilter_mem+C*COMBFILTER_MAXPERIOD;
739    /*_overlap_mem = st->in_mem+C*(st->overlap);*/
740    oldBandE = (celt_word16*)(st->in_mem+C*(2*st->overlap+COMBFILTER_MAXPERIOD));
741
742    if (enc==NULL)
743    {
744       ec_byte_writeinit_buffer(&buf, compressed, nbCompressedBytes);
745       ec_enc_init(&_enc,&buf);
746       enc = &_enc;
747       nbFilledBytes=0;
748    } else {
749       nbFilledBytes=(ec_enc_tell(enc, 0)+4)>>3;
750    }
751    nbAvailableBytes = nbCompressedBytes - nbFilledBytes;
752
753    if (st->vbr_rate_norm>0)
754       effectiveBytes = st->vbr_rate_norm>>BITRES<<LM>>3;
755    else
756       effectiveBytes = nbCompressedBytes;
757
758    effEnd = st->end;
759    if (effEnd > st->mode->effEBands)
760       effEnd = st->mode->effEBands;
761
762    N = M*st->mode->shortMdctSize;
763    ALLOC(in, C*(N+st->overlap), celt_sig);
764
765    /* Find pitch period and gain */
766    {
767       VARDECL(celt_sig, _pre);
768       celt_sig *pre[2];
769       SAVE_STACK;
770       c = 0;
771       ALLOC(_pre, C*(N+COMBFILTER_MAXPERIOD), celt_sig);
772
773       pre[0] = _pre;
774       pre[1] = _pre + (N+COMBFILTER_MAXPERIOD);
775
776       c=0; do {
777          const celt_word16 * restrict pcmp = pcm+c;
778          celt_sig * restrict inp = in+c*(N+st->overlap)+st->overlap;
779
780          for (i=0;i<N;i++)
781          {
782             /* Apply pre-emphasis */
783             celt_sig tmp = MULT16_16(st->mode->preemph[2], SCALEIN(*pcmp));
784             *inp = tmp + st->preemph_memE[c];
785             st->preemph_memE[c] = MULT16_32_Q15(st->mode->preemph[1], *inp)
786                                    - MULT16_32_Q15(st->mode->preemph[0], tmp);
787             inp++;
788             pcmp+=C;
789          }
790          CELT_COPY(pre[c], prefilter_mem+c*COMBFILTER_MAXPERIOD, COMBFILTER_MAXPERIOD);
791          CELT_COPY(pre[c]+COMBFILTER_MAXPERIOD, in+c*(N+st->overlap)+st->overlap, N);
792       } while (++c<C);
793
794 #ifdef ENABLE_POSTFILTER
795       {
796          VARDECL(celt_word16, pitch_buf);
797          ALLOC(pitch_buf, (COMBFILTER_MAXPERIOD+N)>>1, celt_word16);
798          celt_word32 tmp=0;
799          celt_word32 mem0[2]={0,0};
800          celt_word16 mem1[2]={0,0};
801
802          pitch_downsample(pre, pitch_buf, COMBFILTER_MAXPERIOD+N, COMBFILTER_MAXPERIOD+N,
803                           C, mem0, mem1);
804          pitch_search(st->mode, pitch_buf+(COMBFILTER_MAXPERIOD>>1), pitch_buf, N,
805                COMBFILTER_MAXPERIOD-COMBFILTER_MINPERIOD, &pitch_index, &tmp, 1<<LM);
806          pitch_index = COMBFILTER_MAXPERIOD-pitch_index;
807
808          gain1 = remove_doubling(pitch_buf, COMBFILTER_MAXPERIOD, COMBFILTER_MINPERIOD,
809                N, &pitch_index, st->prefilter_period, st->prefilter_gain);
810       }
811       if (pitch_index > COMBFILTER_MAXPERIOD)
812          pitch_index = COMBFILTER_MAXPERIOD;
813       gain1 = MULT16_16_Q15(QCONST16(.7f,15),gain1);
814       if (gain1 > QCONST16(.6f,15))
815          gain1 = QCONST16(.6f,15);
816       if (ABS16(gain1-st->prefilter_gain)<QCONST16(.1,15))
817          gain1=st->prefilter_gain;
818       if (gain1<QCONST16(.2f,15))
819       {
820          ec_enc_bit_prob(enc, 0, 32768);
821          gain1 = 0;
822       } else {
823          int qg;
824          int octave;
825 #ifdef FIXED_POINT
826          qg = ((gain1+2048)>>12)-2;
827 #else
828          qg = floor(.5+gain1*8)-2;
829 #endif
830          ec_enc_bit_prob(enc, 1, 32768);
831          octave = EC_ILOG(pitch_index)-5;
832          ec_enc_uint(enc, octave, 6);
833          ec_enc_bits(enc, pitch_index-(16<<octave), 4+octave);
834          ec_enc_bits(enc, qg, 2);
835          gain1 = QCONST16(.125f,15)*(qg+2);
836       }
837       /*printf("%d %f\n", pitch_index, gain1);*/
838 #else /* ENABLE_POSTFILTER */
839       ec_enc_bit_prob(enc, 0, 32768);
840 #endif /* ENABLE_POSTFILTER */
841
842       c=0; do {
843          CELT_COPY(in+c*(N+st->overlap), st->in_mem+c*(st->overlap), st->overlap);
844 #ifdef ENABLE_POSTFILTER
845          comb_filter(in+c*(N+st->overlap)+st->overlap, pre[c]+COMBFILTER_MAXPERIOD,
846                st->prefilter_period, pitch_index, N, C, -st->prefilter_gain, -gain1, st->mode->window, st->mode->overlap);
847 #endif /* ENABLE_POSTFILTER */
848          CELT_COPY(st->in_mem+c*(st->overlap), in+c*(N+st->overlap)+N, st->overlap);
849
850 #ifdef ENABLE_POSTFILTER
851          if (N>COMBFILTER_MAXPERIOD)
852          {
853             CELT_MOVE(prefilter_mem+c*COMBFILTER_MAXPERIOD, pre[c]+N, COMBFILTER_MAXPERIOD);
854          } else {
855             CELT_MOVE(prefilter_mem+c*COMBFILTER_MAXPERIOD, prefilter_mem+c*COMBFILTER_MAXPERIOD+N, COMBFILTER_MAXPERIOD-N);
856             CELT_MOVE(prefilter_mem+c*COMBFILTER_MAXPERIOD+COMBFILTER_MAXPERIOD-N, pre[c]+COMBFILTER_MAXPERIOD, N);
857          }
858 #endif /* ENABLE_POSTFILTER */
859       } while (++c<C);
860
861       RESTORE_STACK;
862    }
863
864 #ifdef RESYNTH
865    resynth = 1;
866 #else
867    resynth = 0;
868 #endif
869
870    if (st->complexity > 1 && LM>0)
871    {
872       isTransient = M > 1 &&
873          transient_analysis(in, N+st->overlap, C, &st->frame_max, st->overlap);
874    } else {
875       isTransient = 0;
876    }
877
878    if (isTransient)
879       shortBlocks = M;
880    else
881       shortBlocks = 0;
882
883    ALLOC(freq, C*N, celt_sig); /**< Interleaved signal MDCTs */
884    ALLOC(bandE,st->mode->nbEBands*C, celt_ener);
885    ALLOC(bandLogE,st->mode->nbEBands*C, celt_word16);
886    /* Compute MDCTs */
887    compute_mdcts(st->mode, shortBlocks, in, freq, C, LM);
888
889    ALLOC(X, C*N, celt_norm);         /**< Interleaved normalised MDCTs */
890
891    compute_band_energies(st->mode, freq, bandE, effEnd, C, M);
892
893    amp2Log2(st->mode, effEnd, st->end, bandE, bandLogE, C);
894
895    /* Band normalisation */
896    normalise_bands(st->mode, freq, X, bandE, effEnd, C, M);
897
898    ALLOC(tf_res, st->mode->nbEBands, int);
899    /* Needs to be before coarse energy quantization because otherwise the energy gets modified */
900    tf_select = tf_analysis(st->mode, bandLogE, oldBandE, effEnd, C, isTransient, tf_res, effectiveBytes, X, N, LM, &tf_sum);
901    for (i=effEnd;i<st->end;i++)
902       tf_res[i] = tf_res[effEnd-1];
903
904    ALLOC(error, C*st->mode->nbEBands, celt_word16);
905    quant_coarse_energy(st->mode, st->start, st->end, effEnd, bandLogE,
906          oldBandE, nbCompressedBytes*8, error, enc,
907          C, LM, nbAvailableBytes, st->force_intra,
908          &st->delayedIntra, st->complexity >= 4);
909
910    if (LM > 0)
911       ec_enc_bit_prob(enc, shortBlocks!=0, 8192);
912
913    tf_encode(st->start, st->end, isTransient, tf_res, LM, tf_select, enc);
914
915    if (shortBlocks || st->complexity < 3)
916    {
917       if (st->complexity == 0)
918       {
919          has_fold = 0;
920          st->fold_decision = 3;
921       } else {
922          has_fold = 1;
923          st->fold_decision = 1;
924       }
925    } else {
926       has_fold = folding_decision(st->mode, X, &st->tonal_average, &st->fold_decision, effEnd, C, M);
927    }
928    ec_enc_bit_prob(enc, has_fold>>1, 8192);
929    ec_enc_bit_prob(enc, has_fold&1, (has_fold>>1) ? 32768 : 49152);
930
931    ALLOC(offsets, st->mode->nbEBands, int);
932
933    for (i=0;i<st->mode->nbEBands;i++)
934       offsets[i] = 0;
935    /* Dynamic allocation code */
936    /* Make sure that dynamic allocation can't make us bust the budget */
937    if (effectiveBytes > 50 && LM>=1)
938    {
939       int t1, t2;
940       if (LM <= 1)
941       {
942          t1 = 3;
943          t2 = 5;
944       } else {
945          t1 = 2;
946          t2 = 4;
947       }
948       for (i=1;i<st->mode->nbEBands-1;i++)
949       {
950          celt_word32 d2;
951          d2 = 2*bandLogE[i]-bandLogE[i-1]-bandLogE[i+1];
952          if (C==2)
953             d2 = HALF32(d2 + 2*bandLogE[i+st->mode->nbEBands]-
954                   bandLogE[i-1+st->mode->nbEBands]-bandLogE[i+1+st->mode->nbEBands]);
955          if (d2 > SHL16(t1,DB_SHIFT))
956             offsets[i] += 1;
957          if (d2 > SHL16(t2,DB_SHIFT))
958             offsets[i] += 1;
959       }
960    }
961    for (i=0;i<st->mode->nbEBands;i++)
962    {
963       int j;
964       ec_enc_bit_prob(enc, offsets[i]!=0, 1024);
965       if (offsets[i]!=0)
966       {
967          for (j=0;j<offsets[i]-1;j++)
968             ec_enc_bit_prob(enc, 1, 32768);
969          ec_enc_bit_prob(enc, 0, 32768);
970       }
971       offsets[i] *= (6<<BITRES);
972    }
973    alloc_trim = alloc_trim_analysis(st->mode, X, bandLogE, st->mode->nbEBands, LM, C, N);
974    ec_encode_bin(enc, trim_cdf[alloc_trim], trim_cdf[alloc_trim+1], 7);
975
976    /* Variable bitrate */
977    if (st->vbr_rate_norm>0)
978    {
979      celt_word16 alpha;
980      celt_int32 delta, tell;
981      /* The target rate in 8th bits per frame */
982      celt_int32 vbr_rate;
983      celt_int32 target;
984      celt_int32 vbr_bound, max_allowed, min_allowed;
985
986      target = vbr_rate = M*st->vbr_rate_norm;
987
988      target = target + st->vbr_offset - ((40*C+20)<<BITRES);
989
990      /* Shortblocks get a large boost in bitrate, but since they
991         are uncommon long blocks are not greatly affected */
992      if (shortBlocks || tf_sum < -2*(st->end-st->start))
993         target = 7*target/4;
994      else if (tf_sum < -(st->end-st->start))
995         target = 3*target/2;
996      else if (M > 1)
997         target-=(target+14)/28;
998
999      tell = ec_enc_tell(enc, BITRES);
1000
1001      /* The current offset is removed from the target and the space used
1002         so far is added*/
1003      target=target+tell;
1004      /* By how much did we "miss" the target on that frame */
1005      delta = target - vbr_rate;
1006
1007      /* Computes the max bit-rate allowed in VBR more to avoid violating the target rate and buffering */
1008      vbr_bound = vbr_rate;
1009      max_allowed = IMIN(vbr_rate+vbr_bound-st->vbr_reservoir>>(BITRES+3),nbAvailableBytes);
1010      min_allowed = (tell>>(BITRES+3)) + 2 - nbFilledBytes;
1011
1012      /* In VBR mode the frame size must not be reduced so much that it would result in the encoder running out of bits */
1013      nbAvailableBytes = target+(1<<(BITRES+2))>>(BITRES+3);
1014      nbAvailableBytes=IMAX(min_allowed,IMIN(max_allowed,nbAvailableBytes));
1015      target=nbAvailableBytes<<(BITRES+3);
1016
1017      if (st->vbr_count < 970)
1018      {
1019         st->vbr_count++;
1020         alpha = celt_rcp(SHL32(EXTEND32(st->vbr_count+20),16));
1021      } else
1022         alpha = QCONST16(.001f,15);
1023      /* How many bits have we used in excess of what we're allowed */
1024      st->vbr_reservoir += target - vbr_rate;
1025      /*printf ("%d\n", st->vbr_reservoir);*/
1026
1027      /* Compute the offset we need to apply in order to reach the target */
1028      st->vbr_drift += (celt_int32)MULT16_32_Q15(alpha,delta-st->vbr_offset-st->vbr_drift);
1029      st->vbr_offset = -st->vbr_drift;
1030      /*printf ("%d\n", st->vbr_drift);*/
1031
1032      /* We could use any multiple of vbr_rate as bound (depending on the delay) */
1033      if (st->vbr_reservoir < 0)
1034      {
1035         /* We're under the min value -- increase rate */
1036         int adjust = (-st->vbr_reservoir)/(8<<BITRES);
1037         nbAvailableBytes += adjust;
1038         st->vbr_reservoir = 0;
1039         /*printf ("+%d\n", adjust);*/
1040      }
1041      nbCompressedBytes = IMIN(nbCompressedBytes,nbAvailableBytes+nbFilledBytes);
1042
1043      /* This moves the raw bits to take into account the new compressed size */
1044      ec_byte_shrink(&buf, nbCompressedBytes);
1045    }
1046
1047    if (C==2)
1048    {
1049       dual_stereo = stereo_analysis(st->mode, X, st->mode->nbEBands, LM, C, N);
1050       ec_enc_bit_prob(enc, dual_stereo, 32768);
1051    }
1052    if (C==2)
1053    {
1054       int effectiveRate;
1055
1056       /* Account for coarse energy */
1057       effectiveRate = (8*effectiveBytes - 80)>>LM;
1058
1059       /* effectiveRate in kb/s */
1060       effectiveRate = 2*effectiveRate/5;
1061       if (effectiveRate<35)
1062          intensity = 6;
1063       else if (effectiveRate<50)
1064          intensity = 12;
1065       else if (effectiveRate<68)
1066          intensity = 16;
1067       else if (effectiveRate<84)
1068          intensity = 18;
1069       else if (effectiveRate<102)
1070          intensity = 19;
1071       else if (effectiveRate<130)
1072          intensity = 20;
1073       else
1074          intensity = 100;
1075       intensity = IMIN(st->end,IMAX(st->start, intensity));
1076       ec_enc_uint(enc, intensity, 1+st->end-st->start);
1077    }
1078
1079    /* Bit allocation */
1080    ALLOC(fine_quant, st->mode->nbEBands, int);
1081    ALLOC(pulses, st->mode->nbEBands, int);
1082    ALLOC(fine_priority, st->mode->nbEBands, int);
1083
1084    bits = nbCompressedBytes*8 - ec_enc_tell(enc, 0) - 1;
1085    codedBands = compute_allocation(st->mode, st->start, st->end, offsets, alloc_trim, bits, pulses, fine_quant, fine_priority, C, LM);
1086
1087    quant_fine_energy(st->mode, st->start, st->end, bandE, oldBandE, error, fine_quant, enc, C);
1088
1089 #ifdef MEASURE_NORM_MSE
1090    float X0[3000];
1091    float bandE0[60];
1092    c=0; do 
1093       for (i=0;i<N;i++)
1094          X0[i+c*N] = X[i+c*N];
1095    while (++c<C);
1096    for (i=0;i<C*st->mode->nbEBands;i++)
1097       bandE0[i] = bandE[i];
1098 #endif
1099
1100    /* Residual quantisation */
1101    quant_all_bands(1, st->mode, st->start, st->end, X, C==2 ? X+N : NULL,
1102          bandE, pulses, shortBlocks, has_fold, dual_stereo, intensity, tf_res, resynth,
1103          nbCompressedBytes*8, enc, LM, codedBands);
1104
1105    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);
1106
1107 #ifdef RESYNTH
1108    /* Re-synthesis of the coded audio if required */
1109    if (resynth)
1110    {
1111       celt_sig *out_mem[2];
1112       celt_sig *overlap_mem[2];
1113
1114       log2Amp(st->mode, st->start, st->end, bandE, oldBandE, C);
1115
1116 #ifdef MEASURE_NORM_MSE
1117       measure_norm_mse(st->mode, X, X0, bandE, bandE0, M, N, C);
1118 #endif
1119
1120       /* Synthesis */
1121       denormalise_bands(st->mode, X, freq, bandE, effEnd, C, M);
1122
1123       CELT_MOVE(st->syn_mem[0], st->syn_mem[0]+N, MAX_PERIOD);
1124       if (C==2)
1125          CELT_MOVE(st->syn_mem[1], st->syn_mem[1]+N, MAX_PERIOD);
1126
1127       c=0; do
1128          for (i=0;i<M*st->mode->eBands[st->start];i++)
1129             freq[c*N+i] = 0;
1130       while (++c<C);
1131       c=0; do
1132          for (i=M*st->mode->eBands[st->end];i<N;i++)
1133             freq[c*N+i] = 0;
1134       while (++c<C);
1135
1136       out_mem[0] = st->syn_mem[0]+MAX_PERIOD;
1137       if (C==2)
1138          out_mem[1] = st->syn_mem[1]+MAX_PERIOD;
1139
1140       c=0; do
1141          overlap_mem[c] = _overlap_mem + c*st->overlap;
1142       while (++c<C);
1143
1144       compute_inv_mdcts(st->mode, shortBlocks, freq, out_mem, overlap_mem, C, LM);
1145
1146 #ifdef ENABLE_POSTFILTER
1147       c=0; do {
1148          comb_filter(out_mem[c], out_mem[c], st->prefilter_period, st->prefilter_period, st->overlap, C,
1149                st->prefilter_gain, st->prefilter_gain, NULL, 0);
1150          comb_filter(out_mem[c]+st->overlap, out_mem[c]+st->overlap, st->prefilter_period, pitch_index, N-st->overlap, C,
1151                st->prefilter_gain, gain1, st->mode->window, st->mode->overlap);
1152       } while (++c<C);
1153 #endif /* ENABLE_POSTFILTER */
1154
1155       deemphasis(out_mem, (celt_word16*)pcm, N, C, st->mode->preemph, st->preemph_memD);
1156    }
1157 #endif
1158
1159    st->prefilter_period = pitch_index;
1160    st->prefilter_gain = gain1;
1161
1162    /* If there's any room left (can only happen for very high rates),
1163       fill it with zeros */
1164    while (ec_enc_tell(enc,0) + 8 <= nbCompressedBytes*8)
1165       ec_enc_bits(enc, 0, 8);
1166    ec_enc_done(enc);
1167    
1168    RESTORE_STACK;
1169    if (ec_enc_get_error(enc))
1170       return CELT_CORRUPTED_DATA;
1171    else
1172       return nbCompressedBytes;
1173 }
1174
1175 #ifdef FIXED_POINT
1176 #ifndef DISABLE_FLOAT_API
1177 int celt_encode_with_ec_float(CELTEncoder * restrict st, const float * pcm, int frame_size, unsigned char *compressed, int nbCompressedBytes, ec_enc *enc)
1178 {
1179    int j, ret, C, N, LM, M;
1180    VARDECL(celt_int16, in);
1181    SAVE_STACK;
1182
1183    if (pcm==NULL)
1184       return CELT_BAD_ARG;
1185
1186    for (LM=0;LM<4;LM++)
1187       if (st->mode->shortMdctSize<<LM==frame_size)
1188          break;
1189    if (LM>=MAX_CONFIG_SIZES)
1190       return CELT_BAD_ARG;
1191    M=1<<LM;
1192
1193    C = CHANNELS(st->channels);
1194    N = M*st->mode->shortMdctSize;
1195    ALLOC(in, C*N, celt_int16);
1196
1197    for (j=0;j<C*N;j++)
1198      in[j] = FLOAT2INT16(pcm[j]);
1199
1200    ret=celt_encode_with_ec(st,in,frame_size,compressed,nbCompressedBytes, enc);
1201 #ifdef RESYNTH
1202    for (j=0;j<C*N;j++)
1203       ((float*)pcm)[j]=in[j]*(1.f/32768.f);
1204 #endif
1205    RESTORE_STACK;
1206    return ret;
1207
1208 }
1209 #endif /*DISABLE_FLOAT_API*/
1210 #else
1211 int celt_encode_with_ec(CELTEncoder * restrict st, const celt_int16 * pcm, int frame_size, unsigned char *compressed, int nbCompressedBytes, ec_enc *enc)
1212 {
1213    int j, ret, C, N, LM, M;
1214    VARDECL(celt_sig, in);
1215    SAVE_STACK;
1216
1217    if (pcm==NULL)
1218       return CELT_BAD_ARG;
1219
1220    for (LM=0;LM<4;LM++)
1221       if (st->mode->shortMdctSize<<LM==frame_size)
1222          break;
1223    if (LM>=MAX_CONFIG_SIZES)
1224       return CELT_BAD_ARG;
1225    M=1<<LM;
1226
1227    C=CHANNELS(st->channels);
1228    N=M*st->mode->shortMdctSize;
1229    ALLOC(in, C*N, celt_sig);
1230    for (j=0;j<C*N;j++) {
1231      in[j] = SCALEOUT(pcm[j]);
1232    }
1233
1234    ret = celt_encode_with_ec_float(st,in,frame_size,compressed,nbCompressedBytes, enc);
1235 #ifdef RESYNTH
1236    for (j=0;j<C*N;j++)
1237       ((celt_int16*)pcm)[j] = FLOAT2INT16(in[j]);
1238 #endif
1239    RESTORE_STACK;
1240    return ret;
1241 }
1242 #endif
1243
1244 int celt_encode(CELTEncoder * restrict st, const celt_int16 * pcm, int frame_size, unsigned char *compressed, int nbCompressedBytes)
1245 {
1246    return celt_encode_with_ec(st, pcm, frame_size, compressed, nbCompressedBytes, NULL);
1247 }
1248
1249 #ifndef DISABLE_FLOAT_API
1250 int celt_encode_float(CELTEncoder * restrict st, const float * pcm, int frame_size, unsigned char *compressed, int nbCompressedBytes)
1251 {
1252    return celt_encode_with_ec_float(st, pcm, frame_size, compressed, nbCompressedBytes, NULL);
1253 }
1254 #endif /* DISABLE_FLOAT_API */
1255
1256 int celt_encoder_ctl(CELTEncoder * restrict st, int request, ...)
1257 {
1258    va_list ap;
1259    
1260    va_start(ap, request);
1261    switch (request)
1262    {
1263       case CELT_GET_MODE_REQUEST:
1264       {
1265          const CELTMode ** value = va_arg(ap, const CELTMode**);
1266          if (value==0)
1267             goto bad_arg;
1268          *value=st->mode;
1269       }
1270       break;
1271       case CELT_SET_COMPLEXITY_REQUEST:
1272       {
1273          int value = va_arg(ap, celt_int32);
1274          if (value<0 || value>10)
1275             goto bad_arg;
1276          st->complexity = value;
1277       }
1278       break;
1279       case CELT_SET_START_BAND_REQUEST:
1280       {
1281          celt_int32 value = va_arg(ap, celt_int32);
1282          if (value<0 || value>=st->mode->nbEBands)
1283             goto bad_arg;
1284          st->start = value;
1285       }
1286       break;
1287       case CELT_SET_END_BAND_REQUEST:
1288       {
1289          celt_int32 value = va_arg(ap, celt_int32);
1290          if (value<0 || value>=st->mode->nbEBands)
1291             goto bad_arg;
1292          st->end = value;
1293       }
1294       break;
1295       case CELT_SET_PREDICTION_REQUEST:
1296       {
1297          int value = va_arg(ap, celt_int32);
1298          if (value<0 || value>2)
1299             goto bad_arg;
1300          if (value==0)
1301          {
1302             st->force_intra   = 1;
1303          } else if (value==1) {
1304             st->force_intra   = 0;
1305          } else {
1306             st->force_intra   = 0;
1307          }   
1308       }
1309       break;
1310       case CELT_SET_VBR_RATE_REQUEST:
1311       {
1312          celt_int32 value = va_arg(ap, celt_int32);
1313          int frame_rate;
1314          int N = st->mode->shortMdctSize;
1315          if (value<0)
1316             goto bad_arg;
1317          if (value>3072000)
1318             value = 3072000;
1319          frame_rate = ((st->mode->Fs<<3)+(N>>1))/N;
1320          st->vbr_rate_norm = ((value<<(BITRES+3))+(frame_rate>>1))/frame_rate;
1321       }
1322       break;
1323       case CELT_RESET_STATE:
1324       {
1325          CELT_MEMSET((char*)&st->ENCODER_RESET_START, 0,
1326                celt_encoder_get_size(st->mode, st->channels)-
1327                ((char*)&st->ENCODER_RESET_START - (char*)st));
1328          st->vbr_offset = 0;
1329          st->delayedIntra = 1;
1330          st->fold_decision = 1;
1331          st->tonal_average = QCONST16(1.f,8);
1332       }
1333       break;
1334       default:
1335          goto bad_request;
1336    }
1337    va_end(ap);
1338    return CELT_OK;
1339 bad_arg:
1340    va_end(ap);
1341    return CELT_BAD_ARG;
1342 bad_request:
1343    va_end(ap);
1344    return CELT_UNIMPLEMENTED;
1345 }
1346
1347 /**********************************************************************/
1348 /*                                                                    */
1349 /*                             DECODER                                */
1350 /*                                                                    */
1351 /**********************************************************************/
1352 #define DECODE_BUFFER_SIZE 2048
1353
1354 /** Decoder state 
1355  @brief Decoder state
1356  */
1357 struct CELTDecoder {
1358    const CELTMode *mode;
1359    int overlap;
1360    int channels;
1361
1362    int start, end;
1363
1364    /* Everything beyond this point gets cleared on a reset */
1365 #define DECODER_RESET_START last_pitch_index
1366
1367    int last_pitch_index;
1368    int loss_count;
1369    int postfilter_period;
1370    celt_word16 postfilter_gain;
1371
1372    celt_sig preemph_memD[2];
1373    
1374    celt_sig _decode_mem[1]; /* Size = channels*(DECODE_BUFFER_SIZE+mode->overlap) */
1375    /* celt_word16 lpc[],  Size = channels*LPC_ORDER */
1376    /* celt_word16 oldEBands[], Size = channels*mode->nbEBands */
1377 };
1378
1379 int celt_decoder_get_size(const CELTMode *mode, int channels)
1380 {
1381    int size = sizeof(struct CELTDecoder)
1382             + (channels*(DECODE_BUFFER_SIZE+mode->overlap)-1)*sizeof(celt_sig)
1383             + channels*LPC_ORDER*sizeof(celt_word16)
1384             + channels*mode->nbEBands*sizeof(celt_word16);
1385    return size;
1386 }
1387
1388 CELTDecoder *celt_decoder_create(const CELTMode *mode, int channels, int *error)
1389 {
1390    return celt_decoder_init(
1391          (CELTDecoder *)celt_alloc(celt_decoder_get_size(mode, channels)),
1392          mode, channels, error);
1393 }
1394
1395 CELTDecoder *celt_decoder_init(CELTDecoder *st, const CELTMode *mode, int channels, int *error)
1396 {
1397    if (channels < 0 || channels > 2)
1398    {
1399       if (error)
1400          *error = CELT_BAD_ARG;
1401       return NULL;
1402    }
1403
1404    if (st==NULL)
1405    {
1406       if (error)
1407          *error = CELT_ALLOC_FAIL;
1408       return NULL;
1409    }
1410
1411    CELT_MEMSET((char*)st, 0, celt_decoder_get_size(mode, channels));
1412
1413    st->mode = mode;
1414    st->overlap = mode->overlap;
1415    st->channels = channels;
1416
1417    st->start = 0;
1418    st->end = st->mode->effEBands;
1419
1420    st->loss_count = 0;
1421
1422    if (error)
1423       *error = CELT_OK;
1424    return st;
1425 }
1426
1427 void celt_decoder_destroy(CELTDecoder *st)
1428 {
1429    celt_free(st);
1430 }
1431
1432 static void celt_decode_lost(CELTDecoder * restrict st, celt_word16 * restrict pcm, int N, int LM)
1433 {
1434    int c;
1435    int pitch_index;
1436    int overlap = st->mode->overlap;
1437    celt_word16 fade = Q15ONE;
1438    int i, len;
1439    const int C = CHANNELS(st->channels);
1440    int offset;
1441    celt_sig *out_mem[2];
1442    celt_sig *decode_mem[2];
1443    celt_sig *overlap_mem[2];
1444    celt_word16 *lpc;
1445    SAVE_STACK;
1446    
1447    c=0; do {
1448       decode_mem[c] = st->_decode_mem + c*(DECODE_BUFFER_SIZE+st->overlap);
1449       out_mem[c] = decode_mem[c]+DECODE_BUFFER_SIZE-MAX_PERIOD;
1450       overlap_mem[c] = decode_mem[c]+DECODE_BUFFER_SIZE;
1451    } while (++c<C);
1452    lpc = (celt_word16*)(st->_decode_mem+(DECODE_BUFFER_SIZE+st->overlap)*C);
1453
1454    len = N+st->mode->overlap;
1455    
1456    if (st->loss_count == 0)
1457    {
1458       celt_word16 pitch_buf[MAX_PERIOD>>1];
1459       celt_word32 tmp=0;
1460       celt_word32 mem0[2]={0,0};
1461       celt_word16 mem1[2]={0,0};
1462       int len2 = len;
1463       /* FIXME: This is a kludge */
1464       if (len2>MAX_PERIOD>>1)
1465          len2 = MAX_PERIOD>>1;
1466       pitch_downsample(out_mem, pitch_buf, MAX_PERIOD, MAX_PERIOD,
1467                        C, mem0, mem1);
1468       pitch_search(st->mode, pitch_buf+((MAX_PERIOD-len2)>>1), pitch_buf, len2,
1469                    MAX_PERIOD-len2-100, &pitch_index, &tmp, 1<<LM);
1470       pitch_index = MAX_PERIOD-len2-pitch_index;
1471       st->last_pitch_index = pitch_index;
1472    } else {
1473       pitch_index = st->last_pitch_index;
1474       if (st->loss_count < 5)
1475          fade = QCONST16(.8f,15);
1476       else
1477          fade = 0;
1478    }
1479
1480    c=0; do {
1481       /* FIXME: This is more memory than necessary */
1482       celt_word32 e[2*MAX_PERIOD];
1483       celt_word16 exc[2*MAX_PERIOD];
1484       celt_word32 ac[LPC_ORDER+1];
1485       celt_word16 decay = 1;
1486       celt_word32 S1=0;
1487       celt_word16 mem[LPC_ORDER]={0};
1488
1489       offset = MAX_PERIOD-pitch_index;
1490       for (i=0;i<MAX_PERIOD;i++)
1491          exc[i] = ROUND16(out_mem[c][i], SIG_SHIFT);
1492
1493       if (st->loss_count == 0)
1494       {
1495          _celt_autocorr(exc, ac, st->mode->window, st->mode->overlap,
1496                         LPC_ORDER, MAX_PERIOD);
1497
1498          /* Noise floor -40 dB */
1499 #ifdef FIXED_POINT
1500          ac[0] += SHR32(ac[0],13);
1501 #else
1502          ac[0] *= 1.0001f;
1503 #endif
1504          /* Lag windowing */
1505          for (i=1;i<=LPC_ORDER;i++)
1506          {
1507             /*ac[i] *= exp(-.5*(2*M_PI*.002*i)*(2*M_PI*.002*i));*/
1508 #ifdef FIXED_POINT
1509             ac[i] -= MULT16_32_Q15(2*i*i, ac[i]);
1510 #else
1511             ac[i] -= ac[i]*(.008f*i)*(.008f*i);
1512 #endif
1513          }
1514
1515          _celt_lpc(lpc+c*LPC_ORDER, ac, LPC_ORDER);
1516       }
1517       for (i=0;i<LPC_ORDER;i++)
1518          mem[i] = ROUND16(out_mem[c][MAX_PERIOD-1-i], SIG_SHIFT);
1519       fir(exc, lpc+c*LPC_ORDER, exc, MAX_PERIOD, LPC_ORDER, mem);
1520       /*for (i=0;i<MAX_PERIOD;i++)printf("%d ", exc[i]); printf("\n");*/
1521       /* Check if the waveform is decaying (and if so how fast) */
1522       {
1523          celt_word32 E1=1, E2=1;
1524          int period;
1525          if (pitch_index <= MAX_PERIOD/2)
1526             period = pitch_index;
1527          else
1528             period = MAX_PERIOD/2;
1529          for (i=0;i<period;i++)
1530          {
1531             E1 += SHR32(MULT16_16(exc[MAX_PERIOD-period+i],exc[MAX_PERIOD-period+i]),8);
1532             E2 += SHR32(MULT16_16(exc[MAX_PERIOD-2*period+i],exc[MAX_PERIOD-2*period+i]),8);
1533          }
1534          if (E1 > E2)
1535             E1 = E2;
1536          decay = celt_sqrt(frac_div32(SHR(E1,1),E2));
1537       }
1538
1539       /* Copy excitation, taking decay into account */
1540       for (i=0;i<len+st->mode->overlap;i++)
1541       {
1542          celt_word16 tmp;
1543          if (offset+i >= MAX_PERIOD)
1544          {
1545             offset -= pitch_index;
1546             decay = MULT16_16_Q15(decay, decay);
1547          }
1548          e[i] = SHL32(EXTEND32(MULT16_16_Q15(decay, exc[offset+i])), SIG_SHIFT);
1549          tmp = ROUND16(out_mem[c][offset+i],SIG_SHIFT);
1550          S1 += SHR32(MULT16_16(tmp,tmp),8);
1551       }
1552       for (i=0;i<LPC_ORDER;i++)
1553          mem[i] = ROUND16(out_mem[c][MAX_PERIOD-1-i], SIG_SHIFT);
1554       for (i=0;i<len+st->mode->overlap;i++)
1555          e[i] = MULT16_32_Q15(fade, e[i]);
1556       iir(e, lpc+c*LPC_ORDER, e, len+st->mode->overlap, LPC_ORDER, mem);
1557
1558       {
1559          celt_word32 S2=0;
1560          for (i=0;i<len+overlap;i++)
1561          {
1562             celt_word16 tmp = ROUND16(e[i],SIG_SHIFT);
1563             S2 += SHR32(MULT16_16(tmp,tmp),8);
1564          }
1565          /* This checks for an "explosion" in the synthesis */
1566 #ifdef FIXED_POINT
1567          if (!(S1 > SHR32(S2,2)))
1568 #else
1569          /* Float test is written this way to catch NaNs at the same time */
1570          if (!(S1 > 0.2f*S2))
1571 #endif
1572          {
1573             for (i=0;i<len+overlap;i++)
1574                e[i] = 0;
1575          } else if (S1 < S2)
1576          {
1577             celt_word16 ratio = celt_sqrt(frac_div32(SHR32(S1,1)+1,S2+1));
1578             for (i=0;i<len+overlap;i++)
1579                e[i] = MULT16_32_Q15(ratio, e[i]);
1580          }
1581       }
1582
1583 #ifdef ENABLE_POSTFILTER
1584       /* Apply post-filter to the MDCT overlap of the previous frame */
1585       comb_filter(out_mem[c]+MAX_PERIOD, out_mem[c]+MAX_PERIOD, st->postfilter_period, st->postfilter_period, st->overlap, C,
1586                   st->postfilter_gain, st->postfilter_gain, NULL, 0);
1587 #endif /* ENABLE_POSTFILTER */
1588
1589       for (i=0;i<MAX_PERIOD+st->mode->overlap-N;i++)
1590          out_mem[c][i] = out_mem[c][N+i];
1591
1592       /* Apply TDAC to the concealed audio so that it blends with the
1593          previous and next frames */
1594       for (i=0;i<overlap/2;i++)
1595       {
1596          celt_word32 tmp;
1597          tmp = MULT16_32_Q15(st->mode->window[i],           e[N+overlap-1-i]) +
1598                MULT16_32_Q15(st->mode->window[overlap-i-1], e[N+i          ]);
1599          out_mem[c][MAX_PERIOD+i] = MULT16_32_Q15(st->mode->window[overlap-i-1], tmp);
1600          out_mem[c][MAX_PERIOD+overlap-i-1] = MULT16_32_Q15(st->mode->window[i], tmp);
1601       }
1602       for (i=0;i<N;i++)
1603          out_mem[c][MAX_PERIOD-N+i] = e[i];
1604
1605 #ifdef ENABLE_POSTFILTER
1606       /* Apply pre-filter to the MDCT overlap for the next frame (post-filter will be applied then) */
1607       comb_filter(e, out_mem[c]+MAX_PERIOD, st->postfilter_period, st->postfilter_period, st->overlap, C,
1608                   -st->postfilter_gain, -st->postfilter_gain, NULL, 0);
1609 #endif /* ENABLE_POSTFILTER */
1610       for (i=0;i<overlap;i++)
1611          out_mem[c][MAX_PERIOD+i] = e[i];
1612    } while (++c<C);
1613
1614    {
1615       celt_word32 *out_syn[2];
1616       out_syn[0] = out_mem[0]+MAX_PERIOD-N;
1617       if (C==2)
1618          out_syn[1] = out_mem[1]+MAX_PERIOD-N;
1619       deemphasis(out_syn, pcm, N, C, st->mode->preemph, st->preemph_memD);
1620    }
1621    
1622    st->loss_count++;
1623
1624    RESTORE_STACK;
1625 }
1626
1627 #ifdef FIXED_POINT
1628 int celt_decode_with_ec(CELTDecoder * restrict st, const unsigned char *data, int len, celt_int16 * restrict pcm, int frame_size, ec_dec *dec)
1629 {
1630 #else
1631 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)
1632 {
1633 #endif
1634    int c, i, N;
1635    int has_fold;
1636    int bits;
1637    ec_dec _dec;
1638    ec_byte_buffer buf;
1639    VARDECL(celt_sig, freq);
1640    VARDECL(celt_norm, X);
1641    VARDECL(celt_ener, bandE);
1642    VARDECL(int, fine_quant);
1643    VARDECL(int, pulses);
1644    VARDECL(int, offsets);
1645    VARDECL(int, fine_priority);
1646    VARDECL(int, tf_res);
1647    celt_sig *out_mem[2];
1648    celt_sig *decode_mem[2];
1649    celt_sig *overlap_mem[2];
1650    celt_sig *out_syn[2];
1651    celt_word16 *lpc;
1652    celt_word16 *oldBandE;
1653
1654    int shortBlocks;
1655    int isTransient;
1656    int intra_ener;
1657    const int C = CHANNELS(st->channels);
1658    int LM, M;
1659    int nbFilledBytes, nbAvailableBytes;
1660    int effEnd;
1661    int codedBands;
1662    int alloc_trim;
1663    int postfilter_pitch;
1664    celt_word16 postfilter_gain;
1665    int intensity=0;
1666    int dual_stereo=0;
1667    SAVE_STACK;
1668
1669    if (pcm==NULL)
1670       return CELT_BAD_ARG;
1671
1672    for (LM=0;LM<4;LM++)
1673       if (st->mode->shortMdctSize<<LM==frame_size)
1674          break;
1675    if (LM>=MAX_CONFIG_SIZES)
1676       return CELT_BAD_ARG;
1677    M=1<<LM;
1678
1679    c=0; do {
1680       decode_mem[c] = st->_decode_mem + c*(DECODE_BUFFER_SIZE+st->overlap);
1681       out_mem[c] = decode_mem[c]+DECODE_BUFFER_SIZE-MAX_PERIOD;
1682       overlap_mem[c] = decode_mem[c]+DECODE_BUFFER_SIZE;
1683    } while (++c<C);
1684    lpc = (celt_word16*)(st->_decode_mem+(DECODE_BUFFER_SIZE+st->overlap)*C);
1685    oldBandE = lpc+C*LPC_ORDER;
1686
1687    N = M*st->mode->shortMdctSize;
1688
1689    effEnd = st->end;
1690    if (effEnd > st->mode->effEBands)
1691       effEnd = st->mode->effEBands;
1692
1693    ALLOC(freq, C*N, celt_sig); /**< Interleaved signal MDCTs */
1694    ALLOC(X, C*N, celt_norm);   /**< Interleaved normalised MDCTs */
1695    ALLOC(bandE, st->mode->nbEBands*C, celt_ener);
1696    c=0; do
1697       for (i=0;i<M*st->mode->eBands[st->start];i++)
1698          X[c*N+i] = 0;
1699    while (++c<C);
1700    c=0; do   
1701       for (i=M*st->mode->eBands[effEnd];i<N;i++)
1702          X[c*N+i] = 0;
1703    while (++c<C);
1704
1705    if (data == NULL)
1706    {
1707       celt_decode_lost(st, pcm, N, LM);
1708       RESTORE_STACK;
1709       return CELT_OK;
1710    }
1711    if (len<0) {
1712      RESTORE_STACK;
1713      return CELT_BAD_ARG;
1714    }
1715    
1716    if (dec == NULL)
1717    {
1718       ec_byte_readinit(&buf,(unsigned char*)data,len);
1719       ec_dec_init(&_dec,&buf);
1720       dec = &_dec;
1721       nbFilledBytes = 0;
1722    } else {
1723       nbFilledBytes = (ec_dec_tell(dec, 0)+4)>>3;
1724    }
1725    nbAvailableBytes = len-nbFilledBytes;
1726
1727    if (ec_dec_bit_prob(dec, 32768))
1728    {
1729 #ifdef ENABLE_POSTFILTER
1730       int qg, octave;
1731       octave = ec_dec_uint(dec, 6);
1732       postfilter_pitch = (16<<octave)+ec_dec_bits(dec, 4+octave);
1733       qg = ec_dec_bits(dec, 2);
1734       postfilter_gain = QCONST16(.125f,15)*(qg+2);
1735 #else /* ENABLE_POSTFILTER */
1736       RESTORE_STACK;
1737       return CELT_CORRUPTED_DATA;
1738 #endif /* ENABLE_POSTFILTER */
1739
1740    } else {
1741       postfilter_gain = 0;
1742       postfilter_pitch = 0;
1743    }
1744
1745    /* Decode the global flags (first symbols in the stream) */
1746    intra_ener = ec_dec_bit_prob(dec, 8192);
1747    /* Get band energies */
1748    unquant_coarse_energy(st->mode, st->start, st->end, bandE, oldBandE,
1749          intra_ener, dec, C, LM);
1750
1751    if (LM > 0)
1752       isTransient = ec_dec_bit_prob(dec, 8192);
1753    else
1754       isTransient = 0;
1755
1756    if (isTransient)
1757       shortBlocks = M;
1758    else
1759       shortBlocks = 0;
1760
1761    ALLOC(tf_res, st->mode->nbEBands, int);
1762    tf_decode(st->start, st->end, C, isTransient, tf_res, LM, dec);
1763
1764    has_fold = ec_dec_bit_prob(dec, 8192)<<1;
1765    has_fold |= ec_dec_bit_prob(dec, (has_fold>>1) ? 32768 : 49152);
1766
1767    ALLOC(pulses, st->mode->nbEBands, int);
1768    ALLOC(offsets, st->mode->nbEBands, int);
1769    ALLOC(fine_priority, st->mode->nbEBands, int);
1770
1771    for (i=0;i<st->mode->nbEBands;i++)
1772       offsets[i] = 0;
1773    for (i=0;i<st->mode->nbEBands;i++)
1774    {
1775       if (ec_dec_bit_prob(dec, 1024))
1776       {
1777          while (ec_dec_bit_prob(dec, 32768))
1778             offsets[i]++;
1779          offsets[i]++;
1780          offsets[i] *= (6<<BITRES);
1781       }
1782    }
1783
1784    ALLOC(fine_quant, st->mode->nbEBands, int);
1785    {
1786       int fl;
1787       alloc_trim = 0;
1788       fl = ec_decode_bin(dec, 7);
1789       while (trim_cdf[alloc_trim+1] <= fl)
1790          alloc_trim++;
1791       ec_dec_update(dec, trim_cdf[alloc_trim], trim_cdf[alloc_trim+1], 128);
1792    }
1793
1794    if (C==2)
1795    {
1796       dual_stereo = ec_dec_bit_prob(dec, 32768);
1797       intensity = ec_dec_uint(dec, 1+st->end-st->start);
1798    }
1799
1800    bits = len*8 - ec_dec_tell(dec, 0) - 1;
1801    codedBands = compute_allocation(st->mode, st->start, st->end, offsets, alloc_trim, bits, pulses, fine_quant, fine_priority, C, LM);
1802    
1803    unquant_fine_energy(st->mode, st->start, st->end, bandE, oldBandE, fine_quant, dec, C);
1804
1805    /* Decode fixed codebook */
1806    quant_all_bands(0, st->mode, st->start, st->end, X, C==2 ? X+N : NULL,
1807          NULL, pulses, shortBlocks, has_fold, dual_stereo, intensity, tf_res, 1,
1808          len*8, dec, LM, codedBands);
1809
1810    unquant_energy_finalise(st->mode, st->start, st->end, bandE, oldBandE,
1811          fine_quant, fine_priority, len*8-ec_dec_tell(dec, 0), dec, C);
1812
1813    log2Amp(st->mode, st->start, st->end, bandE, oldBandE, C);
1814
1815    /* Synthesis */
1816    denormalise_bands(st->mode, X, freq, bandE, effEnd, C, M);
1817
1818    CELT_MOVE(decode_mem[0], decode_mem[0]+N, DECODE_BUFFER_SIZE-N);
1819    if (C==2)
1820       CELT_MOVE(decode_mem[1], decode_mem[1]+N, DECODE_BUFFER_SIZE-N);
1821
1822    c=0; do
1823       for (i=0;i<M*st->mode->eBands[st->start];i++)
1824          freq[c*N+i] = 0;
1825    while (++c<C);
1826    c=0; do
1827       for (i=M*st->mode->eBands[effEnd];i<N;i++)
1828          freq[c*N+i] = 0;
1829    while (++c<C);
1830
1831    out_syn[0] = out_mem[0]+MAX_PERIOD-N;
1832    if (C==2)
1833       out_syn[1] = out_mem[1]+MAX_PERIOD-N;
1834
1835    /* Compute inverse MDCTs */
1836    compute_inv_mdcts(st->mode, shortBlocks, freq, out_syn, overlap_mem, C, LM);
1837
1838 #ifdef ENABLE_POSTFILTER
1839    c=0; do {
1840       comb_filter(out_syn[c], out_syn[c], st->postfilter_period, st->postfilter_period, st->overlap, C,
1841             st->postfilter_gain, st->postfilter_gain, NULL, 0);
1842       comb_filter(out_syn[c]+st->overlap, out_syn[c]+st->overlap, st->postfilter_period, postfilter_pitch, N-st->overlap, C,
1843             st->postfilter_gain, postfilter_gain, st->mode->window, st->mode->overlap);
1844    } while (++c<C);
1845    st->postfilter_period = postfilter_pitch;
1846    st->postfilter_gain = postfilter_gain;
1847 #endif /* ENABLE_POSTFILTER */
1848
1849    deemphasis(out_syn, pcm, N, C, st->mode->preemph, st->preemph_memD);
1850    st->loss_count = 0;
1851    RESTORE_STACK;
1852    if (ec_dec_get_error(dec))
1853       return CELT_CORRUPTED_DATA;
1854    else
1855       return CELT_OK;
1856 }
1857
1858 #ifdef FIXED_POINT
1859 #ifndef DISABLE_FLOAT_API
1860 int celt_decode_with_ec_float(CELTDecoder * restrict st, const unsigned char *data, int len, float * restrict pcm, int frame_size, ec_dec *dec)
1861 {
1862    int j, ret, C, N, LM, M;
1863    VARDECL(celt_int16, out);
1864    SAVE_STACK;
1865
1866    if (pcm==NULL)
1867       return CELT_BAD_ARG;
1868
1869    for (LM=0;LM<4;LM++)
1870       if (st->mode->shortMdctSize<<LM==frame_size)
1871          break;
1872    if (LM>=MAX_CONFIG_SIZES)
1873       return CELT_BAD_ARG;
1874    M=1<<LM;
1875
1876    C = CHANNELS(st->channels);
1877    N = M*st->mode->shortMdctSize;
1878    
1879    ALLOC(out, C*N, celt_int16);
1880    ret=celt_decode_with_ec(st, data, len, out, frame_size, dec);
1881    if (ret==0)
1882       for (j=0;j<C*N;j++)
1883          pcm[j]=out[j]*(1.f/32768.f);
1884      
1885    RESTORE_STACK;
1886    return ret;
1887 }
1888 #endif /*DISABLE_FLOAT_API*/
1889 #else
1890 int celt_decode_with_ec(CELTDecoder * restrict st, const unsigned char *data, int len, celt_int16 * restrict pcm, int frame_size, ec_dec *dec)
1891 {
1892    int j, ret, C, N, LM, M;
1893    VARDECL(celt_sig, out);
1894    SAVE_STACK;
1895
1896    if (pcm==NULL)
1897       return CELT_BAD_ARG;
1898
1899    for (LM=0;LM<4;LM++)
1900       if (st->mode->shortMdctSize<<LM==frame_size)
1901          break;
1902    if (LM>=MAX_CONFIG_SIZES)
1903       return CELT_BAD_ARG;
1904    M=1<<LM;
1905
1906    C = CHANNELS(st->channels);
1907    N = M*st->mode->shortMdctSize;
1908    ALLOC(out, C*N, celt_sig);
1909
1910    ret=celt_decode_with_ec_float(st, data, len, out, frame_size, dec);
1911
1912    if (ret==0)
1913       for (j=0;j<C*N;j++)
1914          pcm[j] = FLOAT2INT16 (out[j]);
1915    
1916    RESTORE_STACK;
1917    return ret;
1918 }
1919 #endif
1920
1921 int celt_decode(CELTDecoder * restrict st, const unsigned char *data, int len, celt_int16 * restrict pcm, int frame_size)
1922 {
1923    return celt_decode_with_ec(st, data, len, pcm, frame_size, NULL);
1924 }
1925
1926 #ifndef DISABLE_FLOAT_API
1927 int celt_decode_float(CELTDecoder * restrict st, const unsigned char *data, int len, float * restrict pcm, int frame_size)
1928 {
1929    return celt_decode_with_ec_float(st, data, len, pcm, frame_size, NULL);
1930 }
1931 #endif /* DISABLE_FLOAT_API */
1932
1933 int celt_decoder_ctl(CELTDecoder * restrict st, int request, ...)
1934 {
1935    va_list ap;
1936
1937    va_start(ap, request);
1938    switch (request)
1939    {
1940       case CELT_GET_MODE_REQUEST:
1941       {
1942          const CELTMode ** value = va_arg(ap, const CELTMode**);
1943          if (value==0)
1944             goto bad_arg;
1945          *value=st->mode;
1946       }
1947       break;
1948       case CELT_SET_START_BAND_REQUEST:
1949       {
1950          celt_int32 value = va_arg(ap, celt_int32);
1951          if (value<0 || value>=st->mode->nbEBands)
1952             goto bad_arg;
1953          st->start = value;
1954       }
1955       break;
1956       case CELT_SET_END_BAND_REQUEST:
1957       {
1958          celt_int32 value = va_arg(ap, celt_int32);
1959          if (value<0 || value>=st->mode->nbEBands)
1960             goto bad_arg;
1961          st->end = value;
1962       }
1963       break;
1964       case CELT_RESET_STATE:
1965       {
1966          CELT_MEMSET((char*)&st->DECODER_RESET_START, 0,
1967                celt_decoder_get_size(st->mode, st->channels)-
1968                ((char*)&st->DECODER_RESET_START - (char*)st));
1969       }
1970       break;
1971       default:
1972          goto bad_request;
1973    }
1974    va_end(ap);
1975    return CELT_OK;
1976 bad_arg:
1977    va_end(ap);
1978    return CELT_BAD_ARG;
1979 bad_request:
1980       va_end(ap);
1981   return CELT_UNIMPLEMENTED;
1982 }
1983
1984 const char *celt_strerror(int error)
1985 {
1986    static const char *error_strings[8] = {
1987       "success",
1988       "invalid argument",
1989       "invalid mode",
1990       "internal error",
1991       "corrupted stream",
1992       "request not implemented",
1993       "invalid state",
1994       "memory allocation failed"
1995    };
1996    if (error > 0 || error < -7)
1997       return "unknown error";
1998    else 
1999       return error_strings[-error];
2000 }
2001