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