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