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