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