Fix a minor, but bitstream-affecting bug
[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 && LM>0)
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    if (LM > 0)
750       ec_enc_bit_prob(enc, shortBlocks!=0, 8192);
751
752    if (shortBlocks)
753    {
754       if (transient_shift)
755       {
756          int max_time = (N+st->mode->overlap)*(celt_int32)8000/st->mode->Fs;
757          ec_enc_uint(enc, transient_shift, 4);
758          ec_enc_uint(enc, transient_time_quant, max_time);
759       } else {
760          ec_enc_uint(enc, mdct_weight_shift, 4);
761          if (mdct_weight_shift && M!=2)
762             ec_enc_uint(enc, mdct_weight_pos, M-1);
763       }
764    }
765
766    tf_encode(st->start, st->end, isTransient, tf_res, nbAvailableBytes, LM, tf_select, enc);
767
768    if (shortBlocks || st->complexity < 3)
769    {
770       if (st->complexity == 0)
771       {
772          has_fold = 0;
773          st->fold_decision = 3;
774       } else {
775          has_fold = 1;
776          st->fold_decision = 1;
777       }
778    } else {
779       has_fold = folding_decision(st->mode, X, &st->tonal_average, &st->fold_decision, effEnd, C, M);
780    }
781    ec_enc_bit_prob(enc, has_fold>>1, 8192);
782    ec_enc_bit_prob(enc, has_fold&1, (has_fold>>1) ? 32768 : 49152);
783
784    /* Variable bitrate */
785    if (st->vbr_rate_norm>0)
786    {
787      celt_word16 alpha;
788      celt_int32 delta;
789      /* The target rate in 16th bits per frame */
790      celt_int32 vbr_rate;
791      celt_int32 target;
792      celt_int32 vbr_bound, max_allowed;
793
794      vbr_rate = M*st->vbr_rate_norm;
795
796      /* Computes the max bit-rate allowed in VBR more to avoid busting the budget */
797      vbr_bound = vbr_rate;
798      max_allowed = (vbr_rate + vbr_bound - st->vbr_reservoir)>>(BITRES+3);
799      if (max_allowed < 4)
800         max_allowed = 4;
801      if (max_allowed < nbAvailableBytes)
802         nbAvailableBytes = max_allowed;
803      target=vbr_rate;
804
805      /* Shortblocks get a large boost in bitrate, but since they 
806         are uncommon long blocks are not greatly effected */
807      if (shortBlocks)
808        target*=2;
809      else if (M > 1)
810        target-=(target+14)/28;
811
812      /* The average energy is removed from the target and the actual 
813         energy added*/
814      target=target+st->vbr_offset-588+ec_enc_tell(enc, BITRES);
815
816      /* In VBR mode the frame size must not be reduced so much that it would result in the coarse energy busting its budget */
817      target=IMIN(nbAvailableBytes,target);
818      /* Make the adaptation coef (alpha) higher at the beginning */
819      if (st->vbr_count < 990)
820      {
821         st->vbr_count++;
822         alpha = celt_rcp(SHL32(EXTEND32(st->vbr_count+10),16));
823         /*printf ("%d %d\n", st->vbr_count+10, alpha);*/
824      } else
825         alpha = QCONST16(.001f,15);
826
827      /* By how much did we "miss" the target on that frame */
828      delta = (8<<BITRES)*(celt_int32)target - vbr_rate;
829      /* How many bits have we used in excess of what we're allowed */
830      st->vbr_reservoir += delta;
831      /*printf ("%d\n", st->vbr_reservoir);*/
832
833      /* Compute the offset we need to apply in order to reach the target */
834      st->vbr_drift += (celt_int32)MULT16_32_Q15(alpha,delta-st->vbr_offset-st->vbr_drift);
835      st->vbr_offset = -st->vbr_drift;
836      /*printf ("%d\n", st->vbr_drift);*/
837
838      /* We could use any multiple of vbr_rate as bound (depending on the delay) */
839      if (st->vbr_reservoir < 0)
840      {
841         /* We're under the min value -- increase rate */
842         int adjust = 1-(st->vbr_reservoir-1)/(8<<BITRES);
843         st->vbr_reservoir += adjust*(8<<BITRES);
844         target += adjust;
845         /*printf ("+%d\n", adjust);*/
846      }
847      if (target < nbAvailableBytes)
848         nbAvailableBytes = target;
849      nbCompressedBytes = nbAvailableBytes + nbFilledBytes;
850
851      /* This moves the raw bits to take into account the new compressed size */
852      ec_byte_shrink(&buf, nbCompressedBytes);
853    }
854
855    /* Bit allocation */
856    ALLOC(fine_quant, st->mode->nbEBands, int);
857    ALLOC(pulses, st->mode->nbEBands, int);
858    ALLOC(offsets, st->mode->nbEBands, int);
859    ALLOC(fine_priority, st->mode->nbEBands, int);
860
861    for (i=0;i<st->mode->nbEBands;i++)
862       offsets[i] = 0;
863    bits = nbCompressedBytes*8 - ec_enc_tell(enc, 0) - 1;
864    compute_allocation(st->mode, st->start, st->end, offsets, bits, pulses, fine_quant, fine_priority, C, LM);
865
866    quant_fine_energy(st->mode, st->start, st->end, bandE, oldBandE, error, fine_quant, enc, C);
867
868 #ifdef MEASURE_NORM_MSE
869    float X0[3000];
870    float bandE0[60];
871    for (c=0;c<C;c++)
872       for (i=0;i<N;i++)
873          X0[i+c*N] = X[i+c*N];
874    for (i=0;i<C*st->mode->nbEBands;i++)
875       bandE0[i] = bandE[i];
876 #endif
877
878    /* Residual quantisation */
879    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);
880
881    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);
882
883    /* Re-synthesis of the coded audio if required */
884    if (resynth)
885    {
886       VARDECL(celt_sig, _out_mem);
887       celt_sig *out_mem[2];
888       celt_sig *overlap_mem[2];
889
890       log2Amp(st->mode, st->start, st->end, bandE, oldBandE, C);
891
892 #ifdef MEASURE_NORM_MSE
893       measure_norm_mse(st->mode, X, X0, bandE, bandE0, M, N, C);
894 #endif
895
896       if (mdct_weight_shift)
897       {
898          mdct_shape(st->mode, X, 0, mdct_weight_pos+1, N, mdct_weight_shift, effEnd, C, 1, M);
899       }
900
901       /* Synthesis */
902       denormalise_bands(st->mode, X, freq, bandE, effEnd, C, M);
903
904       for (c=0;c<C;c++)
905          for (i=0;i<M*st->mode->eBands[st->start];i++)
906             freq[c*N+i] = 0;
907       for (c=0;c<C;c++)
908          for (i=M*st->mode->eBands[st->end];i<N;i++)
909             freq[c*N+i] = 0;
910
911       ALLOC(_out_mem, C*N, celt_sig);
912
913       for (c=0;c<C;c++)
914       {
915          overlap_mem[c] = _overlap_mem + c*st->overlap;
916          out_mem[c] = _out_mem+c*N;
917       }
918
919       compute_inv_mdcts(st->mode, shortBlocks, freq, transient_time,
920             transient_shift, out_mem, overlap_mem, C, LM);
921
922       /* De-emphasis and put everything back at the right place 
923          in the synthesis history */
924       if (optional_resynthesis != NULL) {
925          deemphasis(out_mem, optional_resynthesis, N, C, st->mode->preemph, st->preemph_memD);
926
927       }
928    }
929
930    /* If there's any room left (can only happen for very high rates),
931       fill it with zeros */
932    while (ec_enc_tell(enc,0) + 8 <= nbCompressedBytes*8)
933       ec_enc_bits(enc, 0, 8);
934    ec_enc_done(enc);
935    
936    RESTORE_STACK;
937    if (ec_enc_get_error(enc))
938       return CELT_CORRUPTED_DATA;
939    else
940       return nbCompressedBytes;
941 }
942
943 #ifdef FIXED_POINT
944 #ifndef DISABLE_FLOAT_API
945 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)
946 {
947    int j, ret, C, N, LM, M;
948    VARDECL(celt_int16, in);
949    SAVE_STACK;
950
951    if (pcm==NULL)
952       return CELT_BAD_ARG;
953
954    for (LM=0;LM<4;LM++)
955       if (st->mode->shortMdctSize<<LM==frame_size)
956          break;
957    if (LM>=MAX_CONFIG_SIZES)
958       return CELT_BAD_ARG;
959    M=1<<LM;
960
961    C = CHANNELS(st->channels);
962    N = M*st->mode->shortMdctSize;
963    ALLOC(in, C*N, celt_int16);
964
965    for (j=0;j<C*N;j++)
966      in[j] = FLOAT2INT16(pcm[j]);
967
968    if (optional_resynthesis != NULL) {
969      ret=celt_encode_with_ec(st,in,in,frame_size,compressed,nbCompressedBytes, enc);
970       for (j=0;j<C*N;j++)
971          optional_resynthesis[j]=in[j]*(1.f/32768.f);
972    } else {
973      ret=celt_encode_with_ec(st,in,NULL,frame_size,compressed,nbCompressedBytes, enc);
974    }
975    RESTORE_STACK;
976    return ret;
977
978 }
979 #endif /*DISABLE_FLOAT_API*/
980 #else
981 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)
982 {
983    int j, ret, C, N, LM, M;
984    VARDECL(celt_sig, in);
985    SAVE_STACK;
986
987    if (pcm==NULL)
988       return CELT_BAD_ARG;
989
990    for (LM=0;LM<4;LM++)
991       if (st->mode->shortMdctSize<<LM==frame_size)
992          break;
993    if (LM>=MAX_CONFIG_SIZES)
994       return CELT_BAD_ARG;
995    M=1<<LM;
996
997    C=CHANNELS(st->channels);
998    N=M*st->mode->shortMdctSize;
999    ALLOC(in, C*N, celt_sig);
1000    for (j=0;j<C*N;j++) {
1001      in[j] = SCALEOUT(pcm[j]);
1002    }
1003
1004    if (optional_resynthesis != NULL) {
1005       ret = celt_encode_with_ec_float(st,in,in,frame_size,compressed,nbCompressedBytes, enc);
1006       for (j=0;j<C*N;j++)
1007          optional_resynthesis[j] = FLOAT2INT16(in[j]);
1008    } else {
1009       ret = celt_encode_with_ec_float(st,in,NULL,frame_size,compressed,nbCompressedBytes, enc);
1010    }
1011    RESTORE_STACK;
1012    return ret;
1013 }
1014 #endif
1015
1016 int celt_encode(CELTEncoder * restrict st, const celt_int16 * pcm, int frame_size, unsigned char *compressed, int nbCompressedBytes)
1017 {
1018    return celt_encode_with_ec(st, pcm, NULL, frame_size, compressed, nbCompressedBytes, NULL);
1019 }
1020
1021 #ifndef DISABLE_FLOAT_API
1022 int celt_encode_float(CELTEncoder * restrict st, const float * pcm, int frame_size, unsigned char *compressed, int nbCompressedBytes)
1023 {
1024    return celt_encode_with_ec_float(st, pcm, NULL, frame_size, compressed, nbCompressedBytes, NULL);
1025 }
1026 #endif /* DISABLE_FLOAT_API */
1027
1028 int celt_encode_resynthesis(CELTEncoder * restrict st, const celt_int16 * pcm, celt_int16 * optional_resynthesis, int frame_size, unsigned char *compressed, int nbCompressedBytes)
1029 {
1030    return celt_encode_with_ec(st, pcm, optional_resynthesis, frame_size, compressed, nbCompressedBytes, NULL);
1031 }
1032
1033 #ifndef DISABLE_FLOAT_API
1034 int celt_encode_resynthesis_float(CELTEncoder * restrict st, const float * pcm, float * optional_resynthesis, int frame_size, unsigned char *compressed, int nbCompressedBytes)
1035 {
1036    return celt_encode_with_ec_float(st, pcm, optional_resynthesis, frame_size, compressed, nbCompressedBytes, NULL);
1037 }
1038 #endif /* DISABLE_FLOAT_API */
1039
1040
1041 int celt_encoder_ctl(CELTEncoder * restrict st, int request, ...)
1042 {
1043    va_list ap;
1044    
1045    va_start(ap, request);
1046    switch (request)
1047    {
1048       case CELT_GET_MODE_REQUEST:
1049       {
1050          const CELTMode ** value = va_arg(ap, const CELTMode**);
1051          if (value==0)
1052             goto bad_arg;
1053          *value=st->mode;
1054       }
1055       break;
1056       case CELT_SET_COMPLEXITY_REQUEST:
1057       {
1058          int value = va_arg(ap, celt_int32);
1059          if (value<0 || value>10)
1060             goto bad_arg;
1061          st->complexity = value;
1062       }
1063       break;
1064       case CELT_SET_START_BAND_REQUEST:
1065       {
1066          celt_int32 value = va_arg(ap, celt_int32);
1067          if (value<0 || value>=st->mode->nbEBands)
1068             goto bad_arg;
1069          st->start = value;
1070       }
1071       break;
1072       case CELT_SET_END_BAND_REQUEST:
1073       {
1074          celt_int32 value = va_arg(ap, celt_int32);
1075          if (value<0 || value>=st->mode->nbEBands)
1076             goto bad_arg;
1077          st->end = value;
1078       }
1079       break;
1080       case CELT_SET_PREDICTION_REQUEST:
1081       {
1082          int value = va_arg(ap, celt_int32);
1083          if (value<0 || value>2)
1084             goto bad_arg;
1085          if (value==0)
1086          {
1087             st->force_intra   = 1;
1088          } else if (value==1) {
1089             st->force_intra   = 0;
1090          } else {
1091             st->force_intra   = 0;
1092          }   
1093       }
1094       break;
1095       case CELT_SET_VBR_RATE_REQUEST:
1096       {
1097          celt_int32 value = va_arg(ap, celt_int32);
1098          int frame_rate;
1099          int N = st->mode->shortMdctSize;
1100          if (value<0)
1101             goto bad_arg;
1102          if (value>3072000)
1103             value = 3072000;
1104          frame_rate = ((st->mode->Fs<<3)+(N>>1))/N;
1105          st->vbr_rate_norm = ((value<<(BITRES+3))+(frame_rate>>1))/frame_rate;
1106       }
1107       break;
1108       case CELT_RESET_STATE:
1109       {
1110          CELT_MEMSET((char*)&st->ENCODER_RESET_START, 0,
1111                celt_encoder_get_size(st->mode, st->channels)-
1112                ((char*)&st->ENCODER_RESET_START - (char*)st));
1113          st->delayedIntra = 1;
1114          st->fold_decision = 1;
1115          st->tonal_average = QCONST16(1.f,8);
1116       }
1117       break;
1118       default:
1119          goto bad_request;
1120    }
1121    va_end(ap);
1122    return CELT_OK;
1123 bad_arg:
1124    va_end(ap);
1125    return CELT_BAD_ARG;
1126 bad_request:
1127    va_end(ap);
1128    return CELT_UNIMPLEMENTED;
1129 }
1130
1131 /**********************************************************************/
1132 /*                                                                    */
1133 /*                             DECODER                                */
1134 /*                                                                    */
1135 /**********************************************************************/
1136 #define DECODE_BUFFER_SIZE 2048
1137
1138 /** Decoder state 
1139  @brief Decoder state
1140  */
1141 struct CELTDecoder {
1142    const CELTMode *mode;
1143    int overlap;
1144    int channels;
1145
1146    int start, end;
1147
1148    /* Everything beyond this point gets cleared on a reset */
1149 #define DECODER_RESET_START last_pitch_index
1150
1151    int last_pitch_index;
1152    int loss_count;
1153
1154    celt_sig preemph_memD[2];
1155    
1156    celt_sig _decode_mem[1]; /* Size = channels*(DECODE_BUFFER_SIZE+mode->overlap) */
1157    /* celt_word16 lpc[],  Size = channels*LPC_ORDER */
1158    /* celt_word16 oldEBands[], Size = channels*mode->nbEBands */
1159 };
1160
1161 int celt_decoder_get_size(const CELTMode *mode, int channels)
1162 {
1163    int size = sizeof(struct CELTDecoder)
1164             + (channels*(DECODE_BUFFER_SIZE+mode->overlap)-1)*sizeof(celt_sig)
1165             + channels*LPC_ORDER*sizeof(celt_word16)
1166             + channels*mode->nbEBands*sizeof(celt_word16);
1167    return size;
1168 }
1169
1170 CELTDecoder *celt_decoder_create(const CELTMode *mode, int channels, int *error)
1171 {
1172    return celt_decoder_init(
1173          (CELTDecoder *)celt_alloc(celt_decoder_get_size(mode, channels)),
1174          mode, channels, error);
1175 }
1176
1177 CELTDecoder *celt_decoder_init(CELTDecoder *st, const CELTMode *mode, int channels, int *error)
1178 {
1179    if (channels < 0 || channels > 2)
1180    {
1181       if (error)
1182          *error = CELT_BAD_ARG;
1183       return NULL;
1184    }
1185
1186    if (st==NULL)
1187    {
1188       if (error)
1189          *error = CELT_ALLOC_FAIL;
1190       return NULL;
1191    }
1192
1193    CELT_MEMSET((char*)st, 0, celt_decoder_get_size(mode, channels));
1194
1195    st->mode = mode;
1196    st->overlap = mode->overlap;
1197    st->channels = channels;
1198
1199    st->start = 0;
1200    st->end = st->mode->effEBands;
1201
1202    st->loss_count = 0;
1203
1204    if (error)
1205       *error = CELT_OK;
1206    return st;
1207 }
1208
1209 void celt_decoder_destroy(CELTDecoder *st)
1210 {
1211    celt_free(st);
1212 }
1213
1214 static void celt_decode_lost(CELTDecoder * restrict st, celt_word16 * restrict pcm, int N, int LM)
1215 {
1216    int c;
1217    int pitch_index;
1218    int overlap = st->mode->overlap;
1219    celt_word16 fade = Q15ONE;
1220    int i, len;
1221    const int C = CHANNELS(st->channels);
1222    int offset;
1223    celt_sig *out_mem[2];
1224    celt_sig *decode_mem[2];
1225    celt_sig *overlap_mem[2];
1226    celt_word16 *lpc;
1227    celt_word16 *oldBandE;
1228    SAVE_STACK;
1229    
1230    for (c=0;c<C;c++)
1231    {
1232       decode_mem[c] = st->_decode_mem + c*(DECODE_BUFFER_SIZE+st->overlap);
1233       out_mem[c] = decode_mem[c]+DECODE_BUFFER_SIZE-MAX_PERIOD;
1234       overlap_mem[c] = decode_mem[c]+DECODE_BUFFER_SIZE;
1235    }
1236    lpc = (celt_word16*)(st->_decode_mem+(DECODE_BUFFER_SIZE+st->overlap)*C);
1237    oldBandE = lpc+C*LPC_ORDER;
1238
1239    len = N+st->mode->overlap;
1240    
1241    if (st->loss_count == 0)
1242    {
1243       celt_word16 pitch_buf[MAX_PERIOD>>1];
1244       celt_word32 tmp=0;
1245       celt_word32 mem0[2]={0,0};
1246       celt_word16 mem1[2]={0,0};
1247       int len2 = len;
1248       /* FIXME: This is a kludge */
1249       if (len2>MAX_PERIOD>>1)
1250          len2 = MAX_PERIOD>>1;
1251       pitch_downsample(out_mem, pitch_buf, MAX_PERIOD, MAX_PERIOD,
1252                        C, mem0, mem1);
1253       pitch_search(st->mode, pitch_buf+((MAX_PERIOD-len2)>>1), pitch_buf, len2,
1254                    MAX_PERIOD-len2-100, &pitch_index, &tmp, 1<<LM);
1255       pitch_index = MAX_PERIOD-len2-pitch_index;
1256       st->last_pitch_index = pitch_index;
1257    } else {
1258       pitch_index = st->last_pitch_index;
1259       if (st->loss_count < 5)
1260          fade = QCONST16(.8f,15);
1261       else
1262          fade = 0;
1263    }
1264
1265    for (c=0;c<C;c++)
1266    {
1267       /* FIXME: This is more memory than necessary */
1268       celt_word32 e[2*MAX_PERIOD];
1269       celt_word16 exc[2*MAX_PERIOD];
1270       celt_word32 ac[LPC_ORDER+1];
1271       celt_word16 decay = 1;
1272       celt_word32 S1=0;
1273       celt_word16 mem[LPC_ORDER]={0};
1274
1275       offset = MAX_PERIOD-pitch_index;
1276       for (i=0;i<MAX_PERIOD;i++)
1277          exc[i] = ROUND16(out_mem[c][i], SIG_SHIFT);
1278
1279       if (st->loss_count == 0)
1280       {
1281          _celt_autocorr(exc, ac, st->mode->window, st->mode->overlap,
1282                         LPC_ORDER, MAX_PERIOD);
1283
1284          /* Noise floor -40 dB */
1285 #ifdef FIXED_POINT
1286          ac[0] += SHR32(ac[0],13);
1287 #else
1288          ac[0] *= 1.0001f;
1289 #endif
1290          /* Lag windowing */
1291          for (i=1;i<=LPC_ORDER;i++)
1292          {
1293             /*ac[i] *= exp(-.5*(2*M_PI*.002*i)*(2*M_PI*.002*i));*/
1294 #ifdef FIXED_POINT
1295             ac[i] -= MULT16_32_Q15(2*i*i, ac[i]);
1296 #else
1297             ac[i] -= ac[i]*(.008f*i)*(.008f*i);
1298 #endif
1299          }
1300
1301          _celt_lpc(lpc+c*LPC_ORDER, ac, LPC_ORDER);
1302       }
1303       fir(exc, lpc+c*LPC_ORDER, exc, MAX_PERIOD, LPC_ORDER, mem);
1304       /*for (i=0;i<MAX_PERIOD;i++)printf("%d ", exc[i]); printf("\n");*/
1305       /* Check if the waveform is decaying (and if so how fast) */
1306       {
1307          celt_word32 E1=1, E2=1;
1308          int period;
1309          if (pitch_index <= MAX_PERIOD/2)
1310             period = pitch_index;
1311          else
1312             period = MAX_PERIOD/2;
1313          for (i=0;i<period;i++)
1314          {
1315             E1 += SHR32(MULT16_16(exc[MAX_PERIOD-period+i],exc[MAX_PERIOD-period+i]),8);
1316             E2 += SHR32(MULT16_16(exc[MAX_PERIOD-2*period+i],exc[MAX_PERIOD-2*period+i]),8);
1317          }
1318          if (E1 > E2)
1319             E1 = E2;
1320          decay = celt_sqrt(frac_div32(SHR(E1,1),E2));
1321       }
1322
1323       /* Copy excitation, taking decay into account */
1324       for (i=0;i<len+st->mode->overlap;i++)
1325       {
1326          if (offset+i >= MAX_PERIOD)
1327          {
1328             offset -= pitch_index;
1329             decay = MULT16_16_Q15(decay, decay);
1330          }
1331          e[i] = SHL32(EXTEND32(MULT16_16_Q15(decay, exc[offset+i])), SIG_SHIFT);
1332          S1 += SHR32(MULT16_16(out_mem[c][offset+i],out_mem[c][offset+i]),8);
1333       }
1334
1335       iir(e, lpc+c*LPC_ORDER, e, len+st->mode->overlap, LPC_ORDER, mem);
1336
1337       {
1338          celt_word32 S2=0;
1339          for (i=0;i<len+overlap;i++)
1340             S2 += SHR32(MULT16_16(e[i],e[i]),8);
1341          /* This checks for an "explosion" in the synthesis */
1342 #ifdef FIXED_POINT
1343          if (!(S1 > SHR32(S2,2)))
1344 #else
1345          /* Float test is written this way to catch NaNs at the same time */
1346          if (!(S1 > 0.2f*S2))
1347 #endif
1348          {
1349             for (i=0;i<len+overlap;i++)
1350                e[i] = 0;
1351          } else if (S1 < S2)
1352          {
1353             celt_word16 ratio = celt_sqrt(frac_div32(SHR32(S1,1)+1,S2+1));
1354             for (i=0;i<len+overlap;i++)
1355                e[i] = MULT16_16_Q15(ratio, e[i]);
1356          }
1357       }
1358
1359       for (i=0;i<MAX_PERIOD+st->mode->overlap-N;i++)
1360          out_mem[c][i] = out_mem[c][N+i];
1361
1362       /* Apply TDAC to the concealed audio so that it blends with the
1363          previous and next frames */
1364       for (i=0;i<overlap/2;i++)
1365       {
1366          celt_word32 tmp1, tmp2;
1367          tmp1 = MULT16_32_Q15(st->mode->window[i          ], e[i          ]) -
1368                 MULT16_32_Q15(st->mode->window[overlap-i-1], e[overlap-i-1]);
1369          tmp2 = MULT16_32_Q15(st->mode->window[i],           e[N+overlap-1-i]) +
1370                 MULT16_32_Q15(st->mode->window[overlap-i-1], e[N+i          ]);
1371          tmp1 = MULT16_32_Q15(fade, tmp1);
1372          tmp2 = MULT16_32_Q15(fade, tmp2);
1373          out_mem[c][MAX_PERIOD+i] = MULT16_32_Q15(st->mode->window[overlap-i-1], tmp2);
1374          out_mem[c][MAX_PERIOD+overlap-i-1] = MULT16_32_Q15(st->mode->window[i], tmp2);
1375          out_mem[c][MAX_PERIOD-N+i] += MULT16_32_Q15(st->mode->window[i], tmp1);
1376          out_mem[c][MAX_PERIOD-N+overlap-i-1] -= MULT16_32_Q15(st->mode->window[overlap-i-1], tmp1);
1377       }
1378       for (i=0;i<N-overlap;i++)
1379          out_mem[c][MAX_PERIOD-N+overlap+i] = MULT16_32_Q15(fade, e[overlap+i]);
1380    }
1381
1382    deemphasis(out_mem, pcm, N, C, st->mode->preemph, st->preemph_memD);
1383    
1384    st->loss_count++;
1385
1386    RESTORE_STACK;
1387 }
1388
1389 #ifdef FIXED_POINT
1390 int celt_decode_with_ec(CELTDecoder * restrict st, const unsigned char *data, int len, celt_int16 * restrict pcm, int frame_size, ec_dec *dec)
1391 {
1392 #else
1393 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)
1394 {
1395 #endif
1396    int c, i, N;
1397    int has_fold;
1398    int bits;
1399    ec_dec _dec;
1400    ec_byte_buffer buf;
1401    VARDECL(celt_sig, freq);
1402    VARDECL(celt_norm, X);
1403    VARDECL(celt_ener, bandE);
1404    VARDECL(int, fine_quant);
1405    VARDECL(int, pulses);
1406    VARDECL(int, offsets);
1407    VARDECL(int, fine_priority);
1408    VARDECL(int, tf_res);
1409    celt_sig *out_mem[2];
1410    celt_sig *decode_mem[2];
1411    celt_sig *overlap_mem[2];
1412    celt_sig *out_syn[2];
1413    celt_word16 *lpc;
1414    celt_word16 *oldBandE;
1415
1416    int shortBlocks;
1417    int isTransient;
1418    int intra_ener;
1419    int transient_time;
1420    int transient_shift;
1421    int mdct_weight_shift=0;
1422    const int C = CHANNELS(st->channels);
1423    int mdct_weight_pos=0;
1424    int LM, M;
1425    int nbFilledBytes, nbAvailableBytes;
1426    int effEnd;
1427    SAVE_STACK;
1428
1429    if (pcm==NULL)
1430       return CELT_BAD_ARG;
1431
1432    for (LM=0;LM<4;LM++)
1433       if (st->mode->shortMdctSize<<LM==frame_size)
1434          break;
1435    if (LM>=MAX_CONFIG_SIZES)
1436       return CELT_BAD_ARG;
1437    M=1<<LM;
1438
1439    for (c=0;c<C;c++)
1440    {
1441       decode_mem[c] = st->_decode_mem + c*(DECODE_BUFFER_SIZE+st->overlap);
1442       out_mem[c] = decode_mem[c]+DECODE_BUFFER_SIZE-MAX_PERIOD;
1443       overlap_mem[c] = decode_mem[c]+DECODE_BUFFER_SIZE;
1444    }
1445    lpc = (celt_word16*)(st->_decode_mem+(DECODE_BUFFER_SIZE+st->overlap)*C);
1446    oldBandE = lpc+C*LPC_ORDER;
1447
1448    N = M*st->mode->shortMdctSize;
1449
1450    effEnd = st->end;
1451    if (effEnd > st->mode->effEBands)
1452       effEnd = st->mode->effEBands;
1453
1454    ALLOC(freq, C*N, celt_sig); /**< Interleaved signal MDCTs */
1455    ALLOC(X, C*N, celt_norm);   /**< Interleaved normalised MDCTs */
1456    ALLOC(bandE, st->mode->nbEBands*C, celt_ener);
1457    for (c=0;c<C;c++)
1458       for (i=0;i<M*st->mode->eBands[st->start];i++)
1459          X[c*N+i] = 0;
1460    for (c=0;c<C;c++)
1461       for (i=M*st->mode->eBands[effEnd];i<N;i++)
1462          X[c*N+i] = 0;
1463
1464    if (data == NULL)
1465    {
1466       celt_decode_lost(st, pcm, N, LM);
1467       RESTORE_STACK;
1468       return CELT_OK;
1469    }
1470    if (len<0) {
1471      RESTORE_STACK;
1472      return CELT_BAD_ARG;
1473    }
1474    
1475    if (dec == NULL)
1476    {
1477       ec_byte_readinit(&buf,(unsigned char*)data,len);
1478       ec_dec_init(&_dec,&buf);
1479       dec = &_dec;
1480       nbFilledBytes = 0;
1481    } else {
1482       nbFilledBytes = (ec_dec_tell(dec, 0)+4)>>3;
1483    }
1484    nbAvailableBytes = len-nbFilledBytes;
1485
1486    /* Decode the global flags (first symbols in the stream) */
1487    intra_ener = ec_dec_bit_prob(dec, 8192);
1488    /* Get band energies */
1489    unquant_coarse_energy(st->mode, st->start, st->end, bandE, oldBandE,
1490          intra_ener, st->mode->prob, dec, C, LM);
1491
1492    if (LM > 0)
1493       isTransient = ec_dec_bit_prob(dec, 8192);
1494    else
1495       isTransient = 0;
1496
1497    if (isTransient)
1498       shortBlocks = M;
1499    else
1500       shortBlocks = 0;
1501
1502    if (isTransient)
1503    {
1504       transient_shift = ec_dec_uint(dec, 4);
1505       if (transient_shift == 3)
1506       {
1507          int transient_time_quant;
1508          int max_time = (N+st->mode->overlap)*(celt_int32)8000/st->mode->Fs;
1509          transient_time_quant = ec_dec_uint(dec, max_time);
1510          transient_time = transient_time_quant*(celt_int32)st->mode->Fs/8000;
1511       } else {
1512          mdct_weight_shift = transient_shift;
1513          if (mdct_weight_shift && M>2)
1514             mdct_weight_pos = ec_dec_uint(dec, M-1);
1515          transient_shift = 0;
1516          transient_time = 0;
1517       }
1518    } else {
1519       transient_time = -1;
1520       transient_shift = 0;
1521    }
1522
1523    ALLOC(tf_res, st->mode->nbEBands, int);
1524    tf_decode(st->start, st->end, C, isTransient, tf_res, nbAvailableBytes, LM, dec);
1525
1526    has_fold = ec_dec_bit_prob(dec, 8192)<<1;
1527    has_fold |= ec_dec_bit_prob(dec, (has_fold>>1) ? 32768 : 49152);
1528
1529    ALLOC(pulses, st->mode->nbEBands, int);
1530    ALLOC(offsets, st->mode->nbEBands, int);
1531    ALLOC(fine_priority, st->mode->nbEBands, int);
1532
1533    for (i=0;i<st->mode->nbEBands;i++)
1534       offsets[i] = 0;
1535
1536    bits = len*8 - ec_dec_tell(dec, 0) - 1;
1537    ALLOC(fine_quant, st->mode->nbEBands, int);
1538    compute_allocation(st->mode, st->start, st->end, offsets, bits, pulses, fine_quant, fine_priority, C, LM);
1539    /*bits = ec_dec_tell(dec, 0);
1540    compute_fine_allocation(st->mode, fine_quant, (20*C+len*8/5-(ec_dec_tell(dec, 0)-bits))/C);*/
1541    
1542    unquant_fine_energy(st->mode, st->start, st->end, bandE, oldBandE, fine_quant, dec, C);
1543
1544    /* Decode fixed codebook */
1545    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);
1546
1547    unquant_energy_finalise(st->mode, st->start, st->end, bandE, oldBandE,
1548          fine_quant, fine_priority, len*8-ec_dec_tell(dec, 0), dec, C);
1549
1550    log2Amp(st->mode, st->start, st->end, bandE, oldBandE, C);
1551
1552    if (mdct_weight_shift)
1553    {
1554       mdct_shape(st->mode, X, 0, mdct_weight_pos+1, N, mdct_weight_shift, effEnd, C, 1, M);
1555    }
1556
1557    /* Synthesis */
1558    denormalise_bands(st->mode, X, freq, bandE, effEnd, C, M);
1559
1560    CELT_MOVE(decode_mem[0], decode_mem[0]+N, DECODE_BUFFER_SIZE-N);
1561    if (C==2)
1562       CELT_MOVE(decode_mem[1], decode_mem[1]+N, DECODE_BUFFER_SIZE-N);
1563
1564    for (c=0;c<C;c++)
1565       for (i=0;i<M*st->mode->eBands[st->start];i++)
1566          freq[c*N+i] = 0;
1567    for (c=0;c<C;c++)
1568       for (i=M*st->mode->eBands[effEnd];i<N;i++)
1569          freq[c*N+i] = 0;
1570
1571    out_syn[0] = out_mem[0]+MAX_PERIOD-N;
1572    if (C==2)
1573       out_syn[1] = out_mem[1]+MAX_PERIOD-N;
1574
1575    /* Compute inverse MDCTs */
1576    compute_inv_mdcts(st->mode, shortBlocks, freq, transient_time,
1577          transient_shift, out_syn, overlap_mem, C, LM);
1578
1579    deemphasis(out_syn, pcm, N, C, st->mode->preemph, st->preemph_memD);
1580    st->loss_count = 0;
1581    RESTORE_STACK;
1582    if (ec_dec_get_error(dec))
1583       return CELT_CORRUPTED_DATA;
1584    else
1585       return CELT_OK;
1586 }
1587
1588 #ifdef FIXED_POINT
1589 #ifndef DISABLE_FLOAT_API
1590 int celt_decode_with_ec_float(CELTDecoder * restrict st, const unsigned char *data, int len, float * restrict pcm, int frame_size, ec_dec *dec)
1591 {
1592    int j, ret, C, N, LM, M;
1593    VARDECL(celt_int16, out);
1594    SAVE_STACK;
1595
1596    if (pcm==NULL)
1597       return CELT_BAD_ARG;
1598
1599    for (LM=0;LM<4;LM++)
1600       if (st->mode->shortMdctSize<<LM==frame_size)
1601          break;
1602    if (LM>=MAX_CONFIG_SIZES)
1603       return CELT_BAD_ARG;
1604    M=1<<LM;
1605
1606    C = CHANNELS(st->channels);
1607    N = M*st->mode->shortMdctSize;
1608    
1609    ALLOC(out, C*N, celt_int16);
1610    ret=celt_decode_with_ec(st, data, len, out, frame_size, dec);
1611    if (ret==0)
1612       for (j=0;j<C*N;j++)
1613          pcm[j]=out[j]*(1.f/32768.f);
1614      
1615    RESTORE_STACK;
1616    return ret;
1617 }
1618 #endif /*DISABLE_FLOAT_API*/
1619 #else
1620 int celt_decode_with_ec(CELTDecoder * restrict st, const unsigned char *data, int len, celt_int16 * restrict pcm, int frame_size, ec_dec *dec)
1621 {
1622    int j, ret, C, N, LM, M;
1623    VARDECL(celt_sig, out);
1624    SAVE_STACK;
1625
1626    if (pcm==NULL)
1627       return CELT_BAD_ARG;
1628
1629    for (LM=0;LM<4;LM++)
1630       if (st->mode->shortMdctSize<<LM==frame_size)
1631          break;
1632    if (LM>=MAX_CONFIG_SIZES)
1633       return CELT_BAD_ARG;
1634    M=1<<LM;
1635
1636    C = CHANNELS(st->channels);
1637    N = M*st->mode->shortMdctSize;
1638    ALLOC(out, C*N, celt_sig);
1639
1640    ret=celt_decode_with_ec_float(st, data, len, out, frame_size, dec);
1641
1642    if (ret==0)
1643       for (j=0;j<C*N;j++)
1644          pcm[j] = FLOAT2INT16 (out[j]);
1645    
1646    RESTORE_STACK;
1647    return ret;
1648 }
1649 #endif
1650
1651 int celt_decode(CELTDecoder * restrict st, const unsigned char *data, int len, celt_int16 * restrict pcm, int frame_size)
1652 {
1653    return celt_decode_with_ec(st, data, len, pcm, frame_size, NULL);
1654 }
1655
1656 #ifndef DISABLE_FLOAT_API
1657 int celt_decode_float(CELTDecoder * restrict st, const unsigned char *data, int len, float * restrict pcm, int frame_size)
1658 {
1659    return celt_decode_with_ec_float(st, data, len, pcm, frame_size, NULL);
1660 }
1661 #endif /* DISABLE_FLOAT_API */
1662
1663 int celt_decoder_ctl(CELTDecoder * restrict st, int request, ...)
1664 {
1665    va_list ap;
1666
1667    va_start(ap, request);
1668    switch (request)
1669    {
1670       case CELT_GET_MODE_REQUEST:
1671       {
1672          const CELTMode ** value = va_arg(ap, const CELTMode**);
1673          if (value==0)
1674             goto bad_arg;
1675          *value=st->mode;
1676       }
1677       break;
1678       case CELT_SET_START_BAND_REQUEST:
1679       {
1680          celt_int32 value = va_arg(ap, celt_int32);
1681          if (value<0 || value>=st->mode->nbEBands)
1682             goto bad_arg;
1683          st->start = value;
1684       }
1685       break;
1686       case CELT_SET_END_BAND_REQUEST:
1687       {
1688          celt_int32 value = va_arg(ap, celt_int32);
1689          if (value<0 || value>=st->mode->nbEBands)
1690             goto bad_arg;
1691          st->end = value;
1692       }
1693       break;
1694       case CELT_RESET_STATE:
1695       {
1696          CELT_MEMSET((char*)&st->DECODER_RESET_START, 0,
1697                celt_decoder_get_size(st->mode, st->channels)-
1698                ((char*)&st->DECODER_RESET_START - (char*)st));
1699       }
1700       break;
1701       default:
1702          goto bad_request;
1703    }
1704    va_end(ap);
1705    return CELT_OK;
1706 bad_arg:
1707    va_end(ap);
1708    return CELT_BAD_ARG;
1709 bad_request:
1710       va_end(ap);
1711   return CELT_UNIMPLEMENTED;
1712 }
1713
1714 const char *celt_strerror(int error)
1715 {
1716    static const char *error_strings[8] = {
1717       "success",
1718       "invalid argument",
1719       "invalid mode",
1720       "internal error",
1721       "corrupted stream",
1722       "request not implemented",
1723       "invalid state",
1724       "memory allocation failed"
1725    };
1726    if (error > 0 || error < -7)
1727       return "unknown error";
1728    else 
1729       return error_strings[-error];
1730 }
1731