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