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