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