Moving energy floor to coarse quantization
[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 16
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          + 2*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(celt_word16, oldLogE);
791    VARDECL(int, fine_quant);
792    VARDECL(celt_word16, error);
793    VARDECL(int, pulses);
794    VARDECL(int, offsets);
795    VARDECL(int, fine_priority);
796    VARDECL(int, tf_res);
797    VARDECL(unsigned char, collapse_masks);
798    celt_sig *_overlap_mem;
799    celt_sig *prefilter_mem;
800    celt_word16 *oldBandE, *oldLogE2;
801    int shortBlocks=0;
802    int isTransient=0;
803    int resynth;
804    const int C = CHANNELS(st->channels);
805    int LM, M;
806    int tf_select;
807    int nbFilledBytes, nbAvailableBytes;
808    int effEnd;
809    int codedBands;
810    int tf_sum;
811    int alloc_trim;
812    int pitch_index=COMBFILTER_MINPERIOD;
813    celt_word16 gain1 = 0;
814    int intensity=0;
815    int dual_stereo=0;
816    int effectiveBytes;
817    celt_word16 pf_threshold;
818    int dynalloc_logp;
819    celt_int32 vbr_rate;
820    celt_int32 total_bits;
821    celt_int32 total_boost;
822    celt_int32 tell;
823    int prefilter_tapset=0;
824    int pf_on;
825    int anti_collapse_rsv;
826    int anti_collapse_on=0;
827    SAVE_STACK;
828
829    if (nbCompressedBytes<0 || pcm==NULL)
830      return CELT_BAD_ARG;
831
832    for (LM=0;LM<4;LM++)
833       if (st->mode->shortMdctSize<<LM==frame_size)
834          break;
835    if (LM>=MAX_CONFIG_SIZES)
836       return CELT_BAD_ARG;
837    M=1<<LM;
838
839    prefilter_mem = st->in_mem+C*(st->overlap);
840    _overlap_mem = prefilter_mem+C*COMBFILTER_MAXPERIOD;
841    /*_overlap_mem = st->in_mem+C*(st->overlap);*/
842    oldBandE = (celt_word16*)(st->in_mem+C*(2*st->overlap+COMBFILTER_MAXPERIOD));
843    oldLogE2 = oldBandE + 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(oldLogE, C*st->mode->nbEBands, celt_word16);
1071    for (i=0;i<C*st->mode->nbEBands;i++)
1072       oldLogE[i] = oldBandE[i];
1073    ALLOC(error, C*st->mode->nbEBands, celt_word16);
1074    quant_coarse_energy(st->mode, st->start, st->end, effEnd, bandLogE,
1075          oldBandE, total_bits, error, enc,
1076          C, LM, nbAvailableBytes, st->force_intra,
1077          &st->delayedIntra, st->complexity >= 4);
1078
1079    tf_encode(st->start, st->end, isTransient, tf_res, LM, tf_select, enc);
1080
1081    st->spread_decision = SPREAD_NORMAL;
1082    if (ec_enc_tell(enc, 0)+4<=total_bits)
1083    {
1084       if (shortBlocks || st->complexity < 3 || nbAvailableBytes < 10*C)
1085       {
1086          if (st->complexity == 0)
1087             st->spread_decision = SPREAD_NONE;
1088       } else {
1089          st->spread_decision = spreading_decision(st->mode, X,
1090                &st->tonal_average, st->spread_decision, &st->hf_average,
1091                &st->tapset_decision, pf_on&&!shortBlocks, effEnd, C, M);
1092       }
1093       ec_enc_icdf(enc, st->spread_decision, spread_icdf, 5);
1094    }
1095
1096    ALLOC(offsets, st->mode->nbEBands, int);
1097
1098    for (i=0;i<st->mode->nbEBands;i++)
1099       offsets[i] = 0;
1100    /* Dynamic allocation code */
1101    /* Make sure that dynamic allocation can't make us bust the budget */
1102    if (effectiveBytes > 50 && LM>=1)
1103    {
1104       int t1, t2;
1105       if (LM <= 1)
1106       {
1107          t1 = 3;
1108          t2 = 5;
1109       } else {
1110          t1 = 2;
1111          t2 = 4;
1112       }
1113       for (i=1;i<st->mode->nbEBands-1;i++)
1114       {
1115          celt_word32 d2;
1116          d2 = 2*bandLogE[i]-bandLogE[i-1]-bandLogE[i+1];
1117          if (C==2)
1118             d2 = HALF32(d2 + 2*bandLogE[i+st->mode->nbEBands]-
1119                   bandLogE[i-1+st->mode->nbEBands]-bandLogE[i+1+st->mode->nbEBands]);
1120          if (d2 > SHL16(t1,DB_SHIFT))
1121             offsets[i] += 1;
1122          if (d2 > SHL16(t2,DB_SHIFT))
1123             offsets[i] += 1;
1124       }
1125    }
1126    dynalloc_logp = 6;
1127    total_bits<<=BITRES;
1128    total_boost = 0;
1129    tell = ec_enc_tell(enc, BITRES);
1130    for (i=st->start;i<st->end;i++)
1131    {
1132       int width, quanta;
1133       int dynalloc_loop_logp;
1134       int boost;
1135       int j;
1136       width = C*(st->mode->eBands[i+1]-st->mode->eBands[i])<<LM;
1137       /* quanta is 6 bits, but no more than 1 bit/sample
1138          and no less than 1/8 bit/sample */
1139       quanta = IMIN(width<<BITRES, IMAX(6<<BITRES, width));
1140       dynalloc_loop_logp = dynalloc_logp;
1141       boost = 0;
1142       for (j = 0; tell+(dynalloc_loop_logp<<BITRES) < total_bits-total_boost
1143             && boost < (64<<LM)*(C<<BITRES); j++)
1144       {
1145          int flag;
1146          flag = j<offsets[i];
1147          ec_enc_bit_logp(enc, flag, dynalloc_loop_logp);
1148          tell = ec_enc_tell(enc, BITRES);
1149          if (!flag)
1150             break;
1151          boost += quanta;
1152          total_boost += quanta;
1153          dynalloc_loop_logp = 1;
1154       }
1155       /* Making dynalloc more likely */
1156       if (j)
1157          dynalloc_logp = IMAX(2, dynalloc_logp-1);
1158       offsets[i] = boost;
1159    }
1160    alloc_trim = 5;
1161    if (tell+(6<<BITRES) <= total_bits - total_boost)
1162    {
1163       alloc_trim = alloc_trim_analysis(st->mode, X, bandLogE,
1164             st->mode->nbEBands, LM, C, N);
1165       ec_enc_icdf(enc, alloc_trim, trim_icdf, 7);
1166       tell = ec_enc_tell(enc, BITRES);
1167    }
1168
1169    /* Variable bitrate */
1170    if (vbr_rate>0)
1171    {
1172      celt_word16 alpha;
1173      celt_int32 delta;
1174      /* The target rate in 8th bits per frame */
1175      celt_int32 target;
1176      celt_int32 min_allowed;
1177
1178      target = vbr_rate + st->vbr_offset - ((40*C+20)<<BITRES);
1179
1180      /* Shortblocks get a large boost in bitrate, but since they
1181         are uncommon long blocks are not greatly affected */
1182      if (shortBlocks || tf_sum < -2*(st->end-st->start))
1183         target = 7*target/4;
1184      else if (tf_sum < -(st->end-st->start))
1185         target = 3*target/2;
1186      else if (M > 1)
1187         target-=(target+14)/28;
1188
1189      /* The current offset is removed from the target and the space used
1190         so far is added*/
1191      target=target+tell;
1192      /* By how much did we "miss" the target on that frame */
1193      delta = target - vbr_rate;
1194
1195      /* In VBR mode the frame size must not be reduced so much that it would
1196          result in the encoder running out of bits.
1197         The margin of 2 bytes ensures that none of the bust-prevention logic
1198          in the decoder will have triggered so far. */
1199      min_allowed = (tell+total_boost+(1<<BITRES+3)-1>>(BITRES+3)) + 2 - nbFilledBytes;
1200
1201      nbAvailableBytes = target+(1<<(BITRES+2))>>(BITRES+3);
1202      nbAvailableBytes = IMAX(min_allowed,nbAvailableBytes);
1203      nbAvailableBytes = IMIN(nbCompressedBytes,nbAvailableBytes+nbFilledBytes) - nbFilledBytes;
1204
1205      target=nbAvailableBytes<<(BITRES+3);
1206
1207      if (st->vbr_count < 970)
1208      {
1209         st->vbr_count++;
1210         alpha = celt_rcp(SHL32(EXTEND32(st->vbr_count+20),16));
1211      } else
1212         alpha = QCONST16(.001f,15);
1213      /* How many bits have we used in excess of what we're allowed */
1214      if (st->constrained_vbr)
1215         st->vbr_reservoir += target - vbr_rate;
1216      /*printf ("%d\n", st->vbr_reservoir);*/
1217
1218      /* Compute the offset we need to apply in order to reach the target */
1219      st->vbr_drift += (celt_int32)MULT16_32_Q15(alpha,delta-st->vbr_offset-st->vbr_drift);
1220      st->vbr_offset = -st->vbr_drift;
1221      /*printf ("%d\n", st->vbr_drift);*/
1222
1223      if (st->constrained_vbr && st->vbr_reservoir < 0)
1224      {
1225         /* We're under the min value -- increase rate */
1226         int adjust = (-st->vbr_reservoir)/(8<<BITRES);
1227         nbAvailableBytes += adjust;
1228         st->vbr_reservoir = 0;
1229         /*printf ("+%d\n", adjust);*/
1230      }
1231      nbCompressedBytes = IMIN(nbCompressedBytes,nbAvailableBytes+nbFilledBytes);
1232
1233      /* This moves the raw bits to take into account the new compressed size */
1234      ec_byte_shrink(&buf, nbCompressedBytes);
1235    }
1236
1237    if (C==2)
1238    {
1239       int effectiveRate;
1240
1241       /* Always use MS for 2.5 ms frames until we can do a better analysis */
1242       if (LM!=0)
1243          dual_stereo = stereo_analysis(st->mode, X, st->mode->nbEBands, LM, C, N);
1244
1245       /* Account for coarse energy */
1246       effectiveRate = (8*effectiveBytes - 80)>>LM;
1247
1248       /* effectiveRate in kb/s */
1249       effectiveRate = 2*effectiveRate/5;
1250       if (effectiveRate<35)
1251          intensity = 8;
1252       else if (effectiveRate<50)
1253          intensity = 12;
1254       else if (effectiveRate<68)
1255          intensity = 16;
1256       else if (effectiveRate<84)
1257          intensity = 18;
1258       else if (effectiveRate<102)
1259          intensity = 19;
1260       else if (effectiveRate<130)
1261          intensity = 20;
1262       else
1263          intensity = 100;
1264       intensity = IMIN(st->end,IMAX(st->start, intensity));
1265    }
1266
1267    /* Bit allocation */
1268    ALLOC(fine_quant, st->mode->nbEBands, int);
1269    ALLOC(pulses, st->mode->nbEBands, int);
1270    ALLOC(fine_priority, st->mode->nbEBands, int);
1271
1272    /* bits =   packet size        -       where we are         - safety*/
1273    bits = (nbCompressedBytes*8<<BITRES) - ec_enc_tell(enc, BITRES) - 1;
1274    anti_collapse_rsv = isTransient&&LM>=2&&bits>=(LM+2<<BITRES) ? (1<<BITRES) : 0;
1275    bits -= anti_collapse_rsv;
1276    codedBands = compute_allocation(st->mode, st->start, st->end, offsets,
1277          alloc_trim, &intensity, &dual_stereo, bits, pulses, fine_quant,
1278          fine_priority, C, LM, enc, 1, st->lastCodedBands);
1279    st->lastCodedBands = codedBands;
1280
1281    quant_fine_energy(st->mode, st->start, st->end, bandE, oldBandE, error, fine_quant, enc, C);
1282
1283 #ifdef MEASURE_NORM_MSE
1284    float X0[3000];
1285    float bandE0[60];
1286    c=0; do 
1287       for (i=0;i<N;i++)
1288          X0[i+c*N] = X[i+c*N];
1289    while (++c<C);
1290    for (i=0;i<C*st->mode->nbEBands;i++)
1291       bandE0[i] = bandE[i];
1292 #endif
1293
1294    /* Residual quantisation */
1295    ALLOC(collapse_masks, st->mode->nbEBands, unsigned char);
1296    quant_all_bands(1, st->mode, st->start, st->end, X, C==2 ? X+N : NULL, collapse_masks,
1297          bandE, pulses, shortBlocks, st->spread_decision, dual_stereo, intensity, tf_res, resynth,
1298          nbCompressedBytes*8, enc, LM, codedBands, &st->rng);
1299
1300    if (anti_collapse_rsv > 0)
1301    {
1302       anti_collapse_on = st->consec_transient<2;
1303       ec_enc_bits(enc, anti_collapse_on, 1);
1304    }
1305    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);
1306
1307 #ifdef RESYNTH
1308    /* Re-synthesis of the coded audio if required */
1309    if (resynth)
1310    {
1311       celt_sig *out_mem[2];
1312       celt_sig *overlap_mem[2];
1313
1314       log2Amp(st->mode, st->start, st->end, bandE, oldBandE, C);
1315
1316 #ifdef MEASURE_NORM_MSE
1317       measure_norm_mse(st->mode, X, X0, bandE, bandE0, M, N, C);
1318 #endif
1319       if (anti_collapse_on)
1320       {
1321          anti_collapse(st->mode, X, collapse_masks, LM, C, N,
1322                st->start, st->end, oldBandE, oldLogE, oldLogE2, pulses, st->rng);
1323       }
1324
1325       /* Synthesis */
1326       denormalise_bands(st->mode, X, freq, bandE, effEnd, C, M);
1327
1328       CELT_MOVE(st->syn_mem[0], st->syn_mem[0]+N, MAX_PERIOD);
1329       if (C==2)
1330          CELT_MOVE(st->syn_mem[1], st->syn_mem[1]+N, MAX_PERIOD);
1331
1332       c=0; do
1333          for (i=0;i<M*st->mode->eBands[st->start];i++)
1334             freq[c*N+i] = 0;
1335       while (++c<C);
1336       c=0; do
1337          for (i=M*st->mode->eBands[st->end];i<N;i++)
1338             freq[c*N+i] = 0;
1339       while (++c<C);
1340
1341       out_mem[0] = st->syn_mem[0]+MAX_PERIOD;
1342       if (C==2)
1343          out_mem[1] = st->syn_mem[1]+MAX_PERIOD;
1344
1345       c=0; do
1346          overlap_mem[c] = _overlap_mem + c*st->overlap;
1347       while (++c<C);
1348
1349       compute_inv_mdcts(st->mode, shortBlocks, freq, out_mem, overlap_mem, C, LM);
1350
1351 #ifdef ENABLE_POSTFILTER
1352       c=0; do {
1353          st->prefilter_period=IMAX(st->prefilter_period, COMBFILTER_MINPERIOD);
1354          st->prefilter_period_old=IMAX(st->prefilter_period_old, COMBFILTER_MINPERIOD);
1355          if (LM!=0)
1356          {
1357             comb_filter(out_mem[c], out_mem[c], st->prefilter_period, st->prefilter_period, st->overlap, C,
1358                   st->prefilter_gain, st->prefilter_gain, st->prefilter_tapset, st->prefilter_tapset,
1359                   NULL, 0);
1360             comb_filter(out_mem[c]+st->overlap, out_mem[c]+st->overlap, st->prefilter_period, pitch_index, N-st->overlap, C,
1361                   st->prefilter_gain, gain1, st->prefilter_tapset, prefilter_tapset,
1362                   st->mode->window, st->mode->overlap);
1363          } else {
1364             comb_filter(out_mem[c], out_mem[c], st->prefilter_period_old, st->prefilter_period, N, C,
1365                   st->prefilter_gain_old, st->prefilter_gain, st->prefilter_tapset_old, st->prefilter_tapset,
1366                   st->mode->window, st->mode->overlap);
1367          }
1368       } while (++c<C);
1369 #endif /* ENABLE_POSTFILTER */
1370
1371       deemphasis(out_mem, (celt_word16*)pcm, N, C, st->mode->preemph, st->preemph_memD);
1372       st->prefilter_period_old = st->prefilter_period;
1373       st->prefilter_gain_old = st->prefilter_gain;
1374       st->prefilter_tapset_old = st->prefilter_tapset;
1375    }
1376 #endif
1377
1378    st->prefilter_period = pitch_index;
1379    st->prefilter_gain = gain1;
1380    st->prefilter_tapset = prefilter_tapset;
1381
1382    /* In case start or end were to change */
1383    c=0; do
1384    {
1385       for (i=0;i<st->start;i++)
1386          oldBandE[c*st->mode->nbEBands+i]=0;
1387       for (i=st->end;i<st->mode->nbEBands;i++)
1388          oldBandE[c*st->mode->nbEBands+i]=0;
1389    } while (++c<C);
1390    for (i=0;i<C*st->mode->nbEBands;i++)
1391       oldLogE2[i] = oldLogE[i];
1392    if (isTransient)
1393       st->consec_transient++;
1394    else
1395       st->consec_transient=0;
1396    st->rng = enc->rng;
1397
1398    /* If there's any room left (can only happen for very high rates),
1399       fill it with zeros */
1400    while (ec_enc_tell(enc,0) + 8 <= nbCompressedBytes*8)
1401       ec_enc_bits(enc, 0, 8);
1402    ec_enc_done(enc);
1403    
1404    RESTORE_STACK;
1405    if (ec_enc_get_error(enc))
1406       return CELT_CORRUPTED_DATA;
1407    else
1408       return nbCompressedBytes;
1409 }
1410
1411 #ifdef FIXED_POINT
1412 #ifndef DISABLE_FLOAT_API
1413 int celt_encode_with_ec_float(CELTEncoder * restrict st, const float * pcm, int frame_size, unsigned char *compressed, int nbCompressedBytes, ec_enc *enc)
1414 {
1415    int j, ret, C, N, LM, M;
1416    VARDECL(celt_int16, in);
1417    SAVE_STACK;
1418
1419    if (pcm==NULL)
1420       return CELT_BAD_ARG;
1421
1422    for (LM=0;LM<4;LM++)
1423       if (st->mode->shortMdctSize<<LM==frame_size)
1424          break;
1425    if (LM>=MAX_CONFIG_SIZES)
1426       return CELT_BAD_ARG;
1427    M=1<<LM;
1428
1429    C = CHANNELS(st->channels);
1430    N = M*st->mode->shortMdctSize;
1431    ALLOC(in, C*N, celt_int16);
1432
1433    for (j=0;j<C*N;j++)
1434      in[j] = FLOAT2INT16(pcm[j]);
1435
1436    ret=celt_encode_with_ec(st,in,frame_size,compressed,nbCompressedBytes, enc);
1437 #ifdef RESYNTH
1438    for (j=0;j<C*N;j++)
1439       ((float*)pcm)[j]=in[j]*(1.f/32768.f);
1440 #endif
1441    RESTORE_STACK;
1442    return ret;
1443
1444 }
1445 #endif /*DISABLE_FLOAT_API*/
1446 #else
1447 int celt_encode_with_ec(CELTEncoder * restrict st, const celt_int16 * pcm, int frame_size, unsigned char *compressed, int nbCompressedBytes, ec_enc *enc)
1448 {
1449    int j, ret, C, N, LM, M;
1450    VARDECL(celt_sig, in);
1451    SAVE_STACK;
1452
1453    if (pcm==NULL)
1454       return CELT_BAD_ARG;
1455
1456    for (LM=0;LM<4;LM++)
1457       if (st->mode->shortMdctSize<<LM==frame_size)
1458          break;
1459    if (LM>=MAX_CONFIG_SIZES)
1460       return CELT_BAD_ARG;
1461    M=1<<LM;
1462
1463    C=CHANNELS(st->channels);
1464    N=M*st->mode->shortMdctSize;
1465    ALLOC(in, C*N, celt_sig);
1466    for (j=0;j<C*N;j++) {
1467      in[j] = SCALEOUT(pcm[j]);
1468    }
1469
1470    ret = celt_encode_with_ec_float(st,in,frame_size,compressed,nbCompressedBytes, enc);
1471 #ifdef RESYNTH
1472    for (j=0;j<C*N;j++)
1473       ((celt_int16*)pcm)[j] = FLOAT2INT16(in[j]);
1474 #endif
1475    RESTORE_STACK;
1476    return ret;
1477 }
1478 #endif
1479
1480 int celt_encode(CELTEncoder * restrict st, const celt_int16 * pcm, int frame_size, unsigned char *compressed, int nbCompressedBytes)
1481 {
1482    return celt_encode_with_ec(st, pcm, frame_size, compressed, nbCompressedBytes, NULL);
1483 }
1484
1485 #ifndef DISABLE_FLOAT_API
1486 int celt_encode_float(CELTEncoder * restrict st, const float * pcm, int frame_size, unsigned char *compressed, int nbCompressedBytes)
1487 {
1488    return celt_encode_with_ec_float(st, pcm, frame_size, compressed, nbCompressedBytes, NULL);
1489 }
1490 #endif /* DISABLE_FLOAT_API */
1491
1492 int celt_encoder_ctl(CELTEncoder * restrict st, int request, ...)
1493 {
1494    va_list ap;
1495    
1496    va_start(ap, request);
1497    switch (request)
1498    {
1499       case CELT_GET_MODE_REQUEST:
1500       {
1501          const CELTMode ** value = va_arg(ap, const CELTMode**);
1502          if (value==0)
1503             goto bad_arg;
1504          *value=st->mode;
1505       }
1506       break;
1507       case CELT_SET_COMPLEXITY_REQUEST:
1508       {
1509          int value = va_arg(ap, celt_int32);
1510          if (value<0 || value>10)
1511             goto bad_arg;
1512          st->complexity = value;
1513       }
1514       break;
1515       case CELT_SET_START_BAND_REQUEST:
1516       {
1517          celt_int32 value = va_arg(ap, celt_int32);
1518          if (value<0 || value>=st->mode->nbEBands)
1519             goto bad_arg;
1520          st->start = value;
1521       }
1522       break;
1523       case CELT_SET_END_BAND_REQUEST:
1524       {
1525          celt_int32 value = va_arg(ap, celt_int32);
1526          if (value<1 || value>st->mode->nbEBands)
1527             goto bad_arg;
1528          st->end = value;
1529       }
1530       break;
1531       case CELT_SET_PREDICTION_REQUEST:
1532       {
1533          int value = va_arg(ap, celt_int32);
1534          if (value<0 || value>2)
1535             goto bad_arg;
1536          if (value==0)
1537          {
1538             st->force_intra   = 1;
1539          } else if (value==1) {
1540             st->force_intra   = 0;
1541          } else {
1542             st->force_intra   = 0;
1543          }   
1544       }
1545       break;
1546       case CELT_SET_VBR_CONSTRAINT_REQUEST:
1547       {
1548          celt_int32 value = va_arg(ap, celt_int32);
1549          st->constrained_vbr = value;
1550       }
1551       break;
1552       case CELT_SET_VBR_RATE_REQUEST:
1553       {
1554          celt_int32 value = va_arg(ap, celt_int32);
1555          int frame_rate;
1556          int N = st->mode->shortMdctSize;
1557          if (value<0)
1558             goto bad_arg;
1559          if (value>3072000)
1560             value = 3072000;
1561          frame_rate = ((st->mode->Fs<<3)+(N>>1))/N;
1562          st->vbr_rate_norm = ((value<<(BITRES+3))+(frame_rate>>1))/frame_rate;
1563       }
1564       break;
1565       case CELT_RESET_STATE:
1566       {
1567          CELT_MEMSET((char*)&st->ENCODER_RESET_START, 0,
1568                celt_encoder_get_size(st->mode, st->channels)-
1569                ((char*)&st->ENCODER_RESET_START - (char*)st));
1570          st->vbr_offset = 0;
1571          st->delayedIntra = 1;
1572          st->spread_decision = SPREAD_NORMAL;
1573          st->tonal_average = QCONST16(1.f,8);
1574       }
1575       break;
1576       default:
1577          goto bad_request;
1578    }
1579    va_end(ap);
1580    return CELT_OK;
1581 bad_arg:
1582    va_end(ap);
1583    return CELT_BAD_ARG;
1584 bad_request:
1585    va_end(ap);
1586    return CELT_UNIMPLEMENTED;
1587 }
1588
1589 /**********************************************************************/
1590 /*                                                                    */
1591 /*                             DECODER                                */
1592 /*                                                                    */
1593 /**********************************************************************/
1594 #define DECODE_BUFFER_SIZE 2048
1595
1596 /** Decoder state 
1597  @brief Decoder state
1598  */
1599 struct CELTDecoder {
1600    const CELTMode *mode;
1601    int overlap;
1602    int channels;
1603
1604    int start, end;
1605
1606    /* Everything beyond this point gets cleared on a reset */
1607 #define DECODER_RESET_START rng
1608
1609    ec_uint32 rng;
1610    int last_pitch_index;
1611    int loss_count;
1612    int postfilter_period;
1613    int postfilter_period_old;
1614    celt_word16 postfilter_gain;
1615    celt_word16 postfilter_gain_old;
1616    int postfilter_tapset;
1617    int postfilter_tapset_old;
1618
1619    celt_sig preemph_memD[2];
1620    
1621    celt_sig _decode_mem[1]; /* Size = channels*(DECODE_BUFFER_SIZE+mode->overlap) */
1622    /* celt_word16 lpc[],  Size = channels*LPC_ORDER */
1623    /* celt_word16 oldEBands[], Size = channels*mode->nbEBands */
1624    /* celt_word16 oldLogE2[], Size = channels*mode->nbEBands */
1625    /* celt_word16 backgroundLogE[], Size = channels*mode->nbEBands */
1626 };
1627
1628 int celt_decoder_get_size(const CELTMode *mode, int channels)
1629 {
1630    int size = sizeof(struct CELTDecoder)
1631             + (channels*(DECODE_BUFFER_SIZE+mode->overlap)-1)*sizeof(celt_sig)
1632             + channels*LPC_ORDER*sizeof(celt_word16)
1633             + 3*channels*mode->nbEBands*sizeof(celt_word16);
1634    return size;
1635 }
1636
1637 CELTDecoder *celt_decoder_create(const CELTMode *mode, int channels, int *error)
1638 {
1639    return celt_decoder_init(
1640          (CELTDecoder *)celt_alloc(celt_decoder_get_size(mode, channels)),
1641          mode, channels, error);
1642 }
1643
1644 CELTDecoder *celt_decoder_init(CELTDecoder *st, const CELTMode *mode, int channels, int *error)
1645 {
1646    if (channels < 0 || channels > 2)
1647    {
1648       if (error)
1649          *error = CELT_BAD_ARG;
1650       return NULL;
1651    }
1652
1653    if (st==NULL)
1654    {
1655       if (error)
1656          *error = CELT_ALLOC_FAIL;
1657       return NULL;
1658    }
1659
1660    CELT_MEMSET((char*)st, 0, celt_decoder_get_size(mode, channels));
1661
1662    st->mode = mode;
1663    st->overlap = mode->overlap;
1664    st->channels = channels;
1665
1666    st->start = 0;
1667    st->end = st->mode->effEBands;
1668
1669    st->loss_count = 0;
1670
1671    if (error)
1672       *error = CELT_OK;
1673    return st;
1674 }
1675
1676 void celt_decoder_destroy(CELTDecoder *st)
1677 {
1678    celt_free(st);
1679 }
1680
1681 static void celt_decode_lost(CELTDecoder * restrict st, celt_word16 * restrict pcm, int N, int LM)
1682 {
1683    int c;
1684    int pitch_index;
1685    int overlap = st->mode->overlap;
1686    celt_word16 fade = Q15ONE;
1687    int i, len;
1688    const int C = CHANNELS(st->channels);
1689    int offset;
1690    celt_sig *out_mem[2];
1691    celt_sig *decode_mem[2];
1692    celt_sig *overlap_mem[2];
1693    celt_word16 *lpc;
1694    celt_word32 *out_syn[2];
1695    celt_word16 *oldBandE, *oldLogE2, *backgroundLogE;
1696    SAVE_STACK;
1697    
1698    c=0; do {
1699       decode_mem[c] = st->_decode_mem + c*(DECODE_BUFFER_SIZE+st->overlap);
1700       out_mem[c] = decode_mem[c]+DECODE_BUFFER_SIZE-MAX_PERIOD;
1701       overlap_mem[c] = decode_mem[c]+DECODE_BUFFER_SIZE;
1702    } while (++c<C);
1703    lpc = (celt_word16*)(st->_decode_mem+(DECODE_BUFFER_SIZE+st->overlap)*C);
1704    oldBandE = lpc+C*LPC_ORDER;
1705    oldLogE2 = oldBandE + C*st->mode->nbEBands;
1706    backgroundLogE = oldLogE2  + C*st->mode->nbEBands;
1707
1708    out_syn[0] = out_mem[0]+MAX_PERIOD-N;
1709    if (C==2)
1710       out_syn[1] = out_mem[1]+MAX_PERIOD-N;
1711
1712    len = N+st->mode->overlap;
1713    
1714    if (st->loss_count >= 5)
1715    {
1716       VARDECL(celt_sig, freq);
1717       VARDECL(celt_norm, X);
1718       VARDECL(celt_ener, bandE);
1719       celt_uint32 seed;
1720
1721       ALLOC(freq, C*N, celt_sig); /**< Interleaved signal MDCTs */
1722       ALLOC(X, C*N, celt_norm);   /**< Interleaved normalised MDCTs */
1723       ALLOC(bandE, st->mode->nbEBands*C, celt_ener);
1724
1725       log2Amp(st->mode, st->start, st->end, bandE, backgroundLogE, C);
1726
1727       seed = st->rng;
1728       for (i=0;i<C*N;i++)
1729       {
1730             seed = lcg_rand(seed);
1731             X[i] = (celt_int32)(seed)>>20;
1732       }
1733       st->rng = seed;
1734       for (c=0;c<C;c++)
1735          for (i=0;i<st->mode->nbEBands;i++)
1736             renormalise_vector(X+N*c+(st->mode->eBands[i]<<LM), (st->mode->eBands[i+1]-st->mode->eBands[i])<<LM, Q15ONE);
1737
1738       denormalise_bands(st->mode, X, freq, bandE, st->mode->nbEBands, C, 1<<LM);
1739
1740       compute_inv_mdcts(st->mode, 0, freq, out_syn, overlap_mem, C, LM);
1741    } else if (st->loss_count == 0)
1742    {
1743       celt_word16 pitch_buf[MAX_PERIOD>>1];
1744       celt_word32 tmp=0;
1745       celt_word32 mem0[2]={0,0};
1746       celt_word16 mem1[2]={0,0};
1747       int len2 = len;
1748       /* FIXME: This is a kludge */
1749       if (len2>MAX_PERIOD>>1)
1750          len2 = MAX_PERIOD>>1;
1751       pitch_downsample(out_mem, pitch_buf, MAX_PERIOD, MAX_PERIOD,
1752                        C, mem0, mem1);
1753       pitch_search(st->mode, pitch_buf+((MAX_PERIOD-len2)>>1), pitch_buf, len2,
1754                    MAX_PERIOD-len2-100, &pitch_index, &tmp, 1<<LM);
1755       pitch_index = MAX_PERIOD-len2-pitch_index;
1756       st->last_pitch_index = pitch_index;
1757    } else {
1758       pitch_index = st->last_pitch_index;
1759       fade = QCONST16(.8f,15);
1760    }
1761
1762    c=0; do {
1763       /* FIXME: This is more memory than necessary */
1764       celt_word32 e[2*MAX_PERIOD];
1765       celt_word16 exc[2*MAX_PERIOD];
1766       celt_word32 ac[LPC_ORDER+1];
1767       celt_word16 decay = 1;
1768       celt_word32 S1=0;
1769       celt_word16 mem[LPC_ORDER]={0};
1770
1771       offset = MAX_PERIOD-pitch_index;
1772       for (i=0;i<MAX_PERIOD;i++)
1773          exc[i] = ROUND16(out_mem[c][i], SIG_SHIFT);
1774
1775       if (st->loss_count == 0)
1776       {
1777          _celt_autocorr(exc, ac, st->mode->window, st->mode->overlap,
1778                         LPC_ORDER, MAX_PERIOD);
1779
1780          /* Noise floor -40 dB */
1781 #ifdef FIXED_POINT
1782          ac[0] += SHR32(ac[0],13);
1783 #else
1784          ac[0] *= 1.0001f;
1785 #endif
1786          /* Lag windowing */
1787          for (i=1;i<=LPC_ORDER;i++)
1788          {
1789             /*ac[i] *= exp(-.5*(2*M_PI*.002*i)*(2*M_PI*.002*i));*/
1790 #ifdef FIXED_POINT
1791             ac[i] -= MULT16_32_Q15(2*i*i, ac[i]);
1792 #else
1793             ac[i] -= ac[i]*(.008f*i)*(.008f*i);
1794 #endif
1795          }
1796
1797          _celt_lpc(lpc+c*LPC_ORDER, ac, LPC_ORDER);
1798       }
1799       for (i=0;i<LPC_ORDER;i++)
1800          mem[i] = ROUND16(out_mem[c][MAX_PERIOD-1-i], SIG_SHIFT);
1801       fir(exc, lpc+c*LPC_ORDER, exc, MAX_PERIOD, LPC_ORDER, mem);
1802       /*for (i=0;i<MAX_PERIOD;i++)printf("%d ", exc[i]); printf("\n");*/
1803       /* Check if the waveform is decaying (and if so how fast) */
1804       {
1805          celt_word32 E1=1, E2=1;
1806          int period;
1807          if (pitch_index <= MAX_PERIOD/2)
1808             period = pitch_index;
1809          else
1810             period = MAX_PERIOD/2;
1811          for (i=0;i<period;i++)
1812          {
1813             E1 += SHR32(MULT16_16(exc[MAX_PERIOD-period+i],exc[MAX_PERIOD-period+i]),8);
1814             E2 += SHR32(MULT16_16(exc[MAX_PERIOD-2*period+i],exc[MAX_PERIOD-2*period+i]),8);
1815          }
1816          if (E1 > E2)
1817             E1 = E2;
1818          decay = celt_sqrt(frac_div32(SHR(E1,1),E2));
1819       }
1820
1821       /* Copy excitation, taking decay into account */
1822       for (i=0;i<len+st->mode->overlap;i++)
1823       {
1824          celt_word16 tmp;
1825          if (offset+i >= MAX_PERIOD)
1826          {
1827             offset -= pitch_index;
1828             decay = MULT16_16_Q15(decay, decay);
1829          }
1830          e[i] = SHL32(EXTEND32(MULT16_16_Q15(decay, exc[offset+i])), SIG_SHIFT);
1831          tmp = ROUND16(out_mem[c][offset+i],SIG_SHIFT);
1832          S1 += SHR32(MULT16_16(tmp,tmp),8);
1833       }
1834       for (i=0;i<LPC_ORDER;i++)
1835          mem[i] = ROUND16(out_mem[c][MAX_PERIOD-1-i], SIG_SHIFT);
1836       for (i=0;i<len+st->mode->overlap;i++)
1837          e[i] = MULT16_32_Q15(fade, e[i]);
1838       iir(e, lpc+c*LPC_ORDER, e, len+st->mode->overlap, LPC_ORDER, mem);
1839
1840       {
1841          celt_word32 S2=0;
1842          for (i=0;i<len+overlap;i++)
1843          {
1844             celt_word16 tmp = ROUND16(e[i],SIG_SHIFT);
1845             S2 += SHR32(MULT16_16(tmp,tmp),8);
1846          }
1847          /* This checks for an "explosion" in the synthesis */
1848 #ifdef FIXED_POINT
1849          if (!(S1 > SHR32(S2,2)))
1850 #else
1851          /* Float test is written this way to catch NaNs at the same time */
1852          if (!(S1 > 0.2f*S2))
1853 #endif
1854          {
1855             for (i=0;i<len+overlap;i++)
1856                e[i] = 0;
1857          } else if (S1 < S2)
1858          {
1859             celt_word16 ratio = celt_sqrt(frac_div32(SHR32(S1,1)+1,S2+1));
1860             for (i=0;i<len+overlap;i++)
1861                e[i] = MULT16_32_Q15(ratio, e[i]);
1862          }
1863       }
1864
1865 #ifdef ENABLE_POSTFILTER
1866       /* Apply post-filter to the MDCT overlap of the previous frame */
1867       comb_filter(out_mem[c]+MAX_PERIOD, out_mem[c]+MAX_PERIOD, st->postfilter_period, st->postfilter_period, st->overlap, C,
1868                   st->postfilter_gain, st->postfilter_gain, st->postfilter_tapset, st->postfilter_tapset,
1869                   NULL, 0);
1870 #endif /* ENABLE_POSTFILTER */
1871
1872       for (i=0;i<MAX_PERIOD+st->mode->overlap-N;i++)
1873          out_mem[c][i] = out_mem[c][N+i];
1874
1875       /* Apply TDAC to the concealed audio so that it blends with the
1876          previous and next frames */
1877       for (i=0;i<overlap/2;i++)
1878       {
1879          celt_word32 tmp;
1880          tmp = MULT16_32_Q15(st->mode->window[i],           e[N+overlap-1-i]) +
1881                MULT16_32_Q15(st->mode->window[overlap-i-1], e[N+i          ]);
1882          out_mem[c][MAX_PERIOD+i] = MULT16_32_Q15(st->mode->window[overlap-i-1], tmp);
1883          out_mem[c][MAX_PERIOD+overlap-i-1] = MULT16_32_Q15(st->mode->window[i], tmp);
1884       }
1885       for (i=0;i<N;i++)
1886          out_mem[c][MAX_PERIOD-N+i] = e[i];
1887
1888 #ifdef ENABLE_POSTFILTER
1889       /* Apply pre-filter to the MDCT overlap for the next frame (post-filter will be applied then) */
1890       comb_filter(e, out_mem[c]+MAX_PERIOD, st->postfilter_period, st->postfilter_period, st->overlap, C,
1891                   -st->postfilter_gain, -st->postfilter_gain, st->postfilter_tapset, st->postfilter_tapset,
1892                   NULL, 0);
1893 #endif /* ENABLE_POSTFILTER */
1894       for (i=0;i<overlap;i++)
1895          out_mem[c][MAX_PERIOD+i] = e[i];
1896    } while (++c<C);
1897
1898    deemphasis(out_syn, pcm, N, C, st->mode->preemph, st->preemph_memD);
1899    
1900    st->loss_count++;
1901
1902    RESTORE_STACK;
1903 }
1904
1905 #ifdef FIXED_POINT
1906 int celt_decode_with_ec(CELTDecoder * restrict st, const unsigned char *data, int len, celt_int16 * restrict pcm, int frame_size, ec_dec *dec)
1907 {
1908 #else
1909 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)
1910 {
1911 #endif
1912    int c, i, N;
1913    int spread_decision;
1914    int bits;
1915    ec_dec _dec;
1916    ec_byte_buffer buf;
1917    VARDECL(celt_sig, freq);
1918    VARDECL(celt_norm, X);
1919    VARDECL(celt_ener, bandE);
1920    VARDECL(celt_word16, oldLogE);
1921    VARDECL(int, fine_quant);
1922    VARDECL(int, pulses);
1923    VARDECL(int, offsets);
1924    VARDECL(int, fine_priority);
1925    VARDECL(int, tf_res);
1926    VARDECL(unsigned char, collapse_masks);
1927    celt_sig *out_mem[2];
1928    celt_sig *decode_mem[2];
1929    celt_sig *overlap_mem[2];
1930    celt_sig *out_syn[2];
1931    celt_word16 *lpc;
1932    celt_word16 *oldBandE, *oldLogE2, *backgroundLogE;
1933
1934    int shortBlocks;
1935    int isTransient;
1936    int intra_ener;
1937    const int C = CHANNELS(st->channels);
1938    int LM, M;
1939    int effEnd;
1940    int codedBands;
1941    int alloc_trim;
1942    int postfilter_pitch;
1943    celt_word16 postfilter_gain;
1944    int intensity=0;
1945    int dual_stereo=0;
1946    celt_int32 total_bits;
1947    celt_int32 tell;
1948    int dynalloc_logp;
1949    int postfilter_tapset;
1950    int anti_collapse_rsv;
1951    int anti_collapse_on=0;
1952    SAVE_STACK;
1953
1954    if (pcm==NULL)
1955       return CELT_BAD_ARG;
1956
1957    for (LM=0;LM<4;LM++)
1958       if (st->mode->shortMdctSize<<LM==frame_size)
1959          break;
1960    if (LM>=MAX_CONFIG_SIZES)
1961       return CELT_BAD_ARG;
1962    M=1<<LM;
1963
1964    c=0; do {
1965       decode_mem[c] = st->_decode_mem + c*(DECODE_BUFFER_SIZE+st->overlap);
1966       out_mem[c] = decode_mem[c]+DECODE_BUFFER_SIZE-MAX_PERIOD;
1967       overlap_mem[c] = decode_mem[c]+DECODE_BUFFER_SIZE;
1968    } while (++c<C);
1969    lpc = (celt_word16*)(st->_decode_mem+(DECODE_BUFFER_SIZE+st->overlap)*C);
1970    oldBandE = lpc+C*LPC_ORDER;
1971    oldLogE2 = oldBandE + C*st->mode->nbEBands;
1972    backgroundLogE = oldLogE2  + C*st->mode->nbEBands;
1973
1974    N = M*st->mode->shortMdctSize;
1975
1976    effEnd = st->end;
1977    if (effEnd > st->mode->effEBands)
1978       effEnd = st->mode->effEBands;
1979
1980    ALLOC(freq, C*N, celt_sig); /**< Interleaved signal MDCTs */
1981    ALLOC(X, C*N, celt_norm);   /**< Interleaved normalised MDCTs */
1982    ALLOC(bandE, st->mode->nbEBands*C, celt_ener);
1983    c=0; do
1984       for (i=0;i<M*st->mode->eBands[st->start];i++)
1985          X[c*N+i] = 0;
1986    while (++c<C);
1987    c=0; do   
1988       for (i=M*st->mode->eBands[effEnd];i<N;i++)
1989          X[c*N+i] = 0;
1990    while (++c<C);
1991
1992    if (data == NULL)
1993    {
1994       celt_decode_lost(st, pcm, N, LM);
1995       RESTORE_STACK;
1996       return CELT_OK;
1997    }
1998    if (len<0) {
1999      RESTORE_STACK;
2000      return CELT_BAD_ARG;
2001    }
2002    
2003    if (dec == NULL)
2004    {
2005       ec_byte_readinit(&buf,(unsigned char*)data,len);
2006       ec_dec_init(&_dec,&buf);
2007       dec = &_dec;
2008    }
2009
2010    total_bits = len*8;
2011    tell = ec_dec_tell(dec, 0);
2012
2013    postfilter_gain = 0;
2014    postfilter_pitch = 0;
2015    postfilter_tapset = 0;
2016    if (st->start==0 && tell+17 <= total_bits)
2017    {
2018       if(ec_dec_bit_logp(dec, 1))
2019       {
2020 #ifdef ENABLE_POSTFILTER
2021          int qg, octave;
2022          octave = ec_dec_uint(dec, 6);
2023          postfilter_pitch = (16<<octave)+ec_dec_bits(dec, 4+octave)-1;
2024          qg = ec_dec_bits(dec, 2);
2025          postfilter_tapset = ec_dec_icdf(dec, tapset_icdf, 2);
2026          postfilter_gain = QCONST16(.125f,15)*(qg+2);
2027 #else /* ENABLE_POSTFILTER */
2028          RESTORE_STACK;
2029          return CELT_CORRUPTED_DATA;
2030 #endif /* ENABLE_POSTFILTER */
2031       }
2032       tell = ec_dec_tell(dec, 0);
2033    }
2034
2035    if (LM > 0 && tell+3 <= total_bits)
2036    {
2037       isTransient = ec_dec_bit_logp(dec, 3);
2038       tell = ec_dec_tell(dec, 0);
2039    }
2040    else
2041       isTransient = 0;
2042
2043    if (isTransient)
2044       shortBlocks = M;
2045    else
2046       shortBlocks = 0;
2047
2048    ALLOC(oldLogE, C*st->mode->nbEBands, celt_word16);
2049    for (i=0;i<C*st->mode->nbEBands;i++)
2050       oldLogE[i] = oldBandE[i];
2051
2052    /* Decode the global flags (first symbols in the stream) */
2053    intra_ener = tell+3<=total_bits ? ec_dec_bit_logp(dec, 3) : 0;
2054    /* Get band energies */
2055    unquant_coarse_energy(st->mode, st->start, st->end, bandE, oldBandE,
2056          intra_ener, dec, C, LM);
2057
2058    ALLOC(tf_res, st->mode->nbEBands, int);
2059    tf_decode(st->start, st->end, C, isTransient, tf_res, LM, dec);
2060
2061    tell = ec_dec_tell(dec, 0);
2062    spread_decision = SPREAD_NORMAL;
2063    if (tell+4 <= total_bits)
2064       spread_decision = ec_dec_icdf(dec, spread_icdf, 5);
2065
2066    ALLOC(pulses, st->mode->nbEBands, int);
2067    ALLOC(offsets, st->mode->nbEBands, int);
2068    ALLOC(fine_priority, st->mode->nbEBands, int);
2069
2070    dynalloc_logp = 6;
2071    total_bits<<=BITRES;
2072    tell = ec_dec_tell(dec, BITRES);
2073    for (i=st->start;i<st->end;i++)
2074    {
2075       int width, quanta;
2076       int dynalloc_loop_logp;
2077       int boost;
2078       width = C*(st->mode->eBands[i+1]-st->mode->eBands[i])<<LM;
2079       /* quanta is 6 bits, but no more than 1 bit/sample
2080          and no less than 1/8 bit/sample */
2081       quanta = IMIN(width<<BITRES, IMAX(6<<BITRES, width));
2082       dynalloc_loop_logp = dynalloc_logp;
2083       boost = 0;
2084       while (tell+(dynalloc_loop_logp<<BITRES) < total_bits &&
2085             boost < (64<<LM)*(C<<BITRES))
2086       {
2087          int flag;
2088          flag = ec_dec_bit_logp(dec, dynalloc_loop_logp);
2089          tell = ec_dec_tell(dec, BITRES);
2090          if (!flag)
2091             break;
2092          boost += quanta;
2093          total_bits -= quanta;
2094          dynalloc_loop_logp = 1;
2095       }
2096       offsets[i] = boost;
2097       /* Making dynalloc more likely */
2098       if (boost>0)
2099          dynalloc_logp = IMAX(2, dynalloc_logp-1);
2100    }
2101
2102    ALLOC(fine_quant, st->mode->nbEBands, int);
2103    alloc_trim = tell+(6<<BITRES) <= total_bits ?
2104          ec_dec_icdf(dec, trim_icdf, 7) : 5;
2105
2106    bits = (len*8<<BITRES) - ec_dec_tell(dec, BITRES) - 1;
2107    anti_collapse_rsv = isTransient&&LM>=2&&bits>=(LM+2<<BITRES) ? (1<<BITRES) : 0;
2108    bits -= anti_collapse_rsv;
2109    codedBands = compute_allocation(st->mode, st->start, st->end, offsets,
2110          alloc_trim, &intensity, &dual_stereo, bits, pulses, fine_quant,
2111          fine_priority, C, LM, dec, 0, 0);
2112    
2113    unquant_fine_energy(st->mode, st->start, st->end, bandE, oldBandE, fine_quant, dec, C);
2114
2115    /* Decode fixed codebook */
2116    ALLOC(collapse_masks, st->mode->nbEBands, unsigned char);
2117    quant_all_bands(0, st->mode, st->start, st->end, X, C==2 ? X+N : NULL, collapse_masks,
2118          NULL, pulses, shortBlocks, spread_decision, dual_stereo, intensity, tf_res, 1,
2119          len*8, dec, LM, codedBands, &st->rng);
2120
2121    if (anti_collapse_rsv > 0)
2122    {
2123       anti_collapse_on = ec_dec_bits(dec, 1);
2124    }
2125
2126    unquant_energy_finalise(st->mode, st->start, st->end, bandE, oldBandE,
2127          fine_quant, fine_priority, len*8-ec_dec_tell(dec, 0), dec, C);
2128
2129    if (anti_collapse_on)
2130       anti_collapse(st->mode, X, collapse_masks, LM, C, N,
2131             st->start, st->end, oldBandE, oldLogE, oldLogE2, pulses, st->rng);
2132
2133    log2Amp(st->mode, st->start, st->end, bandE, oldBandE, C);
2134
2135    /* Synthesis */
2136    denormalise_bands(st->mode, X, freq, bandE, effEnd, C, M);
2137
2138    CELT_MOVE(decode_mem[0], decode_mem[0]+N, DECODE_BUFFER_SIZE-N);
2139    if (C==2)
2140       CELT_MOVE(decode_mem[1], decode_mem[1]+N, DECODE_BUFFER_SIZE-N);
2141
2142    c=0; do
2143       for (i=0;i<M*st->mode->eBands[st->start];i++)
2144          freq[c*N+i] = 0;
2145    while (++c<C);
2146    c=0; do
2147       for (i=M*st->mode->eBands[effEnd];i<N;i++)
2148          freq[c*N+i] = 0;
2149    while (++c<C);
2150
2151    out_syn[0] = out_mem[0]+MAX_PERIOD-N;
2152    if (C==2)
2153       out_syn[1] = out_mem[1]+MAX_PERIOD-N;
2154
2155    /* Compute inverse MDCTs */
2156    compute_inv_mdcts(st->mode, shortBlocks, freq, out_syn, overlap_mem, C, LM);
2157
2158 #ifdef ENABLE_POSTFILTER
2159    c=0; do {
2160       st->postfilter_period=IMAX(st->postfilter_period, COMBFILTER_MINPERIOD);
2161       st->postfilter_period_old=IMAX(st->postfilter_period_old, COMBFILTER_MINPERIOD);
2162       if (LM!=0)
2163       {
2164          comb_filter(out_syn[c], out_syn[c], st->postfilter_period, st->postfilter_period, st->overlap, C,
2165                st->postfilter_gain, st->postfilter_gain, st->postfilter_tapset, st->postfilter_tapset,
2166                NULL, 0);
2167          comb_filter(out_syn[c]+st->overlap, out_syn[c]+st->overlap, st->postfilter_period, postfilter_pitch, N-st->overlap, C,
2168                st->postfilter_gain, postfilter_gain, st->postfilter_tapset, postfilter_tapset,
2169                st->mode->window, st->mode->overlap);
2170       } else {
2171          comb_filter(out_syn[c], out_syn[c], st->postfilter_period_old, st->postfilter_period, N-st->overlap, C,
2172                st->postfilter_gain_old, st->postfilter_gain, st->postfilter_tapset_old, st->postfilter_tapset,
2173                st->mode->window, st->mode->overlap);
2174       }
2175    } while (++c<C);
2176    st->postfilter_period_old = st->postfilter_period;
2177    st->postfilter_gain_old = st->postfilter_gain;
2178    st->postfilter_tapset_old = st->postfilter_tapset;
2179    st->postfilter_period = postfilter_pitch;
2180    st->postfilter_gain = postfilter_gain;
2181    st->postfilter_tapset = postfilter_tapset;
2182 #endif /* ENABLE_POSTFILTER */
2183
2184    /* In case start or end were to change */
2185    c=0; do
2186    {
2187       for (i=0;i<st->start;i++)
2188          oldBandE[c*st->mode->nbEBands+i]=0;
2189       for (i=st->end;i<st->mode->nbEBands;i++)
2190          oldBandE[c*st->mode->nbEBands+i]=0;
2191    } while (++c<C);
2192    for (i=0;i<C*st->mode->nbEBands;i++)
2193       oldLogE2[i] = oldLogE[i];
2194    for (i=0;i<C*st->mode->nbEBands;i++)
2195       backgroundLogE[i] = MIN16(backgroundLogE[i] + M*QCONST16(0.001f,DB_SHIFT), oldLogE[i]);
2196
2197    st->rng = dec->rng;
2198
2199    deemphasis(out_syn, pcm, N, C, st->mode->preemph, st->preemph_memD);
2200    st->loss_count = 0;
2201    RESTORE_STACK;
2202    if (ec_dec_tell(dec,0) > 8*len || ec_dec_get_error(dec))
2203       return CELT_CORRUPTED_DATA;
2204    else
2205       return CELT_OK;
2206 }
2207
2208 #ifdef FIXED_POINT
2209 #ifndef DISABLE_FLOAT_API
2210 int celt_decode_with_ec_float(CELTDecoder * restrict st, const unsigned char *data, int len, float * restrict pcm, int frame_size, ec_dec *dec)
2211 {
2212    int j, ret, C, N, LM, M;
2213    VARDECL(celt_int16, out);
2214    SAVE_STACK;
2215
2216    if (pcm==NULL)
2217       return CELT_BAD_ARG;
2218
2219    for (LM=0;LM<4;LM++)
2220       if (st->mode->shortMdctSize<<LM==frame_size)
2221          break;
2222    if (LM>=MAX_CONFIG_SIZES)
2223       return CELT_BAD_ARG;
2224    M=1<<LM;
2225
2226    C = CHANNELS(st->channels);
2227    N = M*st->mode->shortMdctSize;
2228    
2229    ALLOC(out, C*N, celt_int16);
2230    ret=celt_decode_with_ec(st, data, len, out, frame_size, dec);
2231    if (ret==0)
2232       for (j=0;j<C*N;j++)
2233          pcm[j]=out[j]*(1.f/32768.f);
2234      
2235    RESTORE_STACK;
2236    return ret;
2237 }
2238 #endif /*DISABLE_FLOAT_API*/
2239 #else
2240 int celt_decode_with_ec(CELTDecoder * restrict st, const unsigned char *data, int len, celt_int16 * restrict pcm, int frame_size, ec_dec *dec)
2241 {
2242    int j, ret, C, N, LM, M;
2243    VARDECL(celt_sig, out);
2244    SAVE_STACK;
2245
2246    if (pcm==NULL)
2247       return CELT_BAD_ARG;
2248
2249    for (LM=0;LM<4;LM++)
2250       if (st->mode->shortMdctSize<<LM==frame_size)
2251          break;
2252    if (LM>=MAX_CONFIG_SIZES)
2253       return CELT_BAD_ARG;
2254    M=1<<LM;
2255
2256    C = CHANNELS(st->channels);
2257    N = M*st->mode->shortMdctSize;
2258    ALLOC(out, C*N, celt_sig);
2259
2260    ret=celt_decode_with_ec_float(st, data, len, out, frame_size, dec);
2261
2262    if (ret==0)
2263       for (j=0;j<C*N;j++)
2264          pcm[j] = FLOAT2INT16 (out[j]);
2265    
2266    RESTORE_STACK;
2267    return ret;
2268 }
2269 #endif
2270
2271 int celt_decode(CELTDecoder * restrict st, const unsigned char *data, int len, celt_int16 * restrict pcm, int frame_size)
2272 {
2273    return celt_decode_with_ec(st, data, len, pcm, frame_size, NULL);
2274 }
2275
2276 #ifndef DISABLE_FLOAT_API
2277 int celt_decode_float(CELTDecoder * restrict st, const unsigned char *data, int len, float * restrict pcm, int frame_size)
2278 {
2279    return celt_decode_with_ec_float(st, data, len, pcm, frame_size, NULL);
2280 }
2281 #endif /* DISABLE_FLOAT_API */
2282
2283 int celt_decoder_ctl(CELTDecoder * restrict st, int request, ...)
2284 {
2285    va_list ap;
2286
2287    va_start(ap, request);
2288    switch (request)
2289    {
2290       case CELT_GET_MODE_REQUEST:
2291       {
2292          const CELTMode ** value = va_arg(ap, const CELTMode**);
2293          if (value==0)
2294             goto bad_arg;
2295          *value=st->mode;
2296       }
2297       break;
2298       case CELT_SET_START_BAND_REQUEST:
2299       {
2300          celt_int32 value = va_arg(ap, celt_int32);
2301          if (value<0 || value>=st->mode->nbEBands)
2302             goto bad_arg;
2303          st->start = value;
2304       }
2305       break;
2306       case CELT_SET_END_BAND_REQUEST:
2307       {
2308          celt_int32 value = va_arg(ap, celt_int32);
2309          if (value<0 || value>=st->mode->nbEBands)
2310             goto bad_arg;
2311          st->end = value;
2312       }
2313       break;
2314       case CELT_RESET_STATE:
2315       {
2316          CELT_MEMSET((char*)&st->DECODER_RESET_START, 0,
2317                celt_decoder_get_size(st->mode, st->channels)-
2318                ((char*)&st->DECODER_RESET_START - (char*)st));
2319       }
2320       break;
2321       default:
2322          goto bad_request;
2323    }
2324    va_end(ap);
2325    return CELT_OK;
2326 bad_arg:
2327    va_end(ap);
2328    return CELT_BAD_ARG;
2329 bad_request:
2330       va_end(ap);
2331   return CELT_UNIMPLEMENTED;
2332 }
2333
2334 const char *celt_strerror(int error)
2335 {
2336    static const char *error_strings[8] = {
2337       "success",
2338       "invalid argument",
2339       "invalid mode",
2340       "internal error",
2341       "corrupted stream",
2342       "request not implemented",
2343       "invalid state",
2344       "memory allocation failed"
2345    };
2346    if (error > 0 || error < -7)
2347       return "unknown error";
2348    else 
2349       return error_strings[-error];
2350 }
2351