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