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