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