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