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