6ea83eb1850113809f83f6cf28e593d8114520e1
[opus.git] / libcelt / celt.c
1 /* Copyright (c) 2007-2008 CSIRO
2    Copyright (c) 2007-2010 Xiph.Org Foundation
3    Copyright (c) 2008 Gregory Maxwell 
4    Written by Jean-Marc Valin and Gregory Maxwell */
5 /*
6    Redistribution and use in source and binary forms, with or without
7    modification, are permitted provided that the following conditions
8    are met:
9    
10    - Redistributions of source code must retain the above copyright
11    notice, this list of conditions and the following disclaimer.
12    
13    - Redistributions in binary form must reproduce the above copyright
14    notice, this list of conditions and the following disclaimer in the
15    documentation and/or other materials provided with the distribution.
16    
17    - Neither the name of the Xiph.org Foundation nor the names of its
18    contributors may be used to endorse or promote products derived from
19    this software without specific prior written permission.
20    
21    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
22    ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
23    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
24    A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE FOUNDATION OR
25    CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
26    EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
27    PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
28    PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
29    LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
30    NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
31    SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
32 */
33
34 #ifdef HAVE_CONFIG_H
35 #include "config.h"
36 #endif
37
38 #define CELT_C
39
40 #include "os_support.h"
41 #include "mdct.h"
42 #include <math.h>
43 #include "celt.h"
44 #include "pitch.h"
45 #include "bands.h"
46 #include "modes.h"
47 #include "entcode.h"
48 #include "quant_bands.h"
49 #include "rate.h"
50 #include "stack_alloc.h"
51 #include "mathops.h"
52 #include "float_cast.h"
53 #include <stdarg.h>
54 #include "plc.h"
55
56 static const unsigned char trim_icdf[11] = {126, 124, 119, 109, 87, 41, 19, 9, 4, 2, 0};
57 /* Probs: NONE: 21.875%, LIGHT: 6.25%, NORMAL: 65.625%, AGGRESSIVE: 6.25% */
58 static const unsigned char spread_icdf[4] = {25, 23, 2, 0};
59
60 static const unsigned char tapset_icdf[3]={2,1,0};
61
62 #define COMBFILTER_MAXPERIOD 1024
63 #define COMBFILTER_MINPERIOD 16
64
65 /** Encoder state 
66  @brief Encoder state
67  */
68 struct CELTEncoder {
69    const CELTMode *mode;     /**< Mode used by the encoder */
70    int overlap;
71    int channels;
72    
73    int force_intra;
74    int complexity;
75    int start, end;
76
77    celt_int32 vbr_rate_norm; /* Target number of 8th bits per frame */
78    int constrained_vbr;      /* If zero, VBR can do whatever it likes with the rate */
79
80    /* Everything beyond this point gets cleared on a reset */
81 #define ENCODER_RESET_START frame_max
82
83    celt_word32 frame_max;
84    int spread_decision;
85    int delayedIntra;
86    int tonal_average;
87    int lastCodedBands;
88    int hf_average;
89    int tapset_decision;
90
91    int prefilter_period;
92    celt_word16 prefilter_gain;
93    int prefilter_tapset;
94 #ifdef RESYNTH
95    int prefilter_period_old;
96    celt_word16 prefilter_gain_old;
97    int prefilter_tapset_old;
98 #endif
99    int consec_transient;
100
101    /* VBR-related parameters */
102    celt_int32 vbr_reservoir;
103    celt_int32 vbr_drift;
104    celt_int32 vbr_offset;
105    celt_int32 vbr_count;
106
107    celt_word32 preemph_memE[2];
108    celt_word32 preemph_memD[2];
109
110 #ifdef RESYNTH
111    celt_sig syn_mem[2][2*MAX_PERIOD];
112 #endif
113
114    celt_sig in_mem[1]; /* Size = channels*mode->overlap */
115    /* celt_sig prefilter_mem[],  Size = channels*COMBFILTER_PERIOD */
116    /* celt_sig overlap_mem[],  Size = channels*mode->overlap */
117    /* celt_word16 oldEBands[], Size = 2*channels*mode->nbEBands */
118 };
119
120 int celt_encoder_get_size(const CELTMode *mode, int channels)
121 {
122    int size = sizeof(struct CELTEncoder)
123          + (2*channels*mode->overlap-1)*sizeof(celt_sig)
124          + channels*COMBFILTER_MAXPERIOD*sizeof(celt_sig)
125          + 2*channels*mode->nbEBands*sizeof(celt_word16);
126    return size;
127 }
128
129 CELTEncoder *celt_encoder_create(const CELTMode *mode, int channels, int *error)
130 {
131    return celt_encoder_init(
132          (CELTEncoder *)celt_alloc(celt_encoder_get_size(mode, channels)),
133          mode, channels, error);
134 }
135
136 CELTEncoder *celt_encoder_init(CELTEncoder *st, const CELTMode *mode, int channels, int *error)
137 {
138    if (channels < 0 || channels > 2)
139    {
140       if (error)
141          *error = CELT_BAD_ARG;
142       return NULL;
143    }
144
145    if (st==NULL)
146    {
147       if (error)
148          *error = CELT_ALLOC_FAIL;
149       return NULL;
150    }
151
152    CELT_MEMSET((char*)st, 0, celt_encoder_get_size(mode, channels));
153    
154    st->mode = mode;
155    st->overlap = mode->overlap;
156    st->channels = channels;
157
158    st->start = 0;
159    st->end = st->mode->effEBands;
160    st->constrained_vbr = 1;
161
162    st->vbr_rate_norm = 0;
163    st->vbr_offset = 0;
164    st->force_intra  = 0;
165    st->delayedIntra = 1;
166    st->tonal_average = 256;
167    st->spread_decision = SPREAD_NORMAL;
168    st->hf_average = 0;
169    st->tapset_decision = 0;
170    st->complexity = 5;
171
172    if (error)
173       *error = CELT_OK;
174    return st;
175 }
176
177 void celt_encoder_destroy(CELTEncoder *st)
178 {
179    celt_free(st);
180 }
181
182 static inline celt_int16 FLOAT2INT16(float x)
183 {
184    x = x*CELT_SIG_SCALE;
185    x = MAX32(x, -32768);
186    x = MIN32(x, 32767);
187    return (celt_int16)float2int(x);
188 }
189
190 static inline celt_word16 SIG2WORD16(celt_sig x)
191 {
192 #ifdef FIXED_POINT
193    x = PSHR32(x, SIG_SHIFT);
194    x = MAX32(x, -32768);
195    x = MIN32(x, 32767);
196    return EXTRACT16(x);
197 #else
198    return (celt_word16)x;
199 #endif
200 }
201
202 static int transient_analysis(const celt_word32 * restrict in, int len, int C,
203                               celt_word32 *frame_max, int overlap)
204 {
205    int i;
206    VARDECL(celt_word16, tmp);
207    celt_word32 mem0=0,mem1=0;
208    int is_transient = 0;
209    int block;
210    int N;
211    /* FIXME: Make that smaller */
212    celt_word16 bins[50];
213    SAVE_STACK;
214    ALLOC(tmp, len, celt_word16);
215
216    block = overlap/2;
217    N=len/block;
218    if (C==1)
219    {
220       for (i=0;i<len;i++)
221          tmp[i] = SHR32(in[i],SIG_SHIFT);
222    } else {
223       for (i=0;i<len;i++)
224          tmp[i] = SHR32(ADD32(in[i],in[i+len]), SIG_SHIFT+1);
225    }
226
227    /* High-pass filter: (1 - 2*z^-1 + z^-2) / (1 - z^-1 + .5*z^-2) */
228    for (i=0;i<len;i++)
229    {
230       celt_word32 x,y;
231       x = tmp[i];
232       y = ADD32(mem0, x);
233 #ifdef FIXED_POINT
234       mem0 = mem1 + y - SHL32(x,1);
235       mem1 = x - SHR32(y,1);
236 #else
237       mem0 = mem1 + y - 2*x;
238       mem1 = x - .5f*y;
239 #endif
240       tmp[i] = EXTRACT16(SHR(y,2));
241    }
242    /* First few samples are bad because we don't propagate the memory */
243    for (i=0;i<12;i++)
244       tmp[i] = 0;
245
246    for (i=0;i<N;i++)
247    {
248       int j;
249       float max_abs=0;
250       for (j=0;j<block;j++)
251          max_abs = MAX32(max_abs, tmp[i*block+j]);
252       bins[i] = max_abs;
253    }
254    for (i=0;i<N;i++)
255    {
256       int j;
257       int conseq=0;
258       celt_word16 t1, t2, t3;
259
260       t1 = MULT16_16_Q15(QCONST16(.15f, 15), bins[i]);
261       t2 = MULT16_16_Q15(QCONST16(.4f, 15), bins[i]);
262       t3 = MULT16_16_Q15(QCONST16(.15f, 15), bins[i]);
263       for (j=0;j<i;j++)
264       {
265          if (bins[j] < t1)
266             conseq++;
267          if (bins[j] < t2)
268             conseq++;
269          else
270             conseq = 0;
271       }
272       if (conseq>=3)
273          is_transient=1;
274       conseq = 0;
275       for (j=i+1;j<N;j++)
276       {
277          if (bins[j] < t3)
278             conseq++;
279          else
280             conseq = 0;
281       }
282       if (conseq>=7)
283          is_transient=1;
284    }
285    RESTORE_STACK;
286    return is_transient;
287 }
288
289 /** Apply window and compute the MDCT for all sub-frames and 
290     all channels in a frame */
291 static void compute_mdcts(const CELTMode *mode, int shortBlocks, celt_sig * restrict in, celt_sig * restrict out, int _C, int LM)
292 {
293    const int C = CHANNELS(_C);
294    if (C==1 && !shortBlocks)
295    {
296       const int overlap = OVERLAP(mode);
297       clt_mdct_forward(&mode->mdct, in, out, mode->window, overlap, mode->maxLM-LM);
298    } else {
299       const int overlap = OVERLAP(mode);
300       int N = mode->shortMdctSize<<LM;
301       int B = 1;
302       int b, c;
303       VARDECL(celt_word32, tmp);
304       SAVE_STACK;
305       if (shortBlocks)
306       {
307          /*lookup = &mode->mdct[0];*/
308          N = mode->shortMdctSize;
309          B = shortBlocks;
310       }
311       ALLOC(tmp, N, celt_word32);
312       c=0; do {
313          for (b=0;b<B;b++)
314          {
315             int j;
316             clt_mdct_forward(&mode->mdct, in+c*(B*N+overlap)+b*N, tmp, mode->window, overlap, shortBlocks ? mode->maxLM : mode->maxLM-LM);
317             /* Interleaving the sub-frames */
318             for (j=0;j<N;j++)
319                out[(j*B+b)+c*N*B] = tmp[j];
320          }
321       } while (++c<C);
322       RESTORE_STACK;
323    }
324 }
325
326 /** Compute the IMDCT and apply window for all sub-frames and 
327     all channels in a frame */
328 static void compute_inv_mdcts(const CELTMode *mode, int shortBlocks, celt_sig *X,
329       celt_sig * restrict out_mem[],
330       celt_sig * restrict overlap_mem[], int _C, int LM)
331 {
332    int c;
333    const int C = CHANNELS(_C);
334    const int N = mode->shortMdctSize<<LM;
335    const int overlap = OVERLAP(mode);
336    c=0; do {
337       int j;
338          VARDECL(celt_word32, x);
339          VARDECL(celt_word32, tmp);
340          int b;
341          int N2 = N;
342          int B = 1;
343          SAVE_STACK;
344          
345          ALLOC(x, N+overlap, celt_word32);
346          ALLOC(tmp, N, celt_word32);
347
348          if (shortBlocks)
349          {
350             N2 = mode->shortMdctSize;
351             B = shortBlocks;
352          }
353          /* Prevents problems from the imdct doing the overlap-add */
354          CELT_MEMSET(x, 0, overlap);
355
356          for (b=0;b<B;b++)
357          {
358             /* De-interleaving the sub-frames */
359             for (j=0;j<N2;j++)
360                tmp[j] = X[(j*B+b)+c*N2*B];
361             clt_mdct_backward(&mode->mdct, tmp, x+N2*b, mode->window, overlap, shortBlocks ? mode->maxLM : mode->maxLM-LM);
362          }
363
364          for (j=0;j<overlap;j++)
365             out_mem[c][j] = x[j] + overlap_mem[c][j];
366          for (;j<N;j++)
367             out_mem[c][j] = x[j];
368          for (j=0;j<overlap;j++)
369             overlap_mem[c][j] = x[N+j];
370          RESTORE_STACK;
371    } while (++c<C);
372 }
373
374 static void deemphasis(celt_sig *in[], celt_word16 *pcm, int N, int _C, const celt_word16 *coef, celt_sig *mem)
375 {
376    const int C = CHANNELS(_C);
377    int c;
378    c=0; do {
379       int j;
380       celt_sig * restrict x;
381       celt_word16  * restrict y;
382       celt_sig m = mem[c];
383       x =in[c];
384       y = pcm+c;
385       for (j=0;j<N;j++)
386       {
387          celt_sig tmp = *x + m;
388          m = MULT16_32_Q15(coef[0], tmp)
389            - MULT16_32_Q15(coef[1], *x);
390          tmp = SHL32(MULT16_32_Q15(coef[3], tmp), 2);
391          *y = SCALEOUT(SIG2WORD16(tmp));
392          x++;
393          y+=C;
394       }
395       mem[c] = m;
396    } while (++c<C);
397 }
398
399 #ifdef ENABLE_POSTFILTER
400 static void comb_filter(celt_word32 *y, celt_word32 *x, int T0, int T1, int N,
401       int C, celt_word16 g0, celt_word16 g1, int tapset0, int tapset1,
402       const celt_word16 *window, int overlap)
403 {
404    int i;
405    /* printf ("%d %d %f %f\n", T0, T1, g0, g1); */
406    celt_word16 g00, g01, g02, g10, g11, g12;
407    static const celt_word16 gains[3][3] = {
408          {QCONST16(0.30690f, 15), QCONST16(0.21701f, 15), QCONST16(0.12954f, 15)},
409          {QCONST16(0.46410f, 15), QCONST16(0.26795f, 15), QCONST16(0.f, 15)},
410          {QCONST16(0.8f, 15),     QCONST16(0.1f, 15),     QCONST16(0.f, 15)}};
411    g00 = MULT16_16_Q15(g0, gains[tapset0][0]);
412    g01 = MULT16_16_Q15(g0, gains[tapset0][1]);
413    g02 = MULT16_16_Q15(g0, gains[tapset0][2]);
414    g10 = MULT16_16_Q15(g1, gains[tapset1][0]);
415    g11 = MULT16_16_Q15(g1, gains[tapset1][1]);
416    g12 = MULT16_16_Q15(g1, gains[tapset1][2]);
417    for (i=0;i<overlap;i++)
418    {
419       celt_word16 f;
420       f = MULT16_16_Q15(window[i],window[i]);
421       y[i] = x[i]
422                + MULT16_32_Q15(MULT16_16_Q15((Q15ONE-f),g00),x[i-T0])
423                + MULT16_32_Q15(MULT16_16_Q15((Q15ONE-f),g01),x[i-T0-1])
424                + MULT16_32_Q15(MULT16_16_Q15((Q15ONE-f),g01),x[i-T0+1])
425                + MULT16_32_Q15(MULT16_16_Q15((Q15ONE-f),g02),x[i-T0-2])
426                + MULT16_32_Q15(MULT16_16_Q15((Q15ONE-f),g02),x[i-T0+2])
427                + MULT16_32_Q15(MULT16_16_Q15(f,g10),x[i-T1])
428                + MULT16_32_Q15(MULT16_16_Q15(f,g11),x[i-T1-1])
429                + MULT16_32_Q15(MULT16_16_Q15(f,g11),x[i-T1+1])
430                + MULT16_32_Q15(MULT16_16_Q15(f,g12),x[i-T1-2])
431                + MULT16_32_Q15(MULT16_16_Q15(f,g12),x[i-T1+2]);
432
433    }
434    for (i=overlap;i<N;i++)
435       y[i] = x[i]
436                + MULT16_32_Q15(g10,x[i-T1])
437                + MULT16_32_Q15(g11,x[i-T1-1])
438                + MULT16_32_Q15(g11,x[i-T1+1])
439                + MULT16_32_Q15(g12,x[i-T1-2])
440                + MULT16_32_Q15(g12,x[i-T1+2]);
441 }
442 #endif /* ENABLE_POSTFILTER */
443
444 static const signed char tf_select_table[4][8] = {
445       {0, -1, 0, -1,    0,-1, 0,-1},
446       {0, -1, 0, -2,    1, 0, 1,-1},
447       {0, -2, 0, -3,    2, 0, 1,-1},
448       {0, -2, 0, -3,    3, 0, 1,-1},
449 };
450
451 static celt_word32 l1_metric(const celt_norm *tmp, int N, int LM, int width)
452 {
453    int i, j;
454    static const celt_word16 sqrtM_1[4] = {Q15ONE, QCONST16(.70710678f,15), QCONST16(0.5f,15), QCONST16(0.35355339f,15)};
455    celt_word32 L1;
456    celt_word16 bias;
457    L1=0;
458    for (i=0;i<1<<LM;i++)
459    {
460       celt_word32 L2 = 0;
461       for (j=0;j<N>>LM;j++)
462          L2 = MAC16_16(L2, tmp[(j<<LM)+i], tmp[(j<<LM)+i]);
463       L1 += celt_sqrt(L2);
464    }
465    L1 = MULT16_32_Q15(sqrtM_1[LM], L1);
466    if (width==1)
467       bias = QCONST16(.12f,15)*LM;
468    else if (width==2)
469       bias = QCONST16(.05f,15)*LM;
470    else
471       bias = QCONST16(.02f,15)*LM;
472    L1 = MAC16_32_Q15(L1, bias, L1);
473    return L1;
474 }
475
476 static int tf_analysis(const CELTMode *m, celt_word16 *bandLogE, celt_word16 *oldBandE,
477       int len, int C, int isTransient, int *tf_res, int nbCompressedBytes, celt_norm *X,
478       int N0, int LM, int *tf_sum)
479 {
480    int i;
481    VARDECL(int, metric);
482    int cost0;
483    int cost1;
484    VARDECL(int, path0);
485    VARDECL(int, path1);
486    VARDECL(celt_norm, tmp);
487    int lambda;
488    int tf_select=0;
489    SAVE_STACK;
490
491    if (nbCompressedBytes<15*C)
492    {
493       *tf_sum = 0;
494       for (i=0;i<len;i++)
495          tf_res[i] = isTransient;
496       return 0;
497    }
498    if (nbCompressedBytes<40)
499       lambda = 12;
500    else if (nbCompressedBytes<60)
501       lambda = 6;
502    else if (nbCompressedBytes<100)
503       lambda = 4;
504    else
505       lambda = 3;
506
507    ALLOC(metric, len, int);
508    ALLOC(tmp, (m->eBands[len]-m->eBands[len-1])<<LM, celt_norm);
509    ALLOC(path0, len, int);
510    ALLOC(path1, len, int);
511
512    *tf_sum = 0;
513    for (i=0;i<len;i++)
514    {
515       int j, k, N;
516       celt_word32 L1, best_L1;
517       int best_level=0;
518       N = (m->eBands[i+1]-m->eBands[i])<<LM;
519       for (j=0;j<N;j++)
520          tmp[j] = X[j+(m->eBands[i]<<LM)];
521       /* Just add the right channel if we're in stereo */
522       if (C==2)
523          for (j=0;j<N;j++)
524             tmp[j] = ADD16(tmp[j],X[N0+j+(m->eBands[i]<<LM)]);
525       L1 = l1_metric(tmp, N, isTransient ? LM : 0, N>>LM);
526       best_L1 = L1;
527       /*printf ("%f ", L1);*/
528       for (k=0;k<LM;k++)
529       {
530          int B;
531
532          if (isTransient)
533             B = (LM-k-1);
534          else
535             B = k+1;
536
537          if (isTransient)
538             haar1(tmp, N>>(LM-k), 1<<(LM-k));
539          else
540             haar1(tmp, N>>k, 1<<k);
541
542          L1 = l1_metric(tmp, N, B, N>>LM);
543
544          if (L1 < best_L1)
545          {
546             best_L1 = L1;
547             best_level = k+1;
548          }
549       }
550       /*printf ("%d ", isTransient ? LM-best_level : best_level);*/
551       if (isTransient)
552          metric[i] = best_level;
553       else
554          metric[i] = -best_level;
555       *tf_sum += metric[i];
556    }
557    /*printf("\n");*/
558    /* FIXME: Figure out how to set this */
559    tf_select = 0;
560
561    cost0 = 0;
562    cost1 = isTransient ? 0 : lambda;
563    /* Viterbi forward pass */
564    for (i=1;i<len;i++)
565    {
566       int curr0, curr1;
567       int from0, from1;
568
569       from0 = cost0;
570       from1 = cost1 + lambda;
571       if (from0 < from1)
572       {
573          curr0 = from0;
574          path0[i]= 0;
575       } else {
576          curr0 = from1;
577          path0[i]= 1;
578       }
579
580       from0 = cost0 + lambda;
581       from1 = cost1;
582       if (from0 < from1)
583       {
584          curr1 = from0;
585          path1[i]= 0;
586       } else {
587          curr1 = from1;
588          path1[i]= 1;
589       }
590       cost0 = curr0 + abs(metric[i]-tf_select_table[LM][4*isTransient+2*tf_select+0]);
591       cost1 = curr1 + abs(metric[i]-tf_select_table[LM][4*isTransient+2*tf_select+1]);
592    }
593    tf_res[len-1] = cost0 < cost1 ? 0 : 1;
594    /* Viterbi backward pass to check the decisions */
595    for (i=len-2;i>=0;i--)
596    {
597       if (tf_res[i+1] == 1)
598          tf_res[i] = path1[i+1];
599       else
600          tf_res[i] = path0[i+1];
601    }
602    RESTORE_STACK;
603    return tf_select;
604 }
605
606 static void tf_encode(int start, int end, int isTransient, int *tf_res, int LM, int tf_select, ec_enc *enc)
607 {
608    int curr, i;
609    int tf_select_rsv;
610    int tf_changed;
611    int logp;
612    ec_uint32 budget;
613    ec_uint32 tell;
614    budget = enc->buf->storage*8;
615    tell = ec_enc_tell(enc, 0);
616    logp = isTransient ? 2 : 4;
617    /* Reserve space to code the tf_select decision. */
618    tf_select_rsv = LM>0 && tell+logp+1 <= budget;
619    budget -= tf_select_rsv;
620    curr = tf_changed = 0;
621    for (i=start;i<end;i++)
622    {
623       if (tell+logp<=budget)
624       {
625          ec_enc_bit_logp(enc, tf_res[i] ^ curr, logp);
626          tell = ec_enc_tell(enc, 0);
627          curr = tf_res[i];
628          tf_changed |= curr;
629       }
630       else
631          tf_res[i] = curr;
632       logp = isTransient ? 4 : 5;
633    }
634    /* Only code tf_select if it would actually make a difference. */
635    if (tf_select_rsv &&
636          tf_select_table[LM][4*isTransient+0+tf_changed]!=
637          tf_select_table[LM][4*isTransient+2+tf_changed])
638       ec_enc_bit_logp(enc, tf_select, 1);
639    else
640       tf_select = 0;
641    for (i=start;i<end;i++)
642       tf_res[i] = tf_select_table[LM][4*isTransient+2*tf_select+tf_res[i]];
643    /*printf("%d %d ", isTransient, tf_select); for(i=0;i<end;i++)printf("%d ", tf_res[i]);printf("\n");*/
644 }
645
646 static void tf_decode(int start, int end, int C, int isTransient, int *tf_res, int LM, ec_dec *dec)
647 {
648    int i, curr, tf_select;
649    int tf_select_rsv;
650    int tf_changed;
651    int logp;
652    ec_uint32 budget;
653    ec_uint32 tell;
654
655    budget = dec->buf->storage*8;
656    tell = ec_dec_tell(dec, 0);
657    logp = isTransient ? 2 : 4;
658    tf_select_rsv = LM>0 && tell+logp+1<=budget;
659    budget -= tf_select_rsv;
660    tf_changed = curr = 0;
661    for (i=start;i<end;i++)
662    {
663       if (tell+logp<=budget)
664       {
665          curr ^= ec_dec_bit_logp(dec, logp);
666          tell = ec_dec_tell(dec, 0);
667          tf_changed |= curr;
668       }
669       tf_res[i] = curr;
670       logp = isTransient ? 4 : 5;
671    }
672    tf_select = 0;
673    if (tf_select_rsv &&
674      tf_select_table[LM][4*isTransient+0+tf_changed] !=
675      tf_select_table[LM][4*isTransient+2+tf_changed])
676    {
677       tf_select = ec_dec_bit_logp(dec, 1);
678    }
679    for (i=start;i<end;i++)
680    {
681       tf_res[i] = tf_select_table[LM][4*isTransient+2*tf_select+tf_res[i]];
682    }
683 }
684
685 static int alloc_trim_analysis(const CELTMode *m, const celt_norm *X,
686       const celt_word16 *bandLogE, int nbEBands, int LM, int C, int N0)
687 {
688    int i;
689    celt_word32 diff=0;
690    int c;
691    int trim_index = 5;
692    if (C==2)
693    {
694       celt_word16 sum = 0; /* Q10 */
695       /* Compute inter-channel correlation for low frequencies */
696       for (i=0;i<8;i++)
697       {
698          int j;
699          celt_word32 partial = 0;
700          for (j=m->eBands[i]<<LM;j<m->eBands[i+1]<<LM;j++)
701             partial = MAC16_16(partial, X[j], X[N0+j]);
702          sum = ADD16(sum, EXTRACT16(SHR32(partial, 18)));
703       }
704       sum = MULT16_16_Q15(QCONST16(1.f/8, 15), sum);
705       /*printf ("%f\n", sum);*/
706       if (sum > QCONST16(.995f,10))
707          trim_index-=4;
708       else if (sum > QCONST16(.92f,10))
709          trim_index-=3;
710       else if (sum > QCONST16(.85f,10))
711          trim_index-=2;
712       else if (sum > QCONST16(.8f,10))
713          trim_index-=1;
714    }
715
716    /* Estimate spectral tilt */
717    c=0; do {
718       for (i=0;i<nbEBands-1;i++)
719       {
720          diff += bandLogE[i+c*nbEBands]*(celt_int32)(2+2*i-nbEBands);
721       }
722    } while (++c<0);
723    diff /= C*(nbEBands-1);
724    /*printf("%f\n", diff);*/
725    if (diff > QCONST16(2.f, DB_SHIFT))
726       trim_index--;
727    if (diff > QCONST16(8.f, DB_SHIFT))
728       trim_index--;
729    if (diff < -QCONST16(4.f, DB_SHIFT))
730       trim_index++;
731    if (diff < -QCONST16(10.f, DB_SHIFT))
732       trim_index++;
733
734    if (trim_index<0)
735       trim_index = 0;
736    if (trim_index>10)
737       trim_index = 10;
738    return trim_index;
739 }
740
741 static int stereo_analysis(const CELTMode *m, const celt_norm *X,
742       int nbEBands, int LM, int C, int N0)
743 {
744    int i;
745    int thetas;
746    celt_word32 sumLR = EPSILON, sumMS = EPSILON;
747
748    /* Use the L1 norm to model the entropy of the L/R signal vs the M/S signal */
749    for (i=0;i<13;i++)
750    {
751       int j;
752       for (j=m->eBands[i]<<LM;j<m->eBands[i+1]<<LM;j++)
753       {
754          celt_word16 L, R, M, S;
755          L = X[j];
756          R = X[N0+j];
757          M = L+R;
758          S = L-R;
759          sumLR += EXTEND32(ABS16(L)) + EXTEND32(ABS16(R));
760          sumMS += EXTEND32(ABS16(M)) + EXTEND32(ABS16(S));
761       }
762    }
763    sumMS = MULT16_32_Q15(QCONST16(0.707107f, 15), sumMS);
764    thetas = 13;
765    /* We don't need thetas for lower bands with LM<=1 */
766    if (LM<=1)
767       thetas -= 8;
768    return MULT16_32_Q15((m->eBands[13]<<(LM+1))+thetas, sumMS)
769          > MULT16_32_Q15(m->eBands[13]<<(LM+1), sumLR);
770 }
771
772 #ifdef FIXED_POINT
773 int celt_encode_with_ec(CELTEncoder * restrict st, const celt_int16 * pcm, int frame_size, unsigned char *compressed, int nbCompressedBytes, ec_enc *enc)
774 {
775 #else
776 int celt_encode_with_ec_float(CELTEncoder * restrict st, const celt_sig * pcm, int frame_size, unsigned char *compressed, int nbCompressedBytes, ec_enc *enc)
777 {
778 #endif
779    int i, c, N;
780    int bits;
781    ec_byte_buffer buf;
782    ec_enc         _enc;
783    VARDECL(celt_sig, in);
784    VARDECL(celt_sig, freq);
785    VARDECL(celt_norm, X);
786    VARDECL(celt_ener, bandE);
787    VARDECL(celt_word16, bandLogE);
788    VARDECL(celt_word16, oldLogE);
789    VARDECL(int, fine_quant);
790    VARDECL(celt_word16, error);
791    VARDECL(int, pulses);
792    VARDECL(int, offsets);
793    VARDECL(int, fine_priority);
794    VARDECL(int, tf_res);
795    celt_sig *_overlap_mem;
796    celt_sig *prefilter_mem;
797    celt_word16 *oldBandE, *oldLogE2;
798    int shortBlocks=0;
799    int isTransient=0;
800    int resynth;
801    const int C = CHANNELS(st->channels);
802    int LM, M;
803    int tf_select;
804    int nbFilledBytes, nbAvailableBytes;
805    int effEnd;
806    int codedBands;
807    int tf_sum;
808    int alloc_trim;
809    int pitch_index=COMBFILTER_MINPERIOD;
810    celt_word16 gain1 = 0;
811    int intensity=0;
812    int dual_stereo=0;
813    int effectiveBytes;
814    celt_word16 pf_threshold;
815    int dynalloc_logp;
816    celt_int32 vbr_rate;
817    celt_int32 total_bits;
818    celt_int32 total_boost;
819    celt_int32 tell;
820    int prefilter_tapset=0;
821    int pf_on;
822    int anti_collapse_on=0;
823    SAVE_STACK;
824
825    if (nbCompressedBytes<0 || pcm==NULL)
826      return CELT_BAD_ARG;
827
828    for (LM=0;LM<4;LM++)
829       if (st->mode->shortMdctSize<<LM==frame_size)
830          break;
831    if (LM>=MAX_CONFIG_SIZES)
832       return CELT_BAD_ARG;
833    M=1<<LM;
834
835    prefilter_mem = st->in_mem+C*(st->overlap);
836    _overlap_mem = prefilter_mem+C*COMBFILTER_MAXPERIOD;
837    /*_overlap_mem = st->in_mem+C*(st->overlap);*/
838    oldBandE = (celt_word16*)(st->in_mem+C*(2*st->overlap+COMBFILTER_MAXPERIOD));
839    oldLogE2 = oldBandE + C*st->mode->nbEBands;
840
841    if (enc==NULL)
842    {
843       ec_byte_writeinit_buffer(&buf, compressed, nbCompressedBytes);
844       ec_enc_init(&_enc,&buf);
845       enc = &_enc;
846       tell=1;
847       nbFilledBytes=0;
848    } else {
849       tell=ec_enc_tell(enc, 0);
850       nbFilledBytes=(tell+4)>>3;
851    }
852    nbAvailableBytes = nbCompressedBytes - nbFilledBytes;
853
854    vbr_rate = st->vbr_rate_norm<<LM;
855    if (vbr_rate>0)
856    {
857       effectiveBytes = st->vbr_rate_norm>>BITRES<<LM>>3;
858       /* Computes the max bit-rate allowed in VBR mode to avoid violating the
859           target rate and buffering.
860          We must do this up front so that bust-prevention logic triggers
861           correctly if we don't have enough bits. */
862       if (st->constrained_vbr)
863       {
864          celt_int32 vbr_bound;
865          celt_int32 max_allowed;
866          /* We could use any multiple of vbr_rate as bound (depending on the
867              delay).
868             This is clamped to ensure we use at least one byte if the encoder
869              was entirely empty, but to allow 0 in hybrid mode. */
870          vbr_bound = vbr_rate;
871          max_allowed = IMIN(IMAX((tell+7>>3)-nbFilledBytes,
872                vbr_rate+vbr_bound-st->vbr_reservoir>>(BITRES+3)),
873                nbAvailableBytes);
874          if(max_allowed < nbAvailableBytes)
875          {
876             nbCompressedBytes = nbFilledBytes+max_allowed;
877             nbAvailableBytes = max_allowed;
878             ec_byte_shrink(&buf, nbCompressedBytes);
879          }
880       }
881    } else
882       effectiveBytes = nbCompressedBytes;
883    total_bits = nbCompressedBytes*8;
884
885    effEnd = st->end;
886    if (effEnd > st->mode->effEBands)
887       effEnd = st->mode->effEBands;
888
889    N = M*st->mode->shortMdctSize;
890    ALLOC(in, C*(N+st->overlap), celt_sig);
891
892    /* Find pitch period and gain */
893    {
894       VARDECL(celt_sig, _pre);
895       celt_sig *pre[2];
896       SAVE_STACK;
897       c = 0;
898       ALLOC(_pre, C*(N+COMBFILTER_MAXPERIOD), celt_sig);
899
900       pre[0] = _pre;
901       pre[1] = _pre + (N+COMBFILTER_MAXPERIOD);
902
903       c=0; do {
904          const celt_word16 * restrict pcmp = pcm+c;
905          celt_sig * restrict inp = in+c*(N+st->overlap)+st->overlap;
906
907          for (i=0;i<N;i++)
908          {
909             /* Apply pre-emphasis */
910             celt_sig tmp = MULT16_16(st->mode->preemph[2], SCALEIN(*pcmp));
911             *inp = tmp + st->preemph_memE[c];
912             st->preemph_memE[c] = MULT16_32_Q15(st->mode->preemph[1], *inp)
913                                    - MULT16_32_Q15(st->mode->preemph[0], tmp);
914             inp++;
915             pcmp+=C;
916          }
917          CELT_COPY(pre[c], prefilter_mem+c*COMBFILTER_MAXPERIOD, COMBFILTER_MAXPERIOD);
918          CELT_COPY(pre[c]+COMBFILTER_MAXPERIOD, in+c*(N+st->overlap)+st->overlap, N);
919       } while (++c<C);
920
921 #ifdef ENABLE_POSTFILTER
922       if (nbAvailableBytes>12*C && st->start==0)
923       {
924          VARDECL(celt_word16, pitch_buf);
925          ALLOC(pitch_buf, (COMBFILTER_MAXPERIOD+N)>>1, celt_word16);
926          celt_word32 tmp=0;
927          celt_word32 mem0[2]={0,0};
928          celt_word16 mem1[2]={0,0};
929
930          pitch_downsample(pre, pitch_buf, COMBFILTER_MAXPERIOD+N, COMBFILTER_MAXPERIOD+N,
931                           C, mem0, mem1);
932          pitch_search(st->mode, pitch_buf+(COMBFILTER_MAXPERIOD>>1), pitch_buf, N,
933                COMBFILTER_MAXPERIOD-COMBFILTER_MINPERIOD, &pitch_index, &tmp, 1<<LM);
934          pitch_index = COMBFILTER_MAXPERIOD-pitch_index;
935
936          gain1 = remove_doubling(pitch_buf, COMBFILTER_MAXPERIOD, COMBFILTER_MINPERIOD,
937                N, &pitch_index, st->prefilter_period, st->prefilter_gain);
938          if (pitch_index > COMBFILTER_MAXPERIOD-2)
939             pitch_index = COMBFILTER_MAXPERIOD-2;
940          gain1 = MULT16_16_Q15(QCONST16(.7f,15),gain1);
941          prefilter_tapset = st->tapset_decision;
942       } else {
943          gain1 = 0;
944       }
945
946       /* Gain threshold for enabling the prefilter/postfilter */
947       pf_threshold = QCONST16(.2f,15);
948
949       /* Adjusting the threshold based on rate and continuity */
950       if (abs(pitch_index-st->prefilter_period)*10>pitch_index)
951          pf_threshold += QCONST16(.2f,15);
952       if (nbAvailableBytes<25)
953          pf_threshold += QCONST16(.1f,15);
954       if (nbAvailableBytes<35)
955          pf_threshold += QCONST16(.1f,15);
956       if (st->prefilter_gain > QCONST16(.4f,15))
957          pf_threshold -= QCONST16(.1f,15);
958       if (st->prefilter_gain > QCONST16(.55f,15))
959          pf_threshold -= QCONST16(.1f,15);
960
961       /* Hard threshold at 0.2 */
962       pf_threshold = MAX16(pf_threshold, QCONST16(.2f,15));
963       if (gain1<pf_threshold)
964       {
965          if(tell+15<=total_bits)
966             ec_enc_bit_logp(enc, 0, 1);
967          gain1 = 0;
968          pf_on = 0;
969       } else {
970          int qg;
971          int octave;
972
973          if (gain1 > QCONST16(.6f,15))
974             gain1 = QCONST16(.6f,15);
975          if (ABS16(gain1-st->prefilter_gain)<QCONST16(.1,15))
976             gain1=st->prefilter_gain;
977
978 #ifdef FIXED_POINT
979          qg = ((gain1+2048)>>12)-2;
980 #else
981          qg = floor(.5+gain1*8)-2;
982 #endif
983          ec_enc_bit_logp(enc, 1, 1);
984          pitch_index += 1;
985          octave = EC_ILOG(pitch_index)-5;
986          ec_enc_uint(enc, octave, 6);
987          ec_enc_bits(enc, pitch_index-(16<<octave), 4+octave);
988          pitch_index -= 1;
989          ec_enc_bits(enc, qg, 2);
990          gain1 = QCONST16(.125f,15)*(qg+2);
991          ec_enc_icdf(enc, prefilter_tapset, tapset_icdf, 2);
992          pf_on = 1;
993       }
994       /*printf("%d %f\n", pitch_index, gain1);*/
995 #else /* ENABLE_POSTFILTER */
996       if(tell+17<=total_bits)
997          ec_enc_bit_logp(enc, 0, 1);
998       pf_on = 0;
999 #endif /* ENABLE_POSTFILTER */
1000
1001       c=0; do {
1002          st->prefilter_period=IMAX(st->prefilter_period, COMBFILTER_MINPERIOD);
1003          CELT_COPY(in+c*(N+st->overlap), st->in_mem+c*(st->overlap), st->overlap);
1004 #ifdef ENABLE_POSTFILTER
1005          comb_filter(in+c*(N+st->overlap)+st->overlap, pre[c]+COMBFILTER_MAXPERIOD,
1006                st->prefilter_period, pitch_index, N, C, -st->prefilter_gain, -gain1,
1007                st->prefilter_tapset, prefilter_tapset, st->mode->window, st->mode->overlap);
1008 #endif /* ENABLE_POSTFILTER */
1009          CELT_COPY(st->in_mem+c*(st->overlap), in+c*(N+st->overlap)+N, st->overlap);
1010
1011 #ifdef ENABLE_POSTFILTER
1012          if (N>COMBFILTER_MAXPERIOD)
1013          {
1014             CELT_MOVE(prefilter_mem+c*COMBFILTER_MAXPERIOD, pre[c]+N, COMBFILTER_MAXPERIOD);
1015          } else {
1016             CELT_MOVE(prefilter_mem+c*COMBFILTER_MAXPERIOD, prefilter_mem+c*COMBFILTER_MAXPERIOD+N, COMBFILTER_MAXPERIOD-N);
1017             CELT_MOVE(prefilter_mem+c*COMBFILTER_MAXPERIOD+COMBFILTER_MAXPERIOD-N, pre[c]+COMBFILTER_MAXPERIOD, N);
1018          }
1019 #endif /* ENABLE_POSTFILTER */
1020       } while (++c<C);
1021
1022       RESTORE_STACK;
1023    }
1024
1025 #ifdef RESYNTH
1026    resynth = 1;
1027 #else
1028    resynth = 0;
1029 #endif
1030
1031    isTransient = 0;
1032    shortBlocks = 0;
1033    if (LM>0 && ec_enc_tell(enc, 0)+3<=total_bits)
1034    {
1035       if (st->complexity > 1)
1036       {
1037          isTransient = transient_analysis(in, N+st->overlap, C,
1038                   &st->frame_max, st->overlap);
1039          if (isTransient)
1040             shortBlocks = M;
1041       }
1042       ec_enc_bit_logp(enc, isTransient, 3);
1043    }
1044
1045    ALLOC(freq, C*N, celt_sig); /**< Interleaved signal MDCTs */
1046    ALLOC(bandE,st->mode->nbEBands*C, celt_ener);
1047    ALLOC(bandLogE,st->mode->nbEBands*C, celt_word16);
1048    /* Compute MDCTs */
1049    compute_mdcts(st->mode, shortBlocks, in, freq, C, LM);
1050
1051    ALLOC(X, C*N, celt_norm);         /**< Interleaved normalised MDCTs */
1052
1053    compute_band_energies(st->mode, freq, bandE, effEnd, C, M);
1054
1055    amp2Log2(st->mode, effEnd, st->end, bandE, bandLogE, C);
1056
1057    /* Band normalisation */
1058    normalise_bands(st->mode, freq, X, bandE, effEnd, C, M);
1059
1060    ALLOC(tf_res, st->mode->nbEBands, int);
1061    /* Needs to be before coarse energy quantization because otherwise the energy gets modified */
1062    tf_select = tf_analysis(st->mode, bandLogE, oldBandE, effEnd, C, isTransient, tf_res, effectiveBytes, X, N, LM, &tf_sum);
1063    for (i=effEnd;i<st->end;i++)
1064       tf_res[i] = tf_res[effEnd-1];
1065
1066    ALLOC(oldLogE, C*st->mode->nbEBands, celt_word16);
1067    for (i=0;i<C*st->mode->nbEBands;i++)
1068       oldLogE[i] = oldBandE[i];
1069    ALLOC(error, C*st->mode->nbEBands, celt_word16);
1070    quant_coarse_energy(st->mode, st->start, st->end, effEnd, bandLogE,
1071          oldBandE, total_bits, error, enc,
1072          C, LM, nbAvailableBytes, st->force_intra,
1073          &st->delayedIntra, st->complexity >= 4);
1074
1075    tf_encode(st->start, st->end, isTransient, tf_res, LM, tf_select, enc);
1076
1077    st->spread_decision = SPREAD_NORMAL;
1078    if (ec_enc_tell(enc, 0)+4<=total_bits)
1079    {
1080       if (shortBlocks || st->complexity < 3 || nbAvailableBytes < 10*C)
1081       {
1082          if (st->complexity == 0)
1083             st->spread_decision = SPREAD_NONE;
1084       } else {
1085          st->spread_decision = spreading_decision(st->mode, X,
1086                &st->tonal_average, st->spread_decision, &st->hf_average,
1087                &st->tapset_decision, pf_on&&!shortBlocks, effEnd, C, M);
1088       }
1089       ec_enc_icdf(enc, st->spread_decision, spread_icdf, 5);
1090    }
1091
1092    ALLOC(offsets, st->mode->nbEBands, int);
1093
1094    for (i=0;i<st->mode->nbEBands;i++)
1095       offsets[i] = 0;
1096    /* Dynamic allocation code */
1097    /* Make sure that dynamic allocation can't make us bust the budget */
1098    if (effectiveBytes > 50 && LM>=1)
1099    {
1100       int t1, t2;
1101       if (LM <= 1)
1102       {
1103          t1 = 3;
1104          t2 = 5;
1105       } else {
1106          t1 = 2;
1107          t2 = 4;
1108       }
1109       for (i=1;i<st->mode->nbEBands-1;i++)
1110       {
1111          celt_word32 d2;
1112          d2 = 2*bandLogE[i]-bandLogE[i-1]-bandLogE[i+1];
1113          if (C==2)
1114             d2 = HALF32(d2 + 2*bandLogE[i+st->mode->nbEBands]-
1115                   bandLogE[i-1+st->mode->nbEBands]-bandLogE[i+1+st->mode->nbEBands]);
1116          if (d2 > SHL16(t1,DB_SHIFT))
1117             offsets[i] += 1;
1118          if (d2 > SHL16(t2,DB_SHIFT))
1119             offsets[i] += 1;
1120       }
1121    }
1122    dynalloc_logp = 6;
1123    total_bits<<=BITRES;
1124    total_boost = 0;
1125    tell = ec_enc_tell(enc, BITRES);
1126    for (i=st->start;i<st->end;i++)
1127    {
1128       int width, quanta;
1129       int dynalloc_loop_logp;
1130       int boost;
1131       int j;
1132       width = C*(st->mode->eBands[i+1]-st->mode->eBands[i])<<LM;
1133       /* quanta is 6 bits, but no more than 1 bit/sample
1134          and no less than 1/8 bit/sample */
1135       quanta = IMIN(width<<BITRES, IMAX(6<<BITRES, width));
1136       dynalloc_loop_logp = dynalloc_logp;
1137       boost = 0;
1138       for (j = 0; tell+(dynalloc_loop_logp<<BITRES) < total_bits-total_boost
1139             && boost < (64<<LM)*(C<<BITRES); j++)
1140       {
1141          int flag;
1142          flag = j<offsets[i];
1143          ec_enc_bit_logp(enc, flag, dynalloc_loop_logp);
1144          tell = ec_enc_tell(enc, BITRES);
1145          if (!flag)
1146             break;
1147          boost += quanta;
1148          total_boost += quanta;
1149          dynalloc_loop_logp = 1;
1150       }
1151       /* Making dynalloc more likely */
1152       if (j)
1153          dynalloc_logp = IMAX(2, dynalloc_logp-1);
1154       offsets[i] = boost;
1155    }
1156    alloc_trim = 5;
1157    if (tell+(6<<BITRES) <= total_bits - total_boost)
1158    {
1159       alloc_trim = alloc_trim_analysis(st->mode, X, bandLogE,
1160             st->mode->nbEBands, LM, C, N);
1161       ec_enc_icdf(enc, alloc_trim, trim_icdf, 7);
1162       tell = ec_enc_tell(enc, BITRES);
1163    }
1164
1165    /* Variable bitrate */
1166    if (vbr_rate>0)
1167    {
1168      celt_word16 alpha;
1169      celt_int32 delta;
1170      /* The target rate in 8th bits per frame */
1171      celt_int32 target;
1172      celt_int32 min_allowed;
1173
1174      target = vbr_rate + st->vbr_offset - ((40*C+20)<<BITRES);
1175
1176      /* Shortblocks get a large boost in bitrate, but since they
1177         are uncommon long blocks are not greatly affected */
1178      if (shortBlocks || tf_sum < -2*(st->end-st->start))
1179         target = 7*target/4;
1180      else if (tf_sum < -(st->end-st->start))
1181         target = 3*target/2;
1182      else if (M > 1)
1183         target-=(target+14)/28;
1184
1185      /* The current offset is removed from the target and the space used
1186         so far is added*/
1187      target=target+tell;
1188      /* By how much did we "miss" the target on that frame */
1189      delta = target - vbr_rate;
1190
1191      /* In VBR mode the frame size must not be reduced so much that it would
1192          result in the encoder running out of bits.
1193         The margin of 2 bytes ensures that none of the bust-prevention logic
1194          in the decoder will have triggered so far. */
1195      min_allowed = (tell+total_boost+(1<<BITRES+3)-1>>(BITRES+3)) + 2 - nbFilledBytes;
1196
1197      nbAvailableBytes = target+(1<<(BITRES+2))>>(BITRES+3);
1198      nbAvailableBytes = IMAX(min_allowed,nbAvailableBytes);
1199      nbAvailableBytes = IMIN(nbCompressedBytes,nbAvailableBytes+nbFilledBytes) - nbFilledBytes;
1200
1201      target=nbAvailableBytes<<(BITRES+3);
1202
1203      if (st->vbr_count < 970)
1204      {
1205         st->vbr_count++;
1206         alpha = celt_rcp(SHL32(EXTEND32(st->vbr_count+20),16));
1207      } else
1208         alpha = QCONST16(.001f,15);
1209      /* How many bits have we used in excess of what we're allowed */
1210      if (st->constrained_vbr)
1211         st->vbr_reservoir += target - vbr_rate;
1212      /*printf ("%d\n", st->vbr_reservoir);*/
1213
1214      /* Compute the offset we need to apply in order to reach the target */
1215      st->vbr_drift += (celt_int32)MULT16_32_Q15(alpha,delta-st->vbr_offset-st->vbr_drift);
1216      st->vbr_offset = -st->vbr_drift;
1217      /*printf ("%d\n", st->vbr_drift);*/
1218
1219      if (st->constrained_vbr && st->vbr_reservoir < 0)
1220      {
1221         /* We're under the min value -- increase rate */
1222         int adjust = (-st->vbr_reservoir)/(8<<BITRES);
1223         nbAvailableBytes += adjust;
1224         st->vbr_reservoir = 0;
1225         /*printf ("+%d\n", adjust);*/
1226      }
1227      nbCompressedBytes = IMIN(nbCompressedBytes,nbAvailableBytes+nbFilledBytes);
1228
1229      /* This moves the raw bits to take into account the new compressed size */
1230      ec_byte_shrink(&buf, nbCompressedBytes);
1231    }
1232
1233    if (C==2)
1234    {
1235       int effectiveRate;
1236
1237       /* Always use MS for 2.5 ms frames until we can do a better analysis */
1238       if (LM!=0)
1239          dual_stereo = stereo_analysis(st->mode, X, st->mode->nbEBands, LM, C, N);
1240
1241       /* Account for coarse energy */
1242       effectiveRate = (8*effectiveBytes - 80)>>LM;
1243
1244       /* effectiveRate in kb/s */
1245       effectiveRate = 2*effectiveRate/5;
1246       if (effectiveRate<35)
1247          intensity = 8;
1248       else if (effectiveRate<50)
1249          intensity = 12;
1250       else if (effectiveRate<68)
1251          intensity = 16;
1252       else if (effectiveRate<84)
1253          intensity = 18;
1254       else if (effectiveRate<102)
1255          intensity = 19;
1256       else if (effectiveRate<130)
1257          intensity = 20;
1258       else
1259          intensity = 100;
1260       intensity = IMIN(st->end,IMAX(st->start, intensity));
1261    }
1262
1263    /* Bit allocation */
1264    ALLOC(fine_quant, st->mode->nbEBands, int);
1265    ALLOC(pulses, st->mode->nbEBands, int);
1266    ALLOC(fine_priority, st->mode->nbEBands, int);
1267
1268    /* bits =   packet size        -       where we are         - safety -  anti-collapse*/
1269    bits = (nbCompressedBytes*8<<BITRES) - ec_enc_tell(enc, BITRES) - 1 - (isTransient&&LM>=2 ? (1<<BITRES) : 0);
1270    codedBands = compute_allocation(st->mode, st->start, st->end, offsets,
1271          alloc_trim, &intensity, &dual_stereo, bits, pulses, fine_quant,
1272          fine_priority, C, LM, enc, 1, st->lastCodedBands);
1273    st->lastCodedBands = codedBands;
1274
1275    quant_fine_energy(st->mode, st->start, st->end, bandE, oldBandE, error, fine_quant, enc, C);
1276
1277 #ifdef MEASURE_NORM_MSE
1278    float X0[3000];
1279    float bandE0[60];
1280    c=0; do 
1281       for (i=0;i<N;i++)
1282          X0[i+c*N] = X[i+c*N];
1283    while (++c<C);
1284    for (i=0;i<C*st->mode->nbEBands;i++)
1285       bandE0[i] = bandE[i];
1286 #endif
1287
1288    /* Residual quantisation */
1289    quant_all_bands(1, st->mode, st->start, st->end, X, C==2 ? X+N : NULL,
1290          bandE, pulses, shortBlocks, st->spread_decision, dual_stereo, intensity, tf_res, resynth,
1291          nbCompressedBytes*8, enc, LM, codedBands);
1292
1293    if (isTransient && LM>=2)
1294    {
1295       anti_collapse_on = st->consec_transient<2;
1296       ec_enc_bits(enc, anti_collapse_on, 1);
1297    }
1298    quant_energy_finalise(st->mode, st->start, st->end, bandE, oldBandE, error, fine_quant, fine_priority, nbCompressedBytes*8-ec_enc_tell(enc, 0), enc, C);
1299
1300 #ifdef RESYNTH
1301    /* Re-synthesis of the coded audio if required */
1302    if (resynth)
1303    {
1304       celt_sig *out_mem[2];
1305       celt_sig *overlap_mem[2];
1306
1307       log2Amp(st->mode, st->start, st->end, bandE, oldBandE, C);
1308
1309 #ifdef MEASURE_NORM_MSE
1310       measure_norm_mse(st->mode, X, X0, bandE, bandE0, M, N, C);
1311 #endif
1312       if (anti_collapse_on)
1313       {
1314          anti_collapse(st->mode, X, LM, C, N,
1315                st->start, st->end, oldBandE, oldLogE, oldLogE2, pulses, enc->rng);
1316       }
1317
1318       /* Synthesis */
1319       denormalise_bands(st->mode, X, freq, bandE, effEnd, C, M);
1320
1321       CELT_MOVE(st->syn_mem[0], st->syn_mem[0]+N, MAX_PERIOD);
1322       if (C==2)
1323          CELT_MOVE(st->syn_mem[1], st->syn_mem[1]+N, MAX_PERIOD);
1324
1325       c=0; do
1326          for (i=0;i<M*st->mode->eBands[st->start];i++)
1327             freq[c*N+i] = 0;
1328       while (++c<C);
1329       c=0; do
1330          for (i=M*st->mode->eBands[st->end];i<N;i++)
1331             freq[c*N+i] = 0;
1332       while (++c<C);
1333
1334       out_mem[0] = st->syn_mem[0]+MAX_PERIOD;
1335       if (C==2)
1336          out_mem[1] = st->syn_mem[1]+MAX_PERIOD;
1337
1338       c=0; do
1339          overlap_mem[c] = _overlap_mem + c*st->overlap;
1340       while (++c<C);
1341
1342       compute_inv_mdcts(st->mode, shortBlocks, freq, out_mem, overlap_mem, C, LM);
1343
1344 #ifdef ENABLE_POSTFILTER
1345       c=0; do {
1346          st->prefilter_period=IMAX(st->prefilter_period, COMBFILTER_MINPERIOD);
1347          st->prefilter_period_old=IMAX(st->prefilter_period_old, COMBFILTER_MINPERIOD);
1348          if (LM!=0)
1349          {
1350             comb_filter(out_mem[c], out_mem[c], st->prefilter_period, st->prefilter_period, st->overlap, C,
1351                   st->prefilter_gain, st->prefilter_gain, st->prefilter_tapset, st->prefilter_tapset,
1352                   NULL, 0);
1353             comb_filter(out_mem[c]+st->overlap, out_mem[c]+st->overlap, st->prefilter_period, pitch_index, N-st->overlap, C,
1354                   st->prefilter_gain, gain1, st->prefilter_tapset, prefilter_tapset,
1355                   st->mode->window, st->mode->overlap);
1356          } else {
1357             comb_filter(out_mem[c], out_mem[c], st->prefilter_period_old, st->prefilter_period, N, C,
1358                   st->prefilter_gain_old, st->prefilter_gain, st->prefilter_tapset_old, st->prefilter_tapset,
1359                   st->mode->window, st->mode->overlap);
1360          }
1361       } while (++c<C);
1362 #endif /* ENABLE_POSTFILTER */
1363
1364       deemphasis(out_mem, (celt_word16*)pcm, N, C, st->mode->preemph, st->preemph_memD);
1365       st->prefilter_period_old = st->prefilter_period;
1366       st->prefilter_gain_old = st->prefilter_gain;
1367       st->prefilter_tapset_old = st->prefilter_tapset;
1368    }
1369 #endif
1370
1371    st->prefilter_period = pitch_index;
1372    st->prefilter_gain = gain1;
1373    st->prefilter_tapset = prefilter_tapset;
1374
1375    /* In case start or end were to change */
1376    for (i=0;i<st->start;i++)
1377       oldBandE[i]=0;
1378    for (i=st->end;i<st->mode->nbEBands;i++)
1379       oldBandE[i]=0;
1380    for (i=0;i<C*st->mode->nbEBands;i++)
1381       oldLogE2[i] = oldLogE[i];
1382    if (isTransient)
1383       st->consec_transient++;
1384    else
1385       st->consec_transient=0;
1386
1387    /* If there's any room left (can only happen for very high rates),
1388       fill it with zeros */
1389    while (ec_enc_tell(enc,0) + 8 <= nbCompressedBytes*8)
1390       ec_enc_bits(enc, 0, 8);
1391    ec_enc_done(enc);
1392    
1393    RESTORE_STACK;
1394    if (ec_enc_get_error(enc))
1395       return CELT_CORRUPTED_DATA;
1396    else
1397       return nbCompressedBytes;
1398 }
1399
1400 #ifdef FIXED_POINT
1401 #ifndef DISABLE_FLOAT_API
1402 int celt_encode_with_ec_float(CELTEncoder * restrict st, const float * pcm, int frame_size, unsigned char *compressed, int nbCompressedBytes, ec_enc *enc)
1403 {
1404    int j, ret, C, N, LM, M;
1405    VARDECL(celt_int16, in);
1406    SAVE_STACK;
1407
1408    if (pcm==NULL)
1409       return CELT_BAD_ARG;
1410
1411    for (LM=0;LM<4;LM++)
1412       if (st->mode->shortMdctSize<<LM==frame_size)
1413          break;
1414    if (LM>=MAX_CONFIG_SIZES)
1415       return CELT_BAD_ARG;
1416    M=1<<LM;
1417
1418    C = CHANNELS(st->channels);
1419    N = M*st->mode->shortMdctSize;
1420    ALLOC(in, C*N, celt_int16);
1421
1422    for (j=0;j<C*N;j++)
1423      in[j] = FLOAT2INT16(pcm[j]);
1424
1425    ret=celt_encode_with_ec(st,in,frame_size,compressed,nbCompressedBytes, enc);
1426 #ifdef RESYNTH
1427    for (j=0;j<C*N;j++)
1428       ((float*)pcm)[j]=in[j]*(1.f/32768.f);
1429 #endif
1430    RESTORE_STACK;
1431    return ret;
1432
1433 }
1434 #endif /*DISABLE_FLOAT_API*/
1435 #else
1436 int celt_encode_with_ec(CELTEncoder * restrict st, const celt_int16 * pcm, int frame_size, unsigned char *compressed, int nbCompressedBytes, ec_enc *enc)
1437 {
1438    int j, ret, C, N, LM, M;
1439    VARDECL(celt_sig, in);
1440    SAVE_STACK;
1441
1442    if (pcm==NULL)
1443       return CELT_BAD_ARG;
1444
1445    for (LM=0;LM<4;LM++)
1446       if (st->mode->shortMdctSize<<LM==frame_size)
1447          break;
1448    if (LM>=MAX_CONFIG_SIZES)
1449       return CELT_BAD_ARG;
1450    M=1<<LM;
1451
1452    C=CHANNELS(st->channels);
1453    N=M*st->mode->shortMdctSize;
1454    ALLOC(in, C*N, celt_sig);
1455    for (j=0;j<C*N;j++) {
1456      in[j] = SCALEOUT(pcm[j]);
1457    }
1458
1459    ret = celt_encode_with_ec_float(st,in,frame_size,compressed,nbCompressedBytes, enc);
1460 #ifdef RESYNTH
1461    for (j=0;j<C*N;j++)
1462       ((celt_int16*)pcm)[j] = FLOAT2INT16(in[j]);
1463 #endif
1464    RESTORE_STACK;
1465    return ret;
1466 }
1467 #endif
1468
1469 int celt_encode(CELTEncoder * restrict st, const celt_int16 * pcm, int frame_size, unsigned char *compressed, int nbCompressedBytes)
1470 {
1471    return celt_encode_with_ec(st, pcm, frame_size, compressed, nbCompressedBytes, NULL);
1472 }
1473
1474 #ifndef DISABLE_FLOAT_API
1475 int celt_encode_float(CELTEncoder * restrict st, const float * pcm, int frame_size, unsigned char *compressed, int nbCompressedBytes)
1476 {
1477    return celt_encode_with_ec_float(st, pcm, frame_size, compressed, nbCompressedBytes, NULL);
1478 }
1479 #endif /* DISABLE_FLOAT_API */
1480
1481 int celt_encoder_ctl(CELTEncoder * restrict st, int request, ...)
1482 {
1483    va_list ap;
1484    
1485    va_start(ap, request);
1486    switch (request)
1487    {
1488       case CELT_GET_MODE_REQUEST:
1489       {
1490          const CELTMode ** value = va_arg(ap, const CELTMode**);
1491          if (value==0)
1492             goto bad_arg;
1493          *value=st->mode;
1494       }
1495       break;
1496       case CELT_SET_COMPLEXITY_REQUEST:
1497       {
1498          int value = va_arg(ap, celt_int32);
1499          if (value<0 || value>10)
1500             goto bad_arg;
1501          st->complexity = value;
1502       }
1503       break;
1504       case CELT_SET_START_BAND_REQUEST:
1505       {
1506          celt_int32 value = va_arg(ap, celt_int32);
1507          if (value<0 || value>=st->mode->nbEBands)
1508             goto bad_arg;
1509          st->start = value;
1510       }
1511       break;
1512       case CELT_SET_END_BAND_REQUEST:
1513       {
1514          celt_int32 value = va_arg(ap, celt_int32);
1515          if (value<1 || value>st->mode->nbEBands)
1516             goto bad_arg;
1517          st->end = value;
1518       }
1519       break;
1520       case CELT_SET_PREDICTION_REQUEST:
1521       {
1522          int value = va_arg(ap, celt_int32);
1523          if (value<0 || value>2)
1524             goto bad_arg;
1525          if (value==0)
1526          {
1527             st->force_intra   = 1;
1528          } else if (value==1) {
1529             st->force_intra   = 0;
1530          } else {
1531             st->force_intra   = 0;
1532          }   
1533       }
1534       break;
1535       case CELT_SET_VBR_CONSTRAINT_REQUEST:
1536       {
1537          celt_int32 value = va_arg(ap, celt_int32);
1538          st->constrained_vbr = value;
1539       }
1540       break;
1541       case CELT_SET_VBR_RATE_REQUEST:
1542       {
1543          celt_int32 value = va_arg(ap, celt_int32);
1544          int frame_rate;
1545          int N = st->mode->shortMdctSize;
1546          if (value<0)
1547             goto bad_arg;
1548          if (value>3072000)
1549             value = 3072000;
1550          frame_rate = ((st->mode->Fs<<3)+(N>>1))/N;
1551          st->vbr_rate_norm = ((value<<(BITRES+3))+(frame_rate>>1))/frame_rate;
1552       }
1553       break;
1554       case CELT_RESET_STATE:
1555       {
1556          CELT_MEMSET((char*)&st->ENCODER_RESET_START, 0,
1557                celt_encoder_get_size(st->mode, st->channels)-
1558                ((char*)&st->ENCODER_RESET_START - (char*)st));
1559          st->vbr_offset = 0;
1560          st->delayedIntra = 1;
1561          st->spread_decision = SPREAD_NORMAL;
1562          st->tonal_average = QCONST16(1.f,8);
1563       }
1564       break;
1565       default:
1566          goto bad_request;
1567    }
1568    va_end(ap);
1569    return CELT_OK;
1570 bad_arg:
1571    va_end(ap);
1572    return CELT_BAD_ARG;
1573 bad_request:
1574    va_end(ap);
1575    return CELT_UNIMPLEMENTED;
1576 }
1577
1578 /**********************************************************************/
1579 /*                                                                    */
1580 /*                             DECODER                                */
1581 /*                                                                    */
1582 /**********************************************************************/
1583 #define DECODE_BUFFER_SIZE 2048
1584
1585 /** Decoder state 
1586  @brief Decoder state
1587  */
1588 struct CELTDecoder {
1589    const CELTMode *mode;
1590    int overlap;
1591    int channels;
1592
1593    int start, end;
1594
1595    /* Everything beyond this point gets cleared on a reset */
1596 #define DECODER_RESET_START last_pitch_index
1597
1598    int last_pitch_index;
1599    int loss_count;
1600    int postfilter_period;
1601    int postfilter_period_old;
1602    celt_word16 postfilter_gain;
1603    celt_word16 postfilter_gain_old;
1604    int postfilter_tapset;
1605    int postfilter_tapset_old;
1606
1607    celt_sig preemph_memD[2];
1608    
1609    celt_sig _decode_mem[1]; /* Size = channels*(DECODE_BUFFER_SIZE+mode->overlap) */
1610    /* celt_word16 lpc[],  Size = channels*LPC_ORDER */
1611    /* celt_word16 oldEBands[], Size = 2*channels*mode->nbEBands */
1612 };
1613
1614 int celt_decoder_get_size(const CELTMode *mode, int channels)
1615 {
1616    int size = sizeof(struct CELTDecoder)
1617             + (channels*(DECODE_BUFFER_SIZE+mode->overlap)-1)*sizeof(celt_sig)
1618             + channels*LPC_ORDER*sizeof(celt_word16)
1619             + 2*channels*mode->nbEBands*sizeof(celt_word16);
1620    return size;
1621 }
1622
1623 CELTDecoder *celt_decoder_create(const CELTMode *mode, int channels, int *error)
1624 {
1625    return celt_decoder_init(
1626          (CELTDecoder *)celt_alloc(celt_decoder_get_size(mode, channels)),
1627          mode, channels, error);
1628 }
1629
1630 CELTDecoder *celt_decoder_init(CELTDecoder *st, const CELTMode *mode, int channels, int *error)
1631 {
1632    if (channels < 0 || channels > 2)
1633    {
1634       if (error)
1635          *error = CELT_BAD_ARG;
1636       return NULL;
1637    }
1638
1639    if (st==NULL)
1640    {
1641       if (error)
1642          *error = CELT_ALLOC_FAIL;
1643       return NULL;
1644    }
1645
1646    CELT_MEMSET((char*)st, 0, celt_decoder_get_size(mode, channels));
1647
1648    st->mode = mode;
1649    st->overlap = mode->overlap;
1650    st->channels = channels;
1651
1652    st->start = 0;
1653    st->end = st->mode->effEBands;
1654
1655    st->loss_count = 0;
1656
1657    if (error)
1658       *error = CELT_OK;
1659    return st;
1660 }
1661
1662 void celt_decoder_destroy(CELTDecoder *st)
1663 {
1664    celt_free(st);
1665 }
1666
1667 static void celt_decode_lost(CELTDecoder * restrict st, celt_word16 * restrict pcm, int N, int LM)
1668 {
1669    int c;
1670    int pitch_index;
1671    int overlap = st->mode->overlap;
1672    celt_word16 fade = Q15ONE;
1673    int i, len;
1674    const int C = CHANNELS(st->channels);
1675    int offset;
1676    celt_sig *out_mem[2];
1677    celt_sig *decode_mem[2];
1678    celt_sig *overlap_mem[2];
1679    celt_word16 *lpc;
1680    SAVE_STACK;
1681    
1682    c=0; do {
1683       decode_mem[c] = st->_decode_mem + c*(DECODE_BUFFER_SIZE+st->overlap);
1684       out_mem[c] = decode_mem[c]+DECODE_BUFFER_SIZE-MAX_PERIOD;
1685       overlap_mem[c] = decode_mem[c]+DECODE_BUFFER_SIZE;
1686    } while (++c<C);
1687    lpc = (celt_word16*)(st->_decode_mem+(DECODE_BUFFER_SIZE+st->overlap)*C);
1688
1689    len = N+st->mode->overlap;
1690    
1691    if (st->loss_count == 0)
1692    {
1693       celt_word16 pitch_buf[MAX_PERIOD>>1];
1694       celt_word32 tmp=0;
1695       celt_word32 mem0[2]={0,0};
1696       celt_word16 mem1[2]={0,0};
1697       int len2 = len;
1698       /* FIXME: This is a kludge */
1699       if (len2>MAX_PERIOD>>1)
1700          len2 = MAX_PERIOD>>1;
1701       pitch_downsample(out_mem, pitch_buf, MAX_PERIOD, MAX_PERIOD,
1702                        C, mem0, mem1);
1703       pitch_search(st->mode, pitch_buf+((MAX_PERIOD-len2)>>1), pitch_buf, len2,
1704                    MAX_PERIOD-len2-100, &pitch_index, &tmp, 1<<LM);
1705       pitch_index = MAX_PERIOD-len2-pitch_index;
1706       st->last_pitch_index = pitch_index;
1707    } else {
1708       pitch_index = st->last_pitch_index;
1709       if (st->loss_count < 5)
1710          fade = QCONST16(.8f,15);
1711       else
1712          fade = 0;
1713    }
1714
1715    c=0; do {
1716       /* FIXME: This is more memory than necessary */
1717       celt_word32 e[2*MAX_PERIOD];
1718       celt_word16 exc[2*MAX_PERIOD];
1719       celt_word32 ac[LPC_ORDER+1];
1720       celt_word16 decay = 1;
1721       celt_word32 S1=0;
1722       celt_word16 mem[LPC_ORDER]={0};
1723
1724       offset = MAX_PERIOD-pitch_index;
1725       for (i=0;i<MAX_PERIOD;i++)
1726          exc[i] = ROUND16(out_mem[c][i], SIG_SHIFT);
1727
1728       if (st->loss_count == 0)
1729       {
1730          _celt_autocorr(exc, ac, st->mode->window, st->mode->overlap,
1731                         LPC_ORDER, MAX_PERIOD);
1732
1733          /* Noise floor -40 dB */
1734 #ifdef FIXED_POINT
1735          ac[0] += SHR32(ac[0],13);
1736 #else
1737          ac[0] *= 1.0001f;
1738 #endif
1739          /* Lag windowing */
1740          for (i=1;i<=LPC_ORDER;i++)
1741          {
1742             /*ac[i] *= exp(-.5*(2*M_PI*.002*i)*(2*M_PI*.002*i));*/
1743 #ifdef FIXED_POINT
1744             ac[i] -= MULT16_32_Q15(2*i*i, ac[i]);
1745 #else
1746             ac[i] -= ac[i]*(.008f*i)*(.008f*i);
1747 #endif
1748          }
1749
1750          _celt_lpc(lpc+c*LPC_ORDER, ac, LPC_ORDER);
1751       }
1752       for (i=0;i<LPC_ORDER;i++)
1753          mem[i] = ROUND16(out_mem[c][MAX_PERIOD-1-i], SIG_SHIFT);
1754       fir(exc, lpc+c*LPC_ORDER, exc, MAX_PERIOD, LPC_ORDER, mem);
1755       /*for (i=0;i<MAX_PERIOD;i++)printf("%d ", exc[i]); printf("\n");*/
1756       /* Check if the waveform is decaying (and if so how fast) */
1757       {
1758          celt_word32 E1=1, E2=1;
1759          int period;
1760          if (pitch_index <= MAX_PERIOD/2)
1761             period = pitch_index;
1762          else
1763             period = MAX_PERIOD/2;
1764          for (i=0;i<period;i++)
1765          {
1766             E1 += SHR32(MULT16_16(exc[MAX_PERIOD-period+i],exc[MAX_PERIOD-period+i]),8);
1767             E2 += SHR32(MULT16_16(exc[MAX_PERIOD-2*period+i],exc[MAX_PERIOD-2*period+i]),8);
1768          }
1769          if (E1 > E2)
1770             E1 = E2;
1771          decay = celt_sqrt(frac_div32(SHR(E1,1),E2));
1772       }
1773
1774       /* Copy excitation, taking decay into account */
1775       for (i=0;i<len+st->mode->overlap;i++)
1776       {
1777          celt_word16 tmp;
1778          if (offset+i >= MAX_PERIOD)
1779          {
1780             offset -= pitch_index;
1781             decay = MULT16_16_Q15(decay, decay);
1782          }
1783          e[i] = SHL32(EXTEND32(MULT16_16_Q15(decay, exc[offset+i])), SIG_SHIFT);
1784          tmp = ROUND16(out_mem[c][offset+i],SIG_SHIFT);
1785          S1 += SHR32(MULT16_16(tmp,tmp),8);
1786       }
1787       for (i=0;i<LPC_ORDER;i++)
1788          mem[i] = ROUND16(out_mem[c][MAX_PERIOD-1-i], SIG_SHIFT);
1789       for (i=0;i<len+st->mode->overlap;i++)
1790          e[i] = MULT16_32_Q15(fade, e[i]);
1791       iir(e, lpc+c*LPC_ORDER, e, len+st->mode->overlap, LPC_ORDER, mem);
1792
1793       {
1794          celt_word32 S2=0;
1795          for (i=0;i<len+overlap;i++)
1796          {
1797             celt_word16 tmp = ROUND16(e[i],SIG_SHIFT);
1798             S2 += SHR32(MULT16_16(tmp,tmp),8);
1799          }
1800          /* This checks for an "explosion" in the synthesis */
1801 #ifdef FIXED_POINT
1802          if (!(S1 > SHR32(S2,2)))
1803 #else
1804          /* Float test is written this way to catch NaNs at the same time */
1805          if (!(S1 > 0.2f*S2))
1806 #endif
1807          {
1808             for (i=0;i<len+overlap;i++)
1809                e[i] = 0;
1810          } else if (S1 < S2)
1811          {
1812             celt_word16 ratio = celt_sqrt(frac_div32(SHR32(S1,1)+1,S2+1));
1813             for (i=0;i<len+overlap;i++)
1814                e[i] = MULT16_32_Q15(ratio, e[i]);
1815          }
1816       }
1817
1818 #ifdef ENABLE_POSTFILTER
1819       /* Apply post-filter to the MDCT overlap of the previous frame */
1820       comb_filter(out_mem[c]+MAX_PERIOD, out_mem[c]+MAX_PERIOD, st->postfilter_period, st->postfilter_period, st->overlap, C,
1821                   st->postfilter_gain, st->postfilter_gain, st->postfilter_tapset, st->postfilter_tapset,
1822                   NULL, 0);
1823 #endif /* ENABLE_POSTFILTER */
1824
1825       for (i=0;i<MAX_PERIOD+st->mode->overlap-N;i++)
1826          out_mem[c][i] = out_mem[c][N+i];
1827
1828       /* Apply TDAC to the concealed audio so that it blends with the
1829          previous and next frames */
1830       for (i=0;i<overlap/2;i++)
1831       {
1832          celt_word32 tmp;
1833          tmp = MULT16_32_Q15(st->mode->window[i],           e[N+overlap-1-i]) +
1834                MULT16_32_Q15(st->mode->window[overlap-i-1], e[N+i          ]);
1835          out_mem[c][MAX_PERIOD+i] = MULT16_32_Q15(st->mode->window[overlap-i-1], tmp);
1836          out_mem[c][MAX_PERIOD+overlap-i-1] = MULT16_32_Q15(st->mode->window[i], tmp);
1837       }
1838       for (i=0;i<N;i++)
1839          out_mem[c][MAX_PERIOD-N+i] = e[i];
1840
1841 #ifdef ENABLE_POSTFILTER
1842       /* Apply pre-filter to the MDCT overlap for the next frame (post-filter will be applied then) */
1843       comb_filter(e, out_mem[c]+MAX_PERIOD, st->postfilter_period, st->postfilter_period, st->overlap, C,
1844                   -st->postfilter_gain, -st->postfilter_gain, st->postfilter_tapset, st->postfilter_tapset,
1845                   NULL, 0);
1846 #endif /* ENABLE_POSTFILTER */
1847       for (i=0;i<overlap;i++)
1848          out_mem[c][MAX_PERIOD+i] = e[i];
1849    } while (++c<C);
1850
1851    {
1852       celt_word32 *out_syn[2];
1853       out_syn[0] = out_mem[0]+MAX_PERIOD-N;
1854       if (C==2)
1855          out_syn[1] = out_mem[1]+MAX_PERIOD-N;
1856       deemphasis(out_syn, pcm, N, C, st->mode->preemph, st->preemph_memD);
1857    }
1858    
1859    st->loss_count++;
1860
1861    RESTORE_STACK;
1862 }
1863
1864 #ifdef FIXED_POINT
1865 int celt_decode_with_ec(CELTDecoder * restrict st, const unsigned char *data, int len, celt_int16 * restrict pcm, int frame_size, ec_dec *dec)
1866 {
1867 #else
1868 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)
1869 {
1870 #endif
1871    int c, i, N;
1872    int spread_decision;
1873    int bits;
1874    ec_dec _dec;
1875    ec_byte_buffer buf;
1876    VARDECL(celt_sig, freq);
1877    VARDECL(celt_norm, X);
1878    VARDECL(celt_ener, bandE);
1879    VARDECL(celt_word16, oldLogE);
1880    VARDECL(int, fine_quant);
1881    VARDECL(int, pulses);
1882    VARDECL(int, offsets);
1883    VARDECL(int, fine_priority);
1884    VARDECL(int, tf_res);
1885    celt_sig *out_mem[2];
1886    celt_sig *decode_mem[2];
1887    celt_sig *overlap_mem[2];
1888    celt_sig *out_syn[2];
1889    celt_word16 *lpc;
1890    celt_word16 *oldBandE, *oldLogE2;
1891
1892    int shortBlocks;
1893    int isTransient;
1894    int intra_ener;
1895    const int C = CHANNELS(st->channels);
1896    int LM, M;
1897    int effEnd;
1898    int codedBands;
1899    int alloc_trim;
1900    int postfilter_pitch;
1901    celt_word16 postfilter_gain;
1902    int intensity=0;
1903    int dual_stereo=0;
1904    celt_int32 total_bits;
1905    celt_int32 tell;
1906    int dynalloc_logp;
1907    int postfilter_tapset;
1908    int anti_collapse_on=0;
1909    SAVE_STACK;
1910
1911    if (pcm==NULL)
1912       return CELT_BAD_ARG;
1913
1914    for (LM=0;LM<4;LM++)
1915       if (st->mode->shortMdctSize<<LM==frame_size)
1916          break;
1917    if (LM>=MAX_CONFIG_SIZES)
1918       return CELT_BAD_ARG;
1919    M=1<<LM;
1920
1921    c=0; do {
1922       decode_mem[c] = st->_decode_mem + c*(DECODE_BUFFER_SIZE+st->overlap);
1923       out_mem[c] = decode_mem[c]+DECODE_BUFFER_SIZE-MAX_PERIOD;
1924       overlap_mem[c] = decode_mem[c]+DECODE_BUFFER_SIZE;
1925    } while (++c<C);
1926    lpc = (celt_word16*)(st->_decode_mem+(DECODE_BUFFER_SIZE+st->overlap)*C);
1927    oldBandE = lpc+C*LPC_ORDER;
1928    oldLogE2 = oldBandE + C*st->mode->nbEBands;
1929
1930    N = M*st->mode->shortMdctSize;
1931
1932    effEnd = st->end;
1933    if (effEnd > st->mode->effEBands)
1934       effEnd = st->mode->effEBands;
1935
1936    ALLOC(freq, C*N, celt_sig); /**< Interleaved signal MDCTs */
1937    ALLOC(X, C*N, celt_norm);   /**< Interleaved normalised MDCTs */
1938    ALLOC(bandE, st->mode->nbEBands*C, celt_ener);
1939    c=0; do
1940       for (i=0;i<M*st->mode->eBands[st->start];i++)
1941          X[c*N+i] = 0;
1942    while (++c<C);
1943    c=0; do   
1944       for (i=M*st->mode->eBands[effEnd];i<N;i++)
1945          X[c*N+i] = 0;
1946    while (++c<C);
1947
1948    if (data == NULL)
1949    {
1950       celt_decode_lost(st, pcm, N, LM);
1951       RESTORE_STACK;
1952       return CELT_OK;
1953    }
1954    if (len<0) {
1955      RESTORE_STACK;
1956      return CELT_BAD_ARG;
1957    }
1958    
1959    if (dec == NULL)
1960    {
1961       ec_byte_readinit(&buf,(unsigned char*)data,len);
1962       ec_dec_init(&_dec,&buf);
1963       dec = &_dec;
1964    }
1965
1966    total_bits = len*8;
1967    tell = ec_dec_tell(dec, 0);
1968
1969    postfilter_gain = 0;
1970    postfilter_pitch = 0;
1971    postfilter_tapset = 0;
1972    if (tell+17 <= total_bits)
1973    {
1974       if(ec_dec_bit_logp(dec, 1))
1975       {
1976 #ifdef ENABLE_POSTFILTER
1977          int qg, octave;
1978          octave = ec_dec_uint(dec, 6);
1979          postfilter_pitch = (16<<octave)+ec_dec_bits(dec, 4+octave)-1;
1980          qg = ec_dec_bits(dec, 2);
1981          postfilter_tapset = ec_dec_icdf(dec, tapset_icdf, 2);
1982          postfilter_gain = QCONST16(.125f,15)*(qg+2);
1983          tell = ec_dec_tell(dec, 0);
1984 #else /* ENABLE_POSTFILTER */
1985          RESTORE_STACK;
1986          return CELT_CORRUPTED_DATA;
1987 #endif /* ENABLE_POSTFILTER */
1988       }
1989       tell = ec_dec_tell(dec, 0);
1990    }
1991
1992    if (LM > 0 && tell+3 <= total_bits)
1993    {
1994       isTransient = ec_dec_bit_logp(dec, 3);
1995       tell = ec_dec_tell(dec, 0);
1996    }
1997    else
1998       isTransient = 0;
1999
2000    if (isTransient)
2001       shortBlocks = M;
2002    else
2003       shortBlocks = 0;
2004
2005    ALLOC(oldLogE, C*st->mode->nbEBands, celt_word16);
2006    for (i=0;i<C*st->mode->nbEBands;i++)
2007       oldLogE[i] = oldBandE[i];
2008
2009    /* Decode the global flags (first symbols in the stream) */
2010    intra_ener = tell+3<=total_bits ? ec_dec_bit_logp(dec, 3) : 0;
2011    /* Get band energies */
2012    unquant_coarse_energy(st->mode, st->start, st->end, bandE, oldBandE,
2013          intra_ener, dec, C, LM);
2014
2015    ALLOC(tf_res, st->mode->nbEBands, int);
2016    tf_decode(st->start, st->end, C, isTransient, tf_res, LM, dec);
2017
2018    tell = ec_dec_tell(dec, 0);
2019    spread_decision = SPREAD_NORMAL;
2020    if (tell+4 <= total_bits)
2021       spread_decision = ec_dec_icdf(dec, spread_icdf, 5);
2022
2023    ALLOC(pulses, st->mode->nbEBands, int);
2024    ALLOC(offsets, st->mode->nbEBands, int);
2025    ALLOC(fine_priority, st->mode->nbEBands, int);
2026
2027    dynalloc_logp = 6;
2028    total_bits<<=BITRES;
2029    tell = ec_dec_tell(dec, BITRES);
2030    for (i=st->start;i<st->end;i++)
2031    {
2032       int width, quanta;
2033       int dynalloc_loop_logp;
2034       int boost;
2035       width = C*(st->mode->eBands[i+1]-st->mode->eBands[i])<<LM;
2036       /* quanta is 6 bits, but no more than 1 bit/sample
2037          and no less than 1/8 bit/sample */
2038       quanta = IMIN(width<<BITRES, IMAX(6<<BITRES, width));
2039       dynalloc_loop_logp = dynalloc_logp;
2040       boost = 0;
2041       while (tell+(dynalloc_loop_logp<<BITRES) < total_bits &&
2042             boost < (64<<LM)*(C<<BITRES))
2043       {
2044          int flag;
2045          flag = ec_dec_bit_logp(dec, dynalloc_loop_logp);
2046          tell = ec_dec_tell(dec, BITRES);
2047          if (!flag)
2048             break;
2049          boost += quanta;
2050          total_bits -= quanta;
2051          dynalloc_loop_logp = 1;
2052       }
2053       offsets[i] = boost;
2054       /* Making dynalloc more likely */
2055       if (boost>0)
2056          dynalloc_logp = IMAX(2, dynalloc_logp-1);
2057    }
2058
2059    ALLOC(fine_quant, st->mode->nbEBands, int);
2060    alloc_trim = tell+(6<<BITRES) <= total_bits ?
2061          ec_dec_icdf(dec, trim_icdf, 7) : 5;
2062
2063    bits = (len*8<<BITRES) - ec_dec_tell(dec, BITRES) - 1 - (isTransient&&LM>=2 ? (1<<BITRES) : 0);
2064    codedBands = compute_allocation(st->mode, st->start, st->end, offsets,
2065          alloc_trim, &intensity, &dual_stereo, bits, pulses, fine_quant,
2066          fine_priority, C, LM, dec, 0, 0);
2067    
2068    unquant_fine_energy(st->mode, st->start, st->end, bandE, oldBandE, fine_quant, dec, C);
2069
2070    /* Decode fixed codebook */
2071    quant_all_bands(0, st->mode, st->start, st->end, X, C==2 ? X+N : NULL,
2072          NULL, pulses, shortBlocks, spread_decision, dual_stereo, intensity, tf_res, 1,
2073          len*8, dec, LM, codedBands);
2074
2075    if (isTransient && LM>=2)
2076    {
2077       anti_collapse_on = ec_dec_bits(dec, 1);
2078    }
2079
2080    unquant_energy_finalise(st->mode, st->start, st->end, bandE, oldBandE,
2081          fine_quant, fine_priority, len*8-ec_dec_tell(dec, 0), dec, C);
2082
2083    if (anti_collapse_on)
2084       anti_collapse(st->mode, X, LM, C, N,
2085             st->start, st->end, oldBandE, oldLogE, oldLogE2, pulses, dec->rng);
2086
2087    log2Amp(st->mode, st->start, st->end, bandE, oldBandE, C);
2088
2089    /* Synthesis */
2090    denormalise_bands(st->mode, X, freq, bandE, effEnd, C, M);
2091
2092    CELT_MOVE(decode_mem[0], decode_mem[0]+N, DECODE_BUFFER_SIZE-N);
2093    if (C==2)
2094       CELT_MOVE(decode_mem[1], decode_mem[1]+N, DECODE_BUFFER_SIZE-N);
2095
2096    c=0; do
2097       for (i=0;i<M*st->mode->eBands[st->start];i++)
2098          freq[c*N+i] = 0;
2099    while (++c<C);
2100    c=0; do
2101       for (i=M*st->mode->eBands[effEnd];i<N;i++)
2102          freq[c*N+i] = 0;
2103    while (++c<C);
2104
2105    out_syn[0] = out_mem[0]+MAX_PERIOD-N;
2106    if (C==2)
2107       out_syn[1] = out_mem[1]+MAX_PERIOD-N;
2108
2109    /* Compute inverse MDCTs */
2110    compute_inv_mdcts(st->mode, shortBlocks, freq, out_syn, overlap_mem, C, LM);
2111
2112 #ifdef ENABLE_POSTFILTER
2113    c=0; do {
2114       st->postfilter_period=IMAX(st->postfilter_period, COMBFILTER_MINPERIOD);
2115       st->postfilter_period_old=IMAX(st->postfilter_period_old, COMBFILTER_MINPERIOD);
2116       if (LM!=0)
2117       {
2118          comb_filter(out_syn[c], out_syn[c], st->postfilter_period, st->postfilter_period, st->overlap, C,
2119                st->postfilter_gain, st->postfilter_gain, st->postfilter_tapset, st->postfilter_tapset,
2120                NULL, 0);
2121          comb_filter(out_syn[c]+st->overlap, out_syn[c]+st->overlap, st->postfilter_period, postfilter_pitch, N-st->overlap, C,
2122                st->postfilter_gain, postfilter_gain, st->postfilter_tapset, postfilter_tapset,
2123                st->mode->window, st->mode->overlap);
2124       } else {
2125          comb_filter(out_syn[c], out_syn[c], st->postfilter_period_old, st->postfilter_period, N-st->overlap, C,
2126                st->postfilter_gain_old, st->postfilter_gain, st->postfilter_tapset_old, st->postfilter_tapset,
2127                st->mode->window, st->mode->overlap);
2128       }
2129    } while (++c<C);
2130    st->postfilter_period_old = st->postfilter_period;
2131    st->postfilter_gain_old = st->postfilter_gain;
2132    st->postfilter_tapset_old = st->postfilter_tapset;
2133    st->postfilter_period = postfilter_pitch;
2134    st->postfilter_gain = postfilter_gain;
2135    st->postfilter_tapset = postfilter_tapset;
2136 #endif /* ENABLE_POSTFILTER */
2137
2138    /* In case start or end were to change */
2139    for (i=0;i<st->start;i++)
2140       oldBandE[i]=0;
2141    for (i=st->end;i<st->mode->nbEBands;i++)
2142       oldBandE[i]=0;
2143    for (i=0;i<C*st->mode->nbEBands;i++)
2144       oldLogE2[i] = oldLogE[i];
2145
2146    deemphasis(out_syn, pcm, N, C, st->mode->preemph, st->preemph_memD);
2147    st->loss_count = 0;
2148    RESTORE_STACK;
2149    if (ec_dec_tell(dec,0) > 8*len || ec_dec_get_error(dec))
2150       return CELT_CORRUPTED_DATA;
2151    else
2152       return CELT_OK;
2153 }
2154
2155 #ifdef FIXED_POINT
2156 #ifndef DISABLE_FLOAT_API
2157 int celt_decode_with_ec_float(CELTDecoder * restrict st, const unsigned char *data, int len, float * restrict pcm, int frame_size, ec_dec *dec)
2158 {
2159    int j, ret, C, N, LM, M;
2160    VARDECL(celt_int16, out);
2161    SAVE_STACK;
2162
2163    if (pcm==NULL)
2164       return CELT_BAD_ARG;
2165
2166    for (LM=0;LM<4;LM++)
2167       if (st->mode->shortMdctSize<<LM==frame_size)
2168          break;
2169    if (LM>=MAX_CONFIG_SIZES)
2170       return CELT_BAD_ARG;
2171    M=1<<LM;
2172
2173    C = CHANNELS(st->channels);
2174    N = M*st->mode->shortMdctSize;
2175    
2176    ALLOC(out, C*N, celt_int16);
2177    ret=celt_decode_with_ec(st, data, len, out, frame_size, dec);
2178    if (ret==0)
2179       for (j=0;j<C*N;j++)
2180          pcm[j]=out[j]*(1.f/32768.f);
2181      
2182    RESTORE_STACK;
2183    return ret;
2184 }
2185 #endif /*DISABLE_FLOAT_API*/
2186 #else
2187 int celt_decode_with_ec(CELTDecoder * restrict st, const unsigned char *data, int len, celt_int16 * restrict pcm, int frame_size, ec_dec *dec)
2188 {
2189    int j, ret, C, N, LM, M;
2190    VARDECL(celt_sig, out);
2191    SAVE_STACK;
2192
2193    if (pcm==NULL)
2194       return CELT_BAD_ARG;
2195
2196    for (LM=0;LM<4;LM++)
2197       if (st->mode->shortMdctSize<<LM==frame_size)
2198          break;
2199    if (LM>=MAX_CONFIG_SIZES)
2200       return CELT_BAD_ARG;
2201    M=1<<LM;
2202
2203    C = CHANNELS(st->channels);
2204    N = M*st->mode->shortMdctSize;
2205    ALLOC(out, C*N, celt_sig);
2206
2207    ret=celt_decode_with_ec_float(st, data, len, out, frame_size, dec);
2208
2209    if (ret==0)
2210       for (j=0;j<C*N;j++)
2211          pcm[j] = FLOAT2INT16 (out[j]);
2212    
2213    RESTORE_STACK;
2214    return ret;
2215 }
2216 #endif
2217
2218 int celt_decode(CELTDecoder * restrict st, const unsigned char *data, int len, celt_int16 * restrict pcm, int frame_size)
2219 {
2220    return celt_decode_with_ec(st, data, len, pcm, frame_size, NULL);
2221 }
2222
2223 #ifndef DISABLE_FLOAT_API
2224 int celt_decode_float(CELTDecoder * restrict st, const unsigned char *data, int len, float * restrict pcm, int frame_size)
2225 {
2226    return celt_decode_with_ec_float(st, data, len, pcm, frame_size, NULL);
2227 }
2228 #endif /* DISABLE_FLOAT_API */
2229
2230 int celt_decoder_ctl(CELTDecoder * restrict st, int request, ...)
2231 {
2232    va_list ap;
2233
2234    va_start(ap, request);
2235    switch (request)
2236    {
2237       case CELT_GET_MODE_REQUEST:
2238       {
2239          const CELTMode ** value = va_arg(ap, const CELTMode**);
2240          if (value==0)
2241             goto bad_arg;
2242          *value=st->mode;
2243       }
2244       break;
2245       case CELT_SET_START_BAND_REQUEST:
2246       {
2247          celt_int32 value = va_arg(ap, celt_int32);
2248          if (value<0 || value>=st->mode->nbEBands)
2249             goto bad_arg;
2250          st->start = value;
2251       }
2252       break;
2253       case CELT_SET_END_BAND_REQUEST:
2254       {
2255          celt_int32 value = va_arg(ap, celt_int32);
2256          if (value<0 || value>=st->mode->nbEBands)
2257             goto bad_arg;
2258          st->end = value;
2259       }
2260       break;
2261       case CELT_RESET_STATE:
2262       {
2263          CELT_MEMSET((char*)&st->DECODER_RESET_START, 0,
2264                celt_decoder_get_size(st->mode, st->channels)-
2265                ((char*)&st->DECODER_RESET_START - (char*)st));
2266       }
2267       break;
2268       default:
2269          goto bad_request;
2270    }
2271    va_end(ap);
2272    return CELT_OK;
2273 bad_arg:
2274    va_end(ap);
2275    return CELT_BAD_ARG;
2276 bad_request:
2277       va_end(ap);
2278   return CELT_UNIMPLEMENTED;
2279 }
2280
2281 const char *celt_strerror(int error)
2282 {
2283    static const char *error_strings[8] = {
2284       "success",
2285       "invalid argument",
2286       "invalid mode",
2287       "internal error",
2288       "corrupted stream",
2289       "request not implemented",
2290       "invalid state",
2291       "memory allocation failed"
2292    };
2293    if (error > 0 || error < -7)
2294       return "unknown error";
2295    else 
2296       return error_strings[-error];
2297 }
2298