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