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