92947557018f451e3cdeac03e6e58339dab10c2d
[opus.git] / libcelt / celt.c
1 /* (C) 2007-2008 Jean-Marc Valin, CSIRO
2 */
3 /*
4    Redistribution and use in source and binary forms, with or without
5    modification, are permitted provided that the following conditions
6    are met:
7    
8    - Redistributions of source code must retain the above copyright
9    notice, this list of conditions and the following disclaimer.
10    
11    - Redistributions in binary form must reproduce the above copyright
12    notice, this list of conditions and the following disclaimer in the
13    documentation and/or other materials provided with the distribution.
14    
15    - Neither the name of the Xiph.org Foundation nor the names of its
16    contributors may be used to endorse or promote products derived from
17    this software without specific prior written permission.
18    
19    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
20    ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
21    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
22    A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE FOUNDATION OR
23    CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
24    EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
25    PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
26    PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
27    LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
28    NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
29    SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
30 */
31
32 #ifdef HAVE_CONFIG_H
33 #include "config.h"
34 #endif
35
36 #define CELT_C
37
38 #include "os_support.h"
39 #include "mdct.h"
40 #include <math.h>
41 #include "celt.h"
42 #include "pitch.h"
43 #include "kiss_fftr.h"
44 #include "bands.h"
45 #include "modes.h"
46 #include "entcode.h"
47 #include "quant_pitch.h"
48 #include "quant_bands.h"
49 #include "psy.h"
50 #include "rate.h"
51 #include "stack_alloc.h"
52 #include "mathops.h"
53 #include "float_cast.h"
54 #include <stdarg.h>
55
56 static const celt_word16_t preemph = QCONST16(0.8f,15);
57
58 #ifdef FIXED_POINT
59 static const celt_word16_t transientWindow[16] = {
60      279,  1106,  2454,  4276,  6510,  9081, 11900, 14872,
61    17896, 20868, 23687, 26258, 28492, 30314, 31662, 32489};
62 #else
63 static const float transientWindow[16] = {
64    0.0085135, 0.0337639, 0.0748914, 0.1304955, 0.1986827, 0.2771308, 0.3631685, 0.4538658,
65    0.5461342, 0.6368315, 0.7228692, 0.8013173, 0.8695045, 0.9251086, 0.9662361, 0.9914865};
66 #endif
67
68    
69 /** Encoder state 
70  @brief Encoder state
71  */
72 struct CELTEncoder {
73    const CELTMode *mode;     /**< Mode used by the encoder */
74    int frame_size;
75    int block_size;
76    int overlap;
77    int channels;
78    
79    int pitch_enabled;
80    
81    ec_byte_buffer buf;
82    ec_enc         enc;
83
84    celt_word16_t * restrict preemph_memE; /* Input is 16-bit, so why bother with 32 */
85    celt_sig_t    * restrict preemph_memD;
86
87    celt_sig_t *in_mem;
88    celt_sig_t *out_mem;
89
90    celt_word16_t *oldBandE;
91 #ifdef EXP_PSY
92    celt_word16_t *psy_mem;
93    struct PsyDecay psy;
94 #endif
95 };
96
97 CELTEncoder *celt_encoder_create(const CELTMode *mode)
98 {
99    int N, C;
100    CELTEncoder *st;
101
102    if (check_mode(mode) != CELT_OK)
103       return NULL;
104
105    N = mode->mdctSize;
106    C = mode->nbChannels;
107    st = celt_alloc(sizeof(CELTEncoder));
108    
109    st->mode = mode;
110    st->frame_size = N;
111    st->block_size = N;
112    st->overlap = mode->overlap;
113
114    st->pitch_enabled = 1;
115    
116    ec_byte_writeinit(&st->buf);
117    ec_enc_init(&st->enc,&st->buf);
118
119    st->in_mem = celt_alloc(st->overlap*C*sizeof(celt_sig_t));
120    st->out_mem = celt_alloc((MAX_PERIOD+st->overlap)*C*sizeof(celt_sig_t));
121
122    st->oldBandE = (celt_word16_t*)celt_alloc(C*mode->nbEBands*sizeof(celt_word16_t));
123
124    st->preemph_memE = (celt_word16_t*)celt_alloc(C*sizeof(celt_word16_t));
125    st->preemph_memD = (celt_sig_t*)celt_alloc(C*sizeof(celt_sig_t));
126
127 #ifdef EXP_PSY
128    st->psy_mem = celt_alloc(MAX_PERIOD*sizeof(celt_word16_t));
129    psydecay_init(&st->psy, MAX_PERIOD/2, st->mode->Fs);
130 #endif
131
132    return st;
133 }
134
135 void celt_encoder_destroy(CELTEncoder *st)
136 {
137    if (st == NULL)
138    {
139       celt_warning("NULL passed to celt_encoder_destroy");
140       return;
141    }
142    if (check_mode(st->mode) != CELT_OK)
143       return;
144
145    ec_byte_writeclear(&st->buf);
146
147    celt_free(st->in_mem);
148    celt_free(st->out_mem);
149    
150    celt_free(st->oldBandE);
151    
152    celt_free(st->preemph_memE);
153    celt_free(st->preemph_memD);
154    
155 #ifdef EXP_PSY
156    celt_free (st->psy_mem);
157    psydecay_clear(&st->psy);
158 #endif
159    
160    celt_free(st);
161 }
162
163 static inline celt_int16_t FLOAT2INT16(float x)
164 {
165    x = x*32768.;
166    x = MAX32(x, -32768);
167    x = MIN32(x, 32767);
168    return (celt_int16_t)float2int(x);
169 }
170
171 static inline celt_word16_t SIG2WORD16(celt_sig_t x)
172 {
173 #ifdef FIXED_POINT
174    x = PSHR32(x, SIG_SHIFT);
175    x = MAX32(x, -32768);
176    x = MIN32(x, 32767);
177    return EXTRACT16(x);
178 #else
179    return (celt_word16_t)x;
180 #endif
181 }
182
183 static int transient_analysis(celt_word32_t *in, int len, int C, int *transient_time, int *transient_shift)
184 {
185    int c, i, n;
186    celt_word32_t ratio;
187    /* FIXME: Remove the floats here */
188    VARDECL(celt_word32_t, begin);
189    SAVE_STACK;
190    ALLOC(begin, len, celt_word32_t);
191    for (i=0;i<len;i++)
192       begin[i] = ABS32(SHR32(in[C*i],SIG_SHIFT));
193    for (c=1;c<C;c++)
194    {
195       for (i=0;i<len;i++)
196          begin[i] = MAX32(begin[i], ABS32(SHR32(in[C*i+c],SIG_SHIFT)));
197    }
198    for (i=1;i<len;i++)
199       begin[i] = MAX32(begin[i-1],begin[i]);
200    n = -1;
201    for (i=8;i<len-8;i++)
202    {
203       if (begin[i] < MULT16_32_Q15(QCONST16(.2f,15),begin[len-1]))
204          n=i;
205    }
206    if (n<32)
207    {
208       n = -1;
209       ratio = 0;
210    } else {
211       ratio = DIV32(begin[len-1],1+begin[n-16]);
212    }
213    /*printf ("%d %f\n", n, ratio*ratio);*/
214    if (ratio < 0)
215       ratio = 0;
216    if (ratio > 1000)
217       ratio = 1000;
218    ratio *= ratio;
219    if (ratio < 50)
220       *transient_shift = 0;
221    else if (ratio < 256)
222       *transient_shift = 1;
223    else if (ratio < 4096)
224       *transient_shift = 2;
225    else
226       *transient_shift = 3;
227    *transient_time = n;
228    
229    RESTORE_STACK;
230    return ratio > 20;
231 }
232
233 /** Apply window and compute the MDCT for all sub-frames and all channels in a frame */
234 static void compute_mdcts(const CELTMode *mode, int shortBlocks, celt_sig_t * restrict in, celt_sig_t * restrict out)
235 {
236    const int C = CHANNELS(mode);
237    if (C==1 && !shortBlocks)
238    {
239       const mdct_lookup *lookup = MDCT(mode);
240       const int overlap = OVERLAP(mode);
241       mdct_forward(lookup, in, out, mode->window, overlap);
242    } else if (!shortBlocks) {
243       const mdct_lookup *lookup = MDCT(mode);
244       const int overlap = OVERLAP(mode);
245       const int N = FRAMESIZE(mode);
246       int c;
247       VARDECL(celt_word32_t, x);
248       VARDECL(celt_word32_t, tmp);
249       SAVE_STACK;
250       ALLOC(x, N+overlap, celt_word32_t);
251       ALLOC(tmp, N, celt_word32_t);
252       for (c=0;c<C;c++)
253       {
254          int j;
255          for (j=0;j<N+overlap;j++)
256             x[j] = in[C*j+c];
257          mdct_forward(lookup, x, tmp, mode->window, overlap);
258          /* Interleaving the sub-frames */
259          for (j=0;j<N;j++)
260             out[C*j+c] = tmp[j];
261       }
262       RESTORE_STACK;
263    } else {
264       const mdct_lookup *lookup = &mode->shortMdct;
265       const int overlap = mode->overlap;
266       const int N = mode->shortMdctSize;
267       int b, c;
268       VARDECL(celt_word32_t, x);
269       VARDECL(celt_word32_t, tmp);
270       SAVE_STACK;
271       ALLOC(x, N+overlap, celt_word32_t);
272       ALLOC(tmp, N, celt_word32_t);
273       for (c=0;c<C;c++)
274       {
275          int B = mode->nbShortMdcts;
276          for (b=0;b<B;b++)
277          {
278             int j;
279             for (j=0;j<N+overlap;j++)
280                x[j] = in[C*(b*N+j)+c];
281             mdct_forward(lookup, x, tmp, mode->window, overlap);
282             /* Interleaving the sub-frames */
283             for (j=0;j<N;j++)
284                out[C*(j*B+b)+c] = tmp[j];
285          }
286       }
287       RESTORE_STACK;
288    }
289 }
290
291 /** Compute the IMDCT and apply window for all sub-frames and all channels in a frame */
292 static void compute_inv_mdcts(const CELTMode *mode, int shortBlocks, celt_sig_t *X, int transient_time, int transient_shift, celt_sig_t * restrict out_mem)
293 {
294    int c, N4;
295    const int C = CHANNELS(mode);
296    const int N = FRAMESIZE(mode);
297    const int overlap = OVERLAP(mode);
298    N4 = (N-overlap)>>1;
299    for (c=0;c<C;c++)
300    {
301       int j;
302       if (transient_shift==0 && C==1 && !shortBlocks) {
303          const mdct_lookup *lookup = MDCT(mode);
304          mdct_backward(lookup, X, out_mem+C*(MAX_PERIOD-N-N4), mode->window, overlap);
305       } else if (!shortBlocks) {
306          const mdct_lookup *lookup = MDCT(mode);
307          VARDECL(celt_word32_t, x);
308          VARDECL(celt_word32_t, tmp);
309          SAVE_STACK;
310          ALLOC(x, 2*N, celt_word32_t);
311          ALLOC(tmp, N, celt_word32_t);
312          /* De-interleaving the sub-frames */
313          for (j=0;j<N;j++)
314             tmp[j] = X[C*j+c];
315          /* Prevents problems from the imdct doing the overlap-add */
316          CELT_MEMSET(x+N4, 0, N);
317          mdct_backward(lookup, tmp, x, mode->window, overlap);
318          celt_assert(transient_shift == 0);
319          /* The first and last part would need to be set to zero if we actually
320             wanted to use them. */
321          for (j=0;j<overlap;j++)
322             out_mem[C*(MAX_PERIOD-N)+C*j+c] += x[j+N4];
323          for (j=0;j<overlap;j++)
324             out_mem[C*(MAX_PERIOD)+C*(overlap-j-1)+c] = x[2*N-j-N4-1];
325          for (j=0;j<2*N4;j++)
326             out_mem[C*(MAX_PERIOD-N)+C*(j+overlap)+c] = x[j+N4+overlap];
327          RESTORE_STACK;
328       } else {
329          int b;
330          const int N2 = mode->shortMdctSize;
331          const int B = mode->nbShortMdcts;
332          const mdct_lookup *lookup = &mode->shortMdct;
333          VARDECL(celt_word32_t, x);
334          VARDECL(celt_word32_t, tmp);
335          SAVE_STACK;
336          ALLOC(x, 2*N, celt_word32_t);
337          ALLOC(tmp, N, celt_word32_t);
338          /* Prevents problems from the imdct doing the overlap-add */
339          CELT_MEMSET(x+N4, 0, N2);
340          for (b=0;b<B;b++)
341          {
342             /* De-interleaving the sub-frames */
343             for (j=0;j<N2;j++)
344                tmp[j] = X[C*(j*B+b)+c];
345             mdct_backward(lookup, tmp, x+N4+N2*b, mode->window, overlap);
346          }
347          if (transient_shift > 0)
348          {
349 #ifdef FIXED_POINT
350             for (j=0;j<16;j++)
351                x[N4+transient_time+j-16] = MULT16_32_Q15(SHR16(Q15_ONE-transientWindow[j],transient_shift)+transientWindow[j], SHL32(x[N4+transient_time+j-16],transient_shift));
352             for (j=transient_time;j<N+overlap;j++)
353                x[N4+j] = SHL32(x[N4+j], transient_shift);
354 #else
355             for (j=0;j<16;j++)
356                x[N4+transient_time+j-16] *= 1+transientWindow[j]*((1<<transient_shift)-1);
357             for (j=transient_time;j<N+overlap;j++)
358                x[N4+j] *= 1<<transient_shift;
359 #endif
360          }
361          /* The first and last part would need to be set to zero if we actually
362          wanted to use them. */
363          for (j=0;j<overlap;j++)
364             out_mem[C*(MAX_PERIOD-N)+C*j+c] += x[j+N4];
365          for (j=0;j<overlap;j++)
366             out_mem[C*(MAX_PERIOD)+C*(overlap-j-1)+c] = x[2*N-j-N4-1];
367          for (j=0;j<2*N4;j++)
368             out_mem[C*(MAX_PERIOD-N)+C*(j+overlap)+c] = x[j+N4+overlap];
369          RESTORE_STACK;
370       }
371    }
372 }
373
374 #ifdef FIXED_POINT
375 int celt_encode(CELTEncoder * restrict st, const celt_int16_t * pcm, celt_int16_t * optional_synthesis, unsigned char *compressed, int nbCompressedBytes)
376 {
377 #else
378 int celt_encode_float(CELTEncoder * restrict st, const celt_sig_t * pcm, celt_sig_t * optional_synthesis, unsigned char *compressed, int nbCompressedBytes)
379 {
380 #endif
381    int i, c, N, N4;
382    int has_pitch;
383    int pitch_index;
384    int bits;
385    celt_word32_t curr_power, pitch_power=0;
386    VARDECL(celt_sig_t, in);
387    VARDECL(celt_sig_t, freq);
388    VARDECL(celt_norm_t, X);
389    VARDECL(celt_norm_t, P);
390    VARDECL(celt_ener_t, bandE);
391    VARDECL(celt_pgain_t, gains);
392    VARDECL(int, stereo_mode);
393    VARDECL(int, fine_quant);
394    VARDECL(celt_word16_t, error);
395    VARDECL(int, pulses);
396    VARDECL(int, offsets);
397 #ifdef EXP_PSY
398    VARDECL(celt_word32_t, mask);
399 #endif
400    int shortBlocks=0;
401    int transient_time;
402    int transient_shift;
403    const int C = CHANNELS(st->mode);
404    SAVE_STACK;
405
406    if (check_mode(st->mode) != CELT_OK)
407       return CELT_INVALID_MODE;
408
409    N = st->block_size;
410    N4 = (N-st->overlap)>>1;
411    ALLOC(in, 2*C*N-2*C*N4, celt_sig_t);
412
413    CELT_COPY(in, st->in_mem, C*st->overlap);
414    for (c=0;c<C;c++)
415    {
416       const celt_word16_t * restrict pcmp = pcm+c;
417       celt_sig_t * restrict inp = in+C*st->overlap+c;
418       for (i=0;i<N;i++)
419       {
420          /* Apply pre-emphasis */
421          celt_sig_t tmp = SCALEIN(SHL32(EXTEND32(*pcmp), SIG_SHIFT));
422          *inp = SUB32(tmp, SHR32(MULT16_16(preemph,st->preemph_memE[c]),3));
423          st->preemph_memE[c] = SCALEIN(*pcmp);
424          inp += C;
425          pcmp += C;
426       }
427    }
428    CELT_COPY(st->in_mem, in+C*(2*N-2*N4-st->overlap), C*st->overlap);
429    
430    if (st->mode->nbShortMdcts > 1)
431    {
432       if (transient_analysis(in, N+st->overlap, C, &transient_time, &transient_shift))
433       {
434 #ifndef FIXED_POINT
435          float gain_1;
436 #endif
437          ec_enc_bits(&st->enc, 0, 1); //Pitch off
438          ec_enc_bits(&st->enc, 1, 1); //Transient on
439          ec_enc_bits(&st->enc, transient_shift, 2);
440          if (transient_shift)
441             ec_enc_uint(&st->enc, transient_time, N+st->overlap);
442          if (transient_shift)
443          {
444 #ifdef FIXED_POINT
445             for (c=0;c<C;c++)
446                for (i=0;i<16;i++)
447                   in[C*(transient_time+i-16)+c] = MULT16_32_Q15(EXTRACT16(SHR32(celt_rcp(Q15ONE+MULT16_16(transientWindow[i],((1<<transient_shift)-1))),1)), in[C*(transient_time+i-16)+c]);
448             for (c=0;c<C;c++)
449                for (i=transient_time;i<N+st->overlap;i++)
450                   in[C*i+c] = SHR32(in[C*i+c], transient_shift);
451 #else
452             for (c=0;c<C;c++)
453                for (i=0;i<16;i++)
454                   in[C*(transient_time+i-16)+c] /= 1+transientWindow[i]*((1<<transient_shift)-1);
455             gain_1 = 1./(1<<transient_shift);
456             for (c=0;c<C;c++)
457                for (i=transient_time;i<N+st->overlap;i++)
458                   in[C*i+c] *= gain_1;
459 #endif
460          }
461          shortBlocks = 1;
462       } else {
463          transient_time = -1;
464          transient_shift = 0;
465          shortBlocks = 0;
466       }
467    } else {
468       transient_time = -1;
469       transient_shift = 0;
470       shortBlocks = 0;
471    }
472    /* Pitch analysis: we do it early to save on the peak stack space */
473    if (st->pitch_enabled && !shortBlocks)
474    {
475 #ifdef EXP_PSY
476       VARDECL(celt_word16_t, X);
477       ALLOC(X, MAX_PERIOD/2, celt_word16_t);
478       find_spectral_pitch(st->mode, st->mode->fft, &st->mode->psy, in, st->out_mem, st->mode->window, X, 2*N-2*N4, MAX_PERIOD-(2*N-2*N4), &pitch_index);
479       compute_tonality(st->mode, X, st->psy_mem, MAX_PERIOD);
480 #else
481       find_spectral_pitch(st->mode, st->mode->fft, &st->mode->psy, in, st->out_mem, st->mode->window, NULL, 2*N-2*N4, MAX_PERIOD-(2*N-2*N4), &pitch_index);
482 #endif
483    }
484    ALLOC(freq, C*N, celt_sig_t); /**< Interleaved signal MDCTs */
485    
486    /*for (i=0;i<(B+1)*C*N;i++) printf ("%f(%d) ", in[i], i); printf ("\n");*/
487    /* Compute MDCTs */
488    compute_mdcts(st->mode, shortBlocks, in, freq);
489
490 #ifdef EXP_PSY
491    /*CELT_MOVE(st->psy_mem, st->out_mem+N, MAX_PERIOD+st->overlap-N);
492    for (i=0;i<N;i++)
493       st->psy_mem[MAX_PERIOD+st->overlap-N+i] = in[C*(st->overlap+i)];
494    for (c=1;c<C;c++)
495       for (i=0;i<N;i++)
496          st->psy_mem[MAX_PERIOD+st->overlap-N+i] += in[C*(st->overlap+i)+c];
497    */
498    ALLOC(mask, N, celt_sig_t);
499    compute_mdct_masking(&st->psy, freq, st->psy_mem, mask, C*N);
500
501    /* Invert and stretch the mask to length of X 
502       For some reason, I get better results by using the sqrt instead,
503       although there's no valid reason to. Must investigate further */
504    for (i=0;i<C*N;i++)
505       mask[i] = 1/(.1+mask[i]);
506 #endif
507    
508    /* Deferred allocation after find_spectral_pitch() to reduce the peak memory usage */
509    ALLOC(X, C*N, celt_norm_t);         /**< Interleaved normalised MDCTs */
510    ALLOC(P, C*N, celt_norm_t);         /**< Interleaved normalised pitch MDCTs*/
511    ALLOC(bandE,st->mode->nbEBands*C, celt_ener_t);
512    ALLOC(gains,st->mode->nbPBands, celt_pgain_t);
513
514    /*printf ("%f %f\n", curr_power, pitch_power);*/
515    /*int j;
516    for (j=0;j<B*N;j++)
517       printf ("%f ", X[j]);
518    for (j=0;j<B*N;j++)
519       printf ("%f ", P[j]);
520    printf ("\n");*/
521
522    /* Band normalisation */
523    compute_band_energies(st->mode, freq, bandE);
524    normalise_bands(st->mode, freq, X, bandE);
525    /*for (i=0;i<st->mode->nbEBands;i++)printf("%f ", bandE[i]);printf("\n");*/
526    /*for (i=0;i<N*B*C;i++)printf("%f ", X[i]);printf("\n");*/
527
528    /* Compute MDCTs of the pitch part */
529    if (st->pitch_enabled && !shortBlocks)
530    {
531       /* Normalise the pitch vector as well (discard the energies) */
532       VARDECL(celt_ener_t, bandEp);
533       
534       compute_mdcts(st->mode, 0, st->out_mem+pitch_index*C, freq);
535       ALLOC(bandEp, st->mode->nbEBands*st->mode->nbChannels, celt_ener_t);
536       compute_band_energies(st->mode, freq, bandEp);
537       normalise_bands(st->mode, freq, P, bandEp);
538       pitch_power = bandEp[0]+bandEp[1]+bandEp[2];
539    }
540    curr_power = bandE[0]+bandE[1]+bandE[2];
541    /* Check if we can safely use the pitch (i.e. effective gain isn't too high) */
542    if (st->pitch_enabled && !shortBlocks && (MULT16_32_Q15(QCONST16(.1f, 15),curr_power) + QCONST32(10.f,ENER_SHIFT) < pitch_power))
543    {
544       /* Simulates intensity stereo */
545       /*for (i=30;i<N*B;i++)
546          X[i*C+1] = P[i*C+1] = 0;*/
547
548       /* Pitch prediction */
549       compute_pitch_gain(st->mode, X, P, gains);
550       has_pitch = quant_pitch(gains, st->mode->nbPBands, &st->enc);
551       if (has_pitch)
552          ec_enc_uint(&st->enc, pitch_index, MAX_PERIOD-(2*N-2*N4));
553       else if (st->mode->nbShortMdcts > 1)
554          ec_enc_bits(&st->enc, 0, 1); //Transient off
555    } else {
556       if (!shortBlocks)
557       {
558          ec_enc_bits(&st->enc, 0, 1); //Pitch off
559          if (st->mode->nbShortMdcts > 1)
560            ec_enc_bits(&st->enc, 0, 1); //Transient off
561       }
562       /* No pitch, so we just pretend we found a gain of zero */
563       for (i=0;i<st->mode->nbPBands;i++)
564          gains[i] = 0;
565       for (i=0;i<C*N;i++)
566          P[i] = 0;
567    }
568
569 #ifdef STDIN_TUNING2
570    static int fine_quant[30];
571    static int pulses[30];
572    static int init=0;
573    if (!init)
574    {
575       for (i=0;i<st->mode->nbEBands;i++)
576          scanf("%d ", &fine_quant[i]);
577       for (i=0;i<st->mode->nbEBands;i++)
578          scanf("%d ", &pulses[i]);
579       init = 1;
580    }
581 #else
582    ALLOC(fine_quant, st->mode->nbEBands, int);
583    ALLOC(pulses, st->mode->nbEBands, int);
584 #endif
585    ALLOC(error, C*st->mode->nbEBands, celt_word16_t);
586    quant_coarse_energy(st->mode, bandE, st->oldBandE, nbCompressedBytes*8/3, st->mode->prob, error, &st->enc);
587    
588    ALLOC(offsets, st->mode->nbEBands, int);
589    ALLOC(stereo_mode, st->mode->nbEBands, int);
590    stereo_decision(st->mode, X, stereo_mode, st->mode->nbEBands);
591
592    for (i=0;i<st->mode->nbEBands;i++)
593       offsets[i] = 0;
594    bits = nbCompressedBytes*8 - ec_enc_tell(&st->enc, 0) - 1;
595 #ifndef STDIN_TUNING
596    compute_allocation(st->mode, offsets, stereo_mode, bits, pulses, fine_quant);
597 #endif
598    /*for (i=0;i<st->mode->nbEBands;i++)
599       printf("%d ", fine_quant[i]);
600    for (i=0;i<st->mode->nbEBands;i++)
601       printf("%d ", pulses[i]);
602    printf ("\n");*/
603    /*bits = ec_enc_tell(&st->enc, 0);
604    compute_fine_allocation(st->mode, fine_quant, (20*C+nbCompressedBytes*8/5-(ec_enc_tell(&st->enc, 0)-bits))/C);*/
605    quant_fine_energy(st->mode, bandE, st->oldBandE, error, fine_quant, &st->enc);
606
607    pitch_quant_bands(st->mode, P, gains);
608
609    /*for (i=0;i<B*N;i++) printf("%f ",P[i]);printf("\n");*/
610
611    /* Residual quantisation */
612    quant_bands(st->mode, X, P, NULL, bandE, stereo_mode, pulses, shortBlocks, nbCompressedBytes*8, &st->enc);
613    
614    if (st->pitch_enabled || optional_synthesis!=NULL)
615    {
616       if (C==2)
617          renormalise_bands(st->mode, X);
618       /* Synthesis */
619       denormalise_bands(st->mode, X, freq, bandE);
620       
621       
622       CELT_MOVE(st->out_mem, st->out_mem+C*N, C*(MAX_PERIOD+st->overlap-N));
623       
624       compute_inv_mdcts(st->mode, shortBlocks, freq, transient_time, transient_shift, st->out_mem);
625       /* De-emphasis and put everything back at the right place in the synthesis history */
626       if (optional_synthesis != NULL) {
627          for (c=0;c<C;c++)
628          {
629             int j;
630             for (j=0;j<N;j++)
631             {
632                celt_sig_t tmp = MAC16_32_Q15(st->out_mem[C*(MAX_PERIOD-N)+C*j+c],
633                                    preemph,st->preemph_memD[c]);
634                st->preemph_memD[c] = tmp;
635                optional_synthesis[C*j+c] = SCALEOUT(SIG2WORD16(tmp));
636             }
637          }
638       }
639    }
640    /*fprintf (stderr, "remaining bits after encode = %d\n", nbCompressedBytes*8-ec_enc_tell(&st->enc, 0));*/
641    /*if (ec_enc_tell(&st->enc, 0) < nbCompressedBytes*8 - 7)
642       celt_warning_int ("many unused bits: ", nbCompressedBytes*8-ec_enc_tell(&st->enc, 0));*/
643    /*printf ("%d\n", ec_enc_tell(&st->enc, 0)-8*nbCompressedBytes);*/
644    /* Finishing the stream with a 0101... pattern so that the decoder can check is everything's right */
645    {
646       int val = 0;
647       while (ec_enc_tell(&st->enc, 0) < nbCompressedBytes*8)
648       {
649          ec_enc_uint(&st->enc, val, 2);
650          val = 1-val;
651       }
652    }
653    ec_enc_done(&st->enc);
654    {
655       unsigned char *data;
656       int nbBytes = ec_byte_bytes(&st->buf);
657       if (nbBytes > nbCompressedBytes)
658       {
659          celt_warning_int ("got too many bytes:", nbBytes);
660          RESTORE_STACK;
661          return CELT_INTERNAL_ERROR;
662       }
663       /*printf ("%d\n", *nbBytes);*/
664       data = ec_byte_get_buffer(&st->buf);
665       for (i=0;i<nbBytes;i++)
666          compressed[i] = data[i];
667       for (;i<nbCompressedBytes;i++)
668          compressed[i] = 0;
669    }
670    /* Reset the packing for the next encoding */
671    ec_byte_reset(&st->buf);
672    ec_enc_init(&st->enc,&st->buf);
673
674    RESTORE_STACK;
675    return nbCompressedBytes;
676 }
677
678 #ifdef FIXED_POINT
679 #ifndef DISABLE_FLOAT_API
680 int celt_encode_float(CELTEncoder * restrict st, const float * pcm, float * optional_synthesis, unsigned char *compressed, int nbCompressedBytes)
681 {
682    int j, ret;
683    const int C = CHANNELS(st->mode);
684    const int N = st->block_size;
685    VARDECL(celt_int16_t, in);
686    SAVE_STACK;
687    ALLOC(in, C*N, celt_int16_t);
688
689    for (j=0;j<C*N;j++)
690      in[j] = FLOAT2INT16(pcm[j]);
691
692    if (optional_synthesis != NULL) {
693      ret=celt_encode(st,in,in,compressed,nbCompressedBytes);
694    /*Converts backwards for inplace operation*/
695       for (j=0;j=C*N;j++)
696          optional_synthesis[j]=in[j]*(1/32768.);
697    } else {
698      ret=celt_encode(st,in,NULL,compressed,nbCompressedBytes);
699    }
700    RESTORE_STACK;
701    return ret;
702
703 }
704 #endif /*DISABLE_FLOAT_API*/
705 #else
706 int celt_encode(CELTEncoder * restrict st, const celt_int16_t * pcm, celt_int16_t * optional_synthesis, unsigned char *compressed, int nbCompressedBytes)
707 {
708    int j, ret;
709    VARDECL(celt_sig_t, in);
710    const int C = CHANNELS(st->mode);
711    const int N = st->block_size;
712    SAVE_STACK;
713    ALLOC(in, C*N, celt_sig_t);
714    for (j=0;j<C*N;j++) {
715      in[j] = SCALEOUT(pcm[j]);
716    }
717
718    if (optional_synthesis != NULL) {
719       ret = celt_encode_float(st,in,in,compressed,nbCompressedBytes);
720       for (j=0;j<C*N;j++)
721          optional_synthesis[j] = FLOAT2INT16(in[j]);
722    } else {
723       ret = celt_encode_float(st,in,NULL,compressed,nbCompressedBytes);
724    }
725    RESTORE_STACK;
726    return ret;
727 }
728 #endif
729
730 int celt_encoder_ctl(CELTEncoder * restrict st, int request, ...)
731 {
732    va_list ap;
733    va_start(ap, request);
734    switch (request)
735    {
736       case CELT_SET_COMPLEXITY_REQUEST:
737       {
738          int value = va_arg(ap, int);
739          if (value<0 || value>10)
740             goto bad_arg;
741          if (value<=2)
742             st->pitch_enabled = 0;
743          else
744             st->pitch_enabled = 1;
745       }
746       break;
747       default:
748          goto bad_request;
749    }
750    va_end(ap);
751    return CELT_OK;
752 bad_arg:
753    va_end(ap);
754    return CELT_BAD_ARG;
755 bad_request:
756    va_end(ap);
757    return CELT_UNIMPLEMENTED;
758 }
759
760 /****************************************************************************/
761 /*                                                                          */
762 /*                                DECODER                                   */
763 /*                                                                          */
764 /****************************************************************************/
765
766
767 /** Decoder state 
768  @brief Decoder state
769  */
770 struct CELTDecoder {
771    const CELTMode *mode;
772    int frame_size;
773    int block_size;
774    int overlap;
775
776    ec_byte_buffer buf;
777    ec_enc         enc;
778
779    celt_sig_t * restrict preemph_memD;
780
781    celt_sig_t *out_mem;
782
783    celt_word16_t *oldBandE;
784    
785    int last_pitch_index;
786 };
787
788 CELTDecoder *celt_decoder_create(const CELTMode *mode)
789 {
790    int N, C;
791    CELTDecoder *st;
792
793    if (check_mode(mode) != CELT_OK)
794       return NULL;
795
796    N = mode->mdctSize;
797    C = CHANNELS(mode);
798    st = celt_alloc(sizeof(CELTDecoder));
799    
800    st->mode = mode;
801    st->frame_size = N;
802    st->block_size = N;
803    st->overlap = mode->overlap;
804
805    st->out_mem = celt_alloc((MAX_PERIOD+st->overlap)*C*sizeof(celt_sig_t));
806    
807    st->oldBandE = (celt_word16_t*)celt_alloc(C*mode->nbEBands*sizeof(celt_word16_t));
808
809    st->preemph_memD = (celt_sig_t*)celt_alloc(C*sizeof(celt_sig_t));
810
811    st->last_pitch_index = 0;
812    return st;
813 }
814
815 void celt_decoder_destroy(CELTDecoder *st)
816 {
817    if (st == NULL)
818    {
819       celt_warning("NULL passed to celt_encoder_destroy");
820       return;
821    }
822    if (check_mode(st->mode) != CELT_OK)
823       return;
824
825
826    celt_free(st->out_mem);
827    
828    celt_free(st->oldBandE);
829    
830    celt_free(st->preemph_memD);
831
832    celt_free(st);
833 }
834
835 /** Handles lost packets by just copying past data with the same offset as the last
836     pitch period */
837 static void celt_decode_lost(CELTDecoder * restrict st, celt_word16_t * restrict pcm)
838 {
839    int c, N;
840    int pitch_index;
841    int i, len;
842    VARDECL(celt_sig_t, freq);
843    const int C = CHANNELS(st->mode);
844    int offset;
845    SAVE_STACK;
846    N = st->block_size;
847    ALLOC(freq,C*N, celt_sig_t);         /**< Interleaved signal MDCTs */
848    
849    len = N+st->mode->overlap;
850 #if 0
851    pitch_index = st->last_pitch_index;
852    
853    /* Use the pitch MDCT as the "guessed" signal */
854    compute_mdcts(st->mode, st->mode->window, st->out_mem+pitch_index*C, freq);
855
856 #else
857    find_spectral_pitch(st->mode, st->mode->fft, &st->mode->psy, st->out_mem+MAX_PERIOD-len, st->out_mem, st->mode->window, NULL, len, MAX_PERIOD-len-100, &pitch_index);
858    pitch_index = MAX_PERIOD-len-pitch_index;
859    offset = MAX_PERIOD-pitch_index;
860    while (offset+len >= MAX_PERIOD)
861       offset -= pitch_index;
862    compute_mdcts(st->mode, 0, st->out_mem+offset*C, freq);
863    for (i=0;i<N;i++)
864       freq[i] = MULT16_32_Q15(QCONST16(.9f,15),freq[i]);
865 #endif
866    
867    
868    
869    CELT_MOVE(st->out_mem, st->out_mem+C*N, C*(MAX_PERIOD+st->mode->overlap-N));
870    /* Compute inverse MDCTs */
871    compute_inv_mdcts(st->mode, 0, freq, -1, 1, st->out_mem);
872
873    for (c=0;c<C;c++)
874    {
875       int j;
876       for (j=0;j<N;j++)
877       {
878          celt_sig_t tmp = MAC16_32_Q15(st->out_mem[C*(MAX_PERIOD-N)+C*j+c],
879                                 preemph,st->preemph_memD[c]);
880          st->preemph_memD[c] = tmp;
881          pcm[C*j+c] = SCALEOUT(SIG2WORD16(tmp));
882       }
883    }
884    RESTORE_STACK;
885 }
886
887 #ifdef FIXED_POINT
888 int celt_decode(CELTDecoder * restrict st, unsigned char *data, int len, celt_int16_t * restrict pcm)
889 {
890 #else
891 int celt_decode_float(CELTDecoder * restrict st, unsigned char *data, int len, celt_sig_t * restrict pcm)
892 {
893 #endif
894    int i, c, N, N4;
895    int has_pitch, has_fold;
896    int pitch_index;
897    int bits;
898    ec_dec dec;
899    ec_byte_buffer buf;
900    VARDECL(celt_sig_t, freq);
901    VARDECL(celt_norm_t, X);
902    VARDECL(celt_norm_t, P);
903    VARDECL(celt_ener_t, bandE);
904    VARDECL(celt_pgain_t, gains);
905    VARDECL(int, stereo_mode);
906    VARDECL(int, fine_quant);
907    VARDECL(int, pulses);
908    VARDECL(int, offsets);
909
910    int shortBlocks;
911    int transient_time;
912    int transient_shift;
913    const int C = CHANNELS(st->mode);
914    SAVE_STACK;
915
916    if (check_mode(st->mode) != CELT_OK)
917       return CELT_INVALID_MODE;
918
919    N = st->block_size;
920    N4 = (N-st->overlap)>>1;
921
922    ALLOC(freq, C*N, celt_sig_t); /**< Interleaved signal MDCTs */
923    ALLOC(X, C*N, celt_norm_t);         /**< Interleaved normalised MDCTs */
924    ALLOC(P, C*N, celt_norm_t);         /**< Interleaved normalised pitch MDCTs*/
925    ALLOC(bandE, st->mode->nbEBands*C, celt_ener_t);
926    ALLOC(gains, st->mode->nbPBands, celt_pgain_t);
927    
928    if (check_mode(st->mode) != CELT_OK)
929    {
930       RESTORE_STACK;
931       return CELT_INVALID_MODE;
932    }
933    if (data == NULL)
934    {
935       celt_decode_lost(st, pcm);
936       RESTORE_STACK;
937       return 0;
938    }
939    
940    ec_byte_readinit(&buf,data,len);
941    ec_dec_init(&dec,&buf);
942    
943    has_pitch = ec_dec_bits(&dec, 1);
944    if (has_pitch)
945    {
946       has_fold = ec_dec_bits(&dec, 1);
947       shortBlocks = 0;
948    } else if (st->mode->nbShortMdcts > 1){
949       shortBlocks = ec_dec_bits(&dec, 1);
950       has_fold = 1;
951    } else {
952       shortBlocks = 0;
953       has_fold = 1;
954    }
955    if (shortBlocks)
956    {
957       transient_shift = ec_dec_bits(&dec, 2);
958       if (transient_shift)
959          transient_time = ec_dec_uint(&dec, N+st->mode->overlap);
960       else
961          transient_time = 0;
962    } else {
963       transient_time = -1;
964       transient_shift = 0;
965    }
966    /* Get the pitch gains */
967    
968    /* Get the pitch index */
969    if (has_pitch)
970    {
971       has_pitch = unquant_pitch(gains, st->mode->nbPBands, &dec);
972       pitch_index = ec_dec_uint(&dec, MAX_PERIOD-(2*N-2*N4));
973       st->last_pitch_index = pitch_index;
974    } else {
975       /* FIXME: We could be more intelligent here and just not compute the MDCT */
976       pitch_index = 0;
977       for (i=0;i<st->mode->nbPBands;i++)
978          gains[i] = 0;
979    }
980
981    ALLOC(fine_quant, st->mode->nbEBands, int);
982    /* Get band energies */
983    unquant_coarse_energy(st->mode, bandE, st->oldBandE, len*8/3, st->mode->prob, &dec);
984    
985    ALLOC(pulses, st->mode->nbEBands, int);
986    ALLOC(offsets, st->mode->nbEBands, int);
987    ALLOC(stereo_mode, st->mode->nbEBands, int);
988    stereo_decision(st->mode, X, stereo_mode, st->mode->nbEBands);
989
990    for (i=0;i<st->mode->nbEBands;i++)
991       offsets[i] = 0;
992
993    bits = len*8 - ec_dec_tell(&dec, 0) - 1;
994    compute_allocation(st->mode, offsets, stereo_mode, bits, pulses, fine_quant);
995    /*bits = ec_dec_tell(&dec, 0);
996    compute_fine_allocation(st->mode, fine_quant, (20*C+len*8/5-(ec_dec_tell(&dec, 0)-bits))/C);*/
997    
998    unquant_fine_energy(st->mode, bandE, st->oldBandE, fine_quant, &dec);
999
1000
1001    if (has_pitch) 
1002    {
1003       VARDECL(celt_ener_t, bandEp);
1004       
1005       /* Pitch MDCT */
1006       compute_mdcts(st->mode, 0, st->out_mem+pitch_index*C, freq);
1007       ALLOC(bandEp, st->mode->nbEBands*C, celt_ener_t);
1008       compute_band_energies(st->mode, freq, bandEp);
1009       normalise_bands(st->mode, freq, P, bandEp);
1010    } else {
1011       for (i=0;i<C*N;i++)
1012          P[i] = 0;
1013    }
1014
1015    /* Apply pitch gains */
1016    pitch_quant_bands(st->mode, P, gains);
1017
1018    /* Decode fixed codebook and merge with pitch */
1019    unquant_bands(st->mode, X, P, bandE, stereo_mode, pulses, shortBlocks, len*8, &dec);
1020
1021    if (C==2)
1022    {
1023       renormalise_bands(st->mode, X);
1024    }
1025    /* Synthesis */
1026    denormalise_bands(st->mode, X, freq, bandE);
1027
1028
1029    CELT_MOVE(st->out_mem, st->out_mem+C*N, C*(MAX_PERIOD+st->overlap-N));
1030    /* Compute inverse MDCTs */
1031    compute_inv_mdcts(st->mode, shortBlocks, freq, transient_time, transient_shift, st->out_mem);
1032
1033    for (c=0;c<C;c++)
1034    {
1035       int j;
1036       for (j=0;j<N;j++)
1037       {
1038          celt_sig_t tmp = MAC16_32_Q15(st->out_mem[C*(MAX_PERIOD-N)+C*j+c],
1039                                 preemph,st->preemph_memD[c]);
1040          st->preemph_memD[c] = tmp;
1041          pcm[C*j+c] = SCALEOUT(SIG2WORD16(tmp));
1042       }
1043    }
1044
1045    {
1046       unsigned int val = 0;
1047       while (ec_dec_tell(&dec, 0) < len*8)
1048       {
1049          if (ec_dec_uint(&dec, 2) != val)
1050          {
1051             celt_warning("decode error");
1052             RESTORE_STACK;
1053             return CELT_CORRUPTED_DATA;
1054          }
1055          val = 1-val;
1056       }
1057    }
1058
1059    RESTORE_STACK;
1060    return 0;
1061    /*printf ("\n");*/
1062 }
1063
1064 #ifdef FIXED_POINT
1065 #ifndef DISABLE_FLOAT_API
1066 int celt_decode_float(CELTDecoder * restrict st, unsigned char *data, int len, float * restrict pcm)
1067 {
1068    int j, ret;
1069    const int C = CHANNELS(st->mode);
1070    const int N = st->block_size;
1071    VARDECL(celt_int16_t, out);
1072    SAVE_STACK;
1073    ALLOC(out, C*N, celt_int16_t);
1074
1075    ret=celt_decode(st, data, len, out);
1076
1077    for (j=0;j<C*N;j++)
1078      pcm[j]=out[j]*(1/32768.);
1079    RESTORE_STACK;
1080    return ret;
1081 }
1082 #endif /*DISABLE_FLOAT_API*/
1083 #else
1084 int celt_decode(CELTDecoder * restrict st, unsigned char *data, int len, celt_int16_t * restrict pcm)
1085 {
1086    int j, ret;
1087    VARDECL(celt_sig_t, out);
1088    const int C = CHANNELS(st->mode);
1089    const int N = st->block_size;
1090    SAVE_STACK;
1091    ALLOC(out, C*N, celt_sig_t);
1092
1093    ret=celt_decode_float(st, data, len, out);
1094
1095    for (j=0;j<C*N;j++)
1096      pcm[j] = FLOAT2INT16 (out[j]);
1097
1098    RESTORE_STACK;
1099    return ret;
1100 }
1101 #endif