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