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