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