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