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