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