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