175922c5704699d9ddc80dcea6fa127f237774fc
[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 SIG2INT16(celt_sig_t x)
157 {
158    x = PSHR32(x, SIG_SHIFT);
159    x = MAX32(x, -32768);
160    x = MIN32(x, 32767);
161 #ifdef FIXED_POINT
162    return EXTRACT16(x);
163 #else
164    return (celt_int16_t)floor(.5+x);
165 #endif
166 }
167
168 static int transient_analysis(celt_word32_t *in, int len, int C, int *transient_time, int *transient_shift)
169 {
170    int c, i, n;
171    celt_word32_t ratio;
172    /* FIXME: Remove the floats here */
173    VARDECL(celt_word32_t, begin);
174    SAVE_STACK;
175    ALLOC(begin, len, celt_word32_t);
176    for (i=0;i<len;i++)
177       begin[i] = EXTEND32(ABS16(SHR32(in[C*i],SIG_SHIFT)));
178    for (c=1;c<C;c++)
179    {
180       for (i=0;i<len;i++)
181          begin[i] = MAX32(begin[i], EXTEND32(ABS16(SHR32(in[C*i+c],SIG_SHIFT))));
182    }
183    for (i=1;i<len;i++)
184       begin[i] = MAX32(begin[i-1],begin[i]);
185    n = -1;
186    for (i=8;i<len-8;i++)
187    {
188       if (begin[i] < MULT16_32_Q15(QCONST16(.2f,15),begin[len-1]))
189          n=i;
190    }
191    if (n<32)
192    {
193       n = -1;
194       ratio = 0;
195    } else {
196       ratio = DIV32(begin[len-1],1+begin[n-16]);
197    }
198    /*printf ("%d %f\n", n, ratio*ratio);*/
199    if (ratio < 0)
200       ratio = 0;
201    if (ratio > 1000)
202       ratio = 1000;
203    ratio *= ratio;
204    if (ratio < 50)
205       *transient_shift = 0;
206    else if (ratio < 256)
207       *transient_shift = 1;
208    else if (ratio < 4096)
209       *transient_shift = 2;
210    else
211       *transient_shift = 3;
212    *transient_time = n;
213    
214    RESTORE_STACK;
215    return ratio > 20;
216 }
217
218 /** Apply window and compute the MDCT for all sub-frames and all channels in a frame */
219 static void compute_mdcts(const CELTMode *mode, int shortBlocks, celt_sig_t * restrict in, celt_sig_t * restrict out)
220 {
221    const int C = CHANNELS(mode);
222    if (C==1 && !shortBlocks)
223    {
224       const mdct_lookup *lookup = MDCT(mode);
225       const int overlap = OVERLAP(mode);
226       mdct_forward(lookup, in, out, mode->window, overlap);
227    } else if (!shortBlocks) {
228       const mdct_lookup *lookup = MDCT(mode);
229       const int overlap = OVERLAP(mode);
230       const int N = FRAMESIZE(mode);
231       int c;
232       VARDECL(celt_word32_t, x);
233       VARDECL(celt_word32_t, tmp);
234       SAVE_STACK;
235       ALLOC(x, N+overlap, celt_word32_t);
236       ALLOC(tmp, N, celt_word32_t);
237       for (c=0;c<C;c++)
238       {
239          int j;
240          for (j=0;j<N+overlap;j++)
241             x[j] = in[C*j+c];
242          mdct_forward(lookup, x, tmp, mode->window, overlap);
243          /* Interleaving the sub-frames */
244          for (j=0;j<N;j++)
245             out[C*j+c] = tmp[j];
246       }
247       RESTORE_STACK;
248    } else {
249       const mdct_lookup *lookup = &mode->shortMdct;
250       const int overlap = mode->shortMdctSize;
251       const int N = mode->shortMdctSize;
252       int b, c;
253       VARDECL(celt_word32_t, x);
254       VARDECL(celt_word32_t, tmp);
255       SAVE_STACK;
256       ALLOC(x, N+overlap, celt_word32_t);
257       ALLOC(tmp, N, celt_word32_t);
258       for (c=0;c<C;c++)
259       {
260          int B = mode->nbShortMdcts;
261          for (b=0;b<B;b++)
262          {
263             int j;
264             for (j=0;j<N+overlap;j++)
265                x[j] = in[C*(b*N+j)+c];
266             mdct_forward(lookup, x, tmp, mode->window, overlap);
267             /* Interleaving the sub-frames */
268             for (j=0;j<N;j++)
269                out[C*(j*B+b)+c] = tmp[j];
270          }
271       }
272       RESTORE_STACK;
273    }
274 }
275
276 /** Compute the IMDCT and apply window for all sub-frames and all channels in a frame */
277 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)
278 {
279    int c, N4;
280    const int C = CHANNELS(mode);
281    const int N = FRAMESIZE(mode);
282    const int overlap = OVERLAP(mode);
283    N4 = (N-overlap)>>1;
284    for (c=0;c<C;c++)
285    {
286       int j;
287       if (transient_shift==0 && C==1 && !shortBlocks) {
288          const mdct_lookup *lookup = MDCT(mode);
289          mdct_backward(lookup, X, out_mem+C*(MAX_PERIOD-N-N4), mode->window, overlap);
290       } else if (!shortBlocks) {
291          const mdct_lookup *lookup = MDCT(mode);
292          VARDECL(celt_word32_t, x);
293          VARDECL(celt_word32_t, tmp);
294          SAVE_STACK;
295          ALLOC(x, 2*N, celt_word32_t);
296          ALLOC(tmp, N, celt_word32_t);
297          /* De-interleaving the sub-frames */
298          for (j=0;j<N;j++)
299             tmp[j] = X[C*j+c];
300          /* Prevents problems from the imdct doing the overlap-add */
301          CELT_MEMSET(x+N4, 0, overlap);
302          mdct_backward(lookup, tmp, x, mode->window, overlap);
303          celt_assert(transient_shift == 0);
304          /* The first and last part would need to be set to zero if we actually
305             wanted to use them. */
306          for (j=0;j<overlap;j++)
307             out_mem[C*(MAX_PERIOD-N)+C*j+c] += x[j+N4];
308          for (j=0;j<overlap;j++)
309             out_mem[C*(MAX_PERIOD)+C*(overlap-j-1)+c] = x[2*N-j-N4-1];
310          for (j=0;j<2*N4;j++)
311             out_mem[C*(MAX_PERIOD-N)+C*(j+overlap)+c] = x[j+N4+overlap];
312          RESTORE_STACK;
313       } else {
314          int b;
315          const int N2 = mode->shortMdctSize;
316          const int B = mode->nbShortMdcts;
317          const mdct_lookup *lookup = &mode->shortMdct;
318          VARDECL(celt_word32_t, x);
319          VARDECL(celt_word32_t, tmp);
320          SAVE_STACK;
321          ALLOC(x, 2*N, celt_word32_t);
322          ALLOC(tmp, N, celt_word32_t);
323          /* Prevents problems from the imdct doing the overlap-add */
324          CELT_MEMSET(x+N4, 0, overlap);
325          for (b=0;b<B;b++)
326          {
327             /* De-interleaving the sub-frames */
328             for (j=0;j<N2;j++)
329                tmp[j] = X[C*(j*B+b)+c];
330             mdct_backward(lookup, tmp, x+N4+N2*b, mode->window, overlap);
331          }
332          if (transient_shift > 0)
333          {
334 #ifdef FIXED_POINT
335             for (j=0;j<16;j++)
336                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));
337             for (j=transient_time;j<N+overlap;j++)
338                x[N4+j] = SHL32(x[N4+j], transient_shift);
339 #else
340             for (j=0;j<16;j++)
341                x[N4+transient_time+j-16] *= 1+transientWindow[j]*((1<<transient_shift)-1);
342             for (j=transient_time;j<N+overlap;j++)
343                x[N4+j] *= 1<<transient_shift;
344 #endif
345          }
346          /* The first and last part would need to be set to zero if we actually
347          wanted to use them. */
348          for (j=0;j<overlap;j++)
349             out_mem[C*(MAX_PERIOD-N)+C*j+c] += x[j+N4];
350          for (j=0;j<overlap;j++)
351             out_mem[C*(MAX_PERIOD)+C*(overlap-j-1)+c] = x[2*N-j-N4-1];
352          for (j=0;j<2*N4;j++)
353             out_mem[C*(MAX_PERIOD-N)+C*(j+overlap)+c] = x[j+N4+overlap];
354          RESTORE_STACK;
355       }
356    }
357 }
358
359 int celt_encode(CELTEncoder * restrict st, celt_int16_t * restrict pcm, unsigned char *compressed, int nbCompressedBytes)
360 {
361    int i, c, N, N4;
362    int has_pitch;
363    int pitch_index;
364    celt_word32_t curr_power, pitch_power;
365    VARDECL(celt_sig_t, in);
366    VARDECL(celt_sig_t, freq);
367    VARDECL(celt_norm_t, X);
368    VARDECL(celt_norm_t, P);
369    VARDECL(celt_ener_t, bandE);
370    VARDECL(celt_pgain_t, gains);
371    VARDECL(int, stereo_mode);
372 #ifdef EXP_PSY
373    VARDECL(celt_word32_t, mask);
374 #endif
375    int shortBlocks=0;
376    int transient_time;
377    int transient_shift;
378    const int C = CHANNELS(st->mode);
379    SAVE_STACK;
380
381    if (check_mode(st->mode) != CELT_OK)
382       return CELT_INVALID_MODE;
383
384    N = st->block_size;
385    N4 = (N-st->overlap)>>1;
386    ALLOC(in, 2*C*N-2*C*N4, celt_sig_t);
387
388    CELT_COPY(in, st->in_mem, C*st->overlap);
389    for (c=0;c<C;c++)
390    {
391       const celt_int16_t * restrict pcmp = pcm+c;
392       celt_sig_t * restrict inp = in+C*st->overlap+c;
393       for (i=0;i<N;i++)
394       {
395          /* Apply pre-emphasis */
396          celt_sig_t tmp = SHL32(EXTEND32(*pcmp), SIG_SHIFT);
397          *inp = SUB32(tmp, SHR32(MULT16_16(preemph,st->preemph_memE[c]),1));
398          st->preemph_memE[c] = *pcmp;
399          inp += C;
400          pcmp += C;
401       }
402    }
403    CELT_COPY(st->in_mem, in+C*(2*N-2*N4-st->overlap), C*st->overlap);
404    
405    if (st->mode->nbShortMdcts > 1)
406    {
407       if (transient_analysis(in, N+st->overlap, C, &transient_time, &transient_shift))
408       {
409 #ifndef FIXED_POINT
410          float gain_1;
411 #endif
412          ec_enc_bits(&st->enc, 1, 1);
413          ec_enc_bits(&st->enc, transient_shift, 2);
414          if (transient_shift)
415             ec_enc_uint(&st->enc, transient_time, N+st->overlap);
416          if (transient_shift)
417          {
418 #ifdef FIXED_POINT
419             for (c=0;c<C;c++)
420                for (i=0;i<16;i++)
421                   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]);
422             for (c=0;c<C;c++)
423                for (i=transient_time;i<N+st->overlap;i++)
424                   in[C*i+c] = SHR32(in[C*i+c], transient_shift);
425 #else
426             for (c=0;c<C;c++)
427                for (i=0;i<16;i++)
428                   in[C*(transient_time+i-16)+c] /= 1+transientWindow[i]*((1<<transient_shift)-1);
429             gain_1 = 1./(1<<transient_shift);
430             for (c=0;c<C;c++)
431                for (i=transient_time;i<N+st->overlap;i++)
432                   in[C*i+c] *= gain_1;
433 #endif
434          }
435          shortBlocks = 1;
436       } else {
437          ec_enc_bits(&st->enc, 0, 1);
438          transient_time = -1;
439          transient_shift = 0;
440          shortBlocks = 0;
441       }
442    } else {
443       transient_time = -1;
444       transient_shift = 0;
445       shortBlocks = 0;
446    }
447    /* Pitch analysis: we do it early to save on the peak stack space */
448    if (!shortBlocks)
449       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);
450
451    ALLOC(freq, C*N, celt_sig_t); /**< Interleaved signal MDCTs */
452    
453    /*for (i=0;i<(B+1)*C*N;i++) printf ("%f(%d) ", in[i], i); printf ("\n");*/
454    /* Compute MDCTs */
455    compute_mdcts(st->mode, shortBlocks, in, freq);
456
457 #ifdef EXP_PSY
458    CELT_MOVE(st->psy_mem, st->out_mem+N, MAX_PERIOD+st->overlap-N);
459    for (i=0;i<N;i++)
460       st->psy_mem[MAX_PERIOD+st->overlap-N+i] = in[C*(st->overlap+i)];
461    for (c=1;c<C;c++)
462       for (i=0;i<N;i++)
463          st->psy_mem[MAX_PERIOD+st->overlap-N+i] += in[C*(st->overlap+i)+c];
464
465    ALLOC(mask, N, celt_sig_t);
466    compute_mdct_masking(&st->psy, freq, st->psy_mem, mask, C*N);
467
468    /* Invert and stretch the mask to length of X 
469       For some reason, I get better results by using the sqrt instead,
470       although there's no valid reason to. Must investigate further */
471    for (i=0;i<C*N;i++)
472       mask[i] = 1/(.1+mask[i]);
473 #endif
474    
475    /* Deferred allocation after find_spectral_pitch() to reduce the peak memory usage */
476    ALLOC(X, C*N, celt_norm_t);         /**< Interleaved normalised MDCTs */
477    ALLOC(P, C*N, celt_norm_t);         /**< Interleaved normalised pitch MDCTs*/
478    ALLOC(bandE,st->mode->nbEBands*C, celt_ener_t);
479    ALLOC(gains,st->mode->nbPBands, celt_pgain_t);
480
481    /*printf ("%f %f\n", curr_power, pitch_power);*/
482    /*int j;
483    for (j=0;j<B*N;j++)
484       printf ("%f ", X[j]);
485    for (j=0;j<B*N;j++)
486       printf ("%f ", P[j]);
487    printf ("\n");*/
488
489    /* Band normalisation */
490    compute_band_energies(st->mode, freq, bandE);
491    normalise_bands(st->mode, freq, X, bandE);
492    /*for (i=0;i<st->mode->nbEBands;i++)printf("%f ", bandE[i]);printf("\n");*/
493    /*for (i=0;i<N*B*C;i++)printf("%f ", X[i]);printf("\n");*/
494
495    /* Compute MDCTs of the pitch part */
496    if (!shortBlocks)
497       compute_mdcts(st->mode, 0, st->out_mem+pitch_index*C, freq);
498
499    {
500       /* Normalise the pitch vector as well (discard the energies) */
501       VARDECL(celt_ener_t, bandEp);
502       ALLOC(bandEp, st->mode->nbEBands*st->mode->nbChannels, celt_ener_t);
503       compute_band_energies(st->mode, freq, bandEp);
504       normalise_bands(st->mode, freq, P, bandEp);
505       pitch_power = bandEp[0]+bandEp[1]+bandEp[2];
506    }
507    curr_power = bandE[0]+bandE[1]+bandE[2];
508    /* Check if we can safely use the pitch (i.e. effective gain isn't too high) */
509    if (!shortBlocks && (MULT16_32_Q15(QCONST16(.1f, 15),curr_power) + QCONST32(10.f,ENER_SHIFT) < pitch_power))
510    {
511       /* Simulates intensity stereo */
512       /*for (i=30;i<N*B;i++)
513          X[i*C+1] = P[i*C+1] = 0;*/
514
515       /* Pitch prediction */
516       compute_pitch_gain(st->mode, X, P, gains);
517       has_pitch = quant_pitch(gains, st->mode->nbPBands, &st->enc);
518       if (has_pitch)
519          ec_enc_uint(&st->enc, pitch_index, MAX_PERIOD-(2*N-2*N4));
520    } else {
521       /* No pitch, so we just pretend we found a gain of zero */
522       for (i=0;i<st->mode->nbPBands;i++)
523          gains[i] = 0;
524       ec_enc_bits(&st->enc, 0, 7);
525       for (i=0;i<C*N;i++)
526          P[i] = 0;
527    }
528    quant_energy(st->mode, bandE, st->oldBandE, 20*C+nbCompressedBytes*8/5, st->mode->prob, &st->enc);
529
530    ALLOC(stereo_mode, st->mode->nbEBands, int);
531    stereo_decision(st->mode, X, stereo_mode, st->mode->nbEBands);
532
533    pitch_quant_bands(st->mode, P, gains);
534
535    /*for (i=0;i<B*N;i++) printf("%f ",P[i]);printf("\n");*/
536
537    /* Residual quantisation */
538    quant_bands(st->mode, X, P, NULL, bandE, stereo_mode, nbCompressedBytes*8, shortBlocks, &st->enc);
539    
540    if (C==2)
541    {
542       renormalise_bands(st->mode, X);
543    }
544    /* Synthesis */
545    denormalise_bands(st->mode, X, freq, bandE);
546
547
548    CELT_MOVE(st->out_mem, st->out_mem+C*N, C*(MAX_PERIOD+st->overlap-N));
549
550    compute_inv_mdcts(st->mode, shortBlocks, freq, transient_time, transient_shift, st->out_mem);
551    /* De-emphasis and put everything back at the right place in the synthesis history */
552 #ifndef SHORTCUTS
553    for (c=0;c<C;c++)
554    {
555       int j;
556       celt_sig_t * restrict outp=st->out_mem+C*(MAX_PERIOD-N)+c;
557       celt_int16_t * restrict pcmp = pcm+c;
558       for (j=0;j<N;j++)
559       {
560          celt_sig_t tmp = ADD32(*outp, MULT16_32_Q15(preemph,st->preemph_memD[c]));
561          st->preemph_memD[c] = tmp;
562          *pcmp = SIG2INT16(tmp);
563          pcmp += C;
564          outp += C;
565       }
566    }
567 #endif
568    /*if (ec_enc_tell(&st->enc, 0) < nbCompressedBytes*8 - 7)
569       celt_warning_int ("many unused bits: ", nbCompressedBytes*8-ec_enc_tell(&st->enc, 0));*/
570    /*printf ("%d\n", ec_enc_tell(&st->enc, 0)-8*nbCompressedBytes);*/
571    /* Finishing the stream with a 0101... pattern so that the decoder can check is everything's right */
572    {
573       int val = 0;
574       while (ec_enc_tell(&st->enc, 0) < nbCompressedBytes*8)
575       {
576          ec_enc_uint(&st->enc, val, 2);
577          val = 1-val;
578       }
579    }
580    ec_enc_done(&st->enc);
581    {
582       unsigned char *data;
583       int nbBytes = ec_byte_bytes(&st->buf);
584       if (nbBytes > nbCompressedBytes)
585       {
586          celt_warning_int ("got too many bytes:", nbBytes);
587          RESTORE_STACK;
588          return CELT_INTERNAL_ERROR;
589       }
590       /*printf ("%d\n", *nbBytes);*/
591       data = ec_byte_get_buffer(&st->buf);
592       for (i=0;i<nbBytes;i++)
593          compressed[i] = data[i];
594       for (;i<nbCompressedBytes;i++)
595          compressed[i] = 0;
596    }
597    /* Reset the packing for the next encoding */
598    ec_byte_reset(&st->buf);
599    ec_enc_init(&st->enc,&st->buf);
600
601    RESTORE_STACK;
602    return nbCompressedBytes;
603 }
604
605
606 /****************************************************************************/
607 /*                                                                          */
608 /*                                DECODER                                   */
609 /*                                                                          */
610 /****************************************************************************/
611
612
613 /** Decoder state 
614  @brief Decoder state
615  */
616 struct CELTDecoder {
617    const CELTMode *mode;
618    int frame_size;
619    int block_size;
620    int overlap;
621
622    ec_byte_buffer buf;
623    ec_enc         enc;
624
625    celt_sig_t * restrict preemph_memD;
626
627    celt_sig_t *out_mem;
628
629    celt_word16_t *oldBandE;
630    
631    int last_pitch_index;
632 };
633
634 CELTDecoder *celt_decoder_create(const CELTMode *mode)
635 {
636    int N, C;
637    CELTDecoder *st;
638
639    if (check_mode(mode) != CELT_OK)
640       return NULL;
641
642    N = mode->mdctSize;
643    C = CHANNELS(mode);
644    st = celt_alloc(sizeof(CELTDecoder));
645    
646    st->mode = mode;
647    st->frame_size = N;
648    st->block_size = N;
649    st->overlap = mode->overlap;
650
651    st->out_mem = celt_alloc((MAX_PERIOD+st->overlap)*C*sizeof(celt_sig_t));
652    
653    st->oldBandE = (celt_word16_t*)celt_alloc(C*mode->nbEBands*sizeof(celt_word16_t));
654
655    st->preemph_memD = (celt_sig_t*)celt_alloc(C*sizeof(celt_sig_t));;
656
657    st->last_pitch_index = 0;
658    return st;
659 }
660
661 void celt_decoder_destroy(CELTDecoder *st)
662 {
663    if (st == NULL)
664    {
665       celt_warning("NULL passed to celt_encoder_destroy");
666       return;
667    }
668    if (check_mode(st->mode) != CELT_OK)
669       return;
670
671
672    celt_free(st->out_mem);
673    
674    celt_free(st->oldBandE);
675    
676    celt_free(st->preemph_memD);
677
678    celt_free(st);
679 }
680
681 /** Handles lost packets by just copying past data with the same offset as the last
682     pitch period */
683 static void celt_decode_lost(CELTDecoder * restrict st, short * restrict pcm)
684 {
685    int c, N;
686    int pitch_index;
687    int i, len;
688    VARDECL(celt_sig_t, freq);
689    const int C = CHANNELS(st->mode);
690    int offset;
691    SAVE_STACK;
692    N = st->block_size;
693    ALLOC(freq,C*N, celt_sig_t);         /**< Interleaved signal MDCTs */
694    
695    len = N+st->mode->overlap;
696 #if 0
697    pitch_index = st->last_pitch_index;
698    
699    /* Use the pitch MDCT as the "guessed" signal */
700    compute_mdcts(st->mode, st->mode->window, st->out_mem+pitch_index*C, freq);
701
702 #else
703    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);
704    pitch_index = MAX_PERIOD-len-pitch_index;
705    offset = MAX_PERIOD-pitch_index;
706    while (offset+len >= MAX_PERIOD)
707       offset -= pitch_index;
708    compute_mdcts(st->mode, 0, st->out_mem+offset*C, freq);
709    for (i=0;i<N;i++)
710       freq[i] = MULT16_32_Q15(QCONST16(.9f,15),freq[i]);
711 #endif
712    
713    
714    
715    CELT_MOVE(st->out_mem, st->out_mem+C*N, C*(MAX_PERIOD+st->mode->overlap-N));
716    /* Compute inverse MDCTs */
717    compute_inv_mdcts(st->mode, 0, freq, -1, 1, st->out_mem);
718
719    for (c=0;c<C;c++)
720    {
721       int j;
722       for (j=0;j<N;j++)
723       {
724          celt_sig_t tmp = ADD32(st->out_mem[C*(MAX_PERIOD-N)+C*j+c],
725                                 MULT16_32_Q15(preemph,st->preemph_memD[c]));
726          st->preemph_memD[c] = tmp;
727          pcm[C*j+c] = SIG2INT16(tmp);
728       }
729    }
730    RESTORE_STACK;
731 }
732
733 int celt_decode(CELTDecoder * restrict st, unsigned char *data, int len, celt_int16_t * restrict pcm)
734 {
735    int c, N, N4;
736    int has_pitch;
737    int pitch_index;
738    ec_dec dec;
739    ec_byte_buffer buf;
740    VARDECL(celt_sig_t, freq);
741    VARDECL(celt_norm_t, X);
742    VARDECL(celt_norm_t, P);
743    VARDECL(celt_ener_t, bandE);
744    VARDECL(celt_pgain_t, gains);
745    VARDECL(int, stereo_mode);
746    int shortBlocks;
747    int transient_time;
748    int transient_shift;
749    const int C = CHANNELS(st->mode);
750    SAVE_STACK;
751
752    if (check_mode(st->mode) != CELT_OK)
753       return CELT_INVALID_MODE;
754
755    N = st->block_size;
756    N4 = (N-st->overlap)>>1;
757
758    ALLOC(freq, C*N, celt_sig_t); /**< Interleaved signal MDCTs */
759    ALLOC(X, C*N, celt_norm_t);         /**< Interleaved normalised MDCTs */
760    ALLOC(P, C*N, celt_norm_t);         /**< Interleaved normalised pitch MDCTs*/
761    ALLOC(bandE, st->mode->nbEBands*C, celt_ener_t);
762    ALLOC(gains, st->mode->nbPBands, celt_pgain_t);
763    
764    if (check_mode(st->mode) != CELT_OK)
765    {
766       RESTORE_STACK;
767       return CELT_INVALID_MODE;
768    }
769    if (data == NULL)
770    {
771       celt_decode_lost(st, pcm);
772       RESTORE_STACK;
773       return 0;
774    }
775    
776    ec_byte_readinit(&buf,data,len);
777    ec_dec_init(&dec,&buf);
778    
779    if (st->mode->nbShortMdcts > 1)
780    {
781       shortBlocks = ec_dec_bits(&dec, 1);
782       if (shortBlocks)
783       {
784          transient_shift = ec_dec_bits(&dec, 2);
785          if (transient_shift)
786             transient_time = ec_dec_uint(&dec, N+st->mode->overlap);
787          else
788             transient_time = 0;
789       } else {
790          transient_time = -1;
791          transient_shift = 0;
792       }
793    } else {
794       shortBlocks = 0;
795       transient_time = -1;
796       transient_shift = 0;
797    }
798    /* Get the pitch gains */
799    has_pitch = unquant_pitch(gains, st->mode->nbPBands, &dec);
800    
801    /* Get the pitch index */
802    if (has_pitch)
803    {
804       pitch_index = ec_dec_uint(&dec, MAX_PERIOD-(2*N-2*N4));
805       st->last_pitch_index = pitch_index;
806    } else {
807       /* FIXME: We could be more intelligent here and just not compute the MDCT */
808       pitch_index = 0;
809    }
810
811    /* Get band energies */
812    unquant_energy(st->mode, bandE, st->oldBandE, 20*C+len*8/5, st->mode->prob, &dec);
813
814    /* Pitch MDCT */
815    compute_mdcts(st->mode, 0, st->out_mem+pitch_index*C, freq);
816
817    {
818       VARDECL(celt_ener_t, bandEp);
819       ALLOC(bandEp, st->mode->nbEBands*C, celt_ener_t);
820       compute_band_energies(st->mode, freq, bandEp);
821       normalise_bands(st->mode, freq, P, bandEp);
822    }
823
824    ALLOC(stereo_mode, st->mode->nbEBands, int);
825    stereo_decision(st->mode, X, stereo_mode, st->mode->nbEBands);
826    /* Apply pitch gains */
827    pitch_quant_bands(st->mode, P, gains);
828
829    /* Decode fixed codebook and merge with pitch */
830    unquant_bands(st->mode, X, P, bandE, stereo_mode, len*8, shortBlocks, &dec);
831
832    if (C==2)
833    {
834       renormalise_bands(st->mode, X);
835    }
836    /* Synthesis */
837    denormalise_bands(st->mode, X, freq, bandE);
838
839
840    CELT_MOVE(st->out_mem, st->out_mem+C*N, C*(MAX_PERIOD+st->overlap-N));
841    /* Compute inverse MDCTs */
842    compute_inv_mdcts(st->mode, shortBlocks, freq, transient_time, transient_shift, st->out_mem);
843
844    for (c=0;c<C;c++)
845    {
846       int j;
847       const celt_sig_t * restrict outp=st->out_mem+C*(MAX_PERIOD-N)+c;
848       celt_int16_t * restrict pcmp = pcm+c;
849       for (j=0;j<N;j++)
850       {
851          celt_sig_t tmp = ADD32(*outp, MULT16_32_Q15(preemph,st->preemph_memD[c]));
852          st->preemph_memD[c] = tmp;
853          *pcmp = SIG2INT16(tmp);
854          pcmp += C;
855          outp += C;
856       }
857    }
858
859    {
860       unsigned int val = 0;
861       while (ec_dec_tell(&dec, 0) < len*8)
862       {
863          if (ec_dec_uint(&dec, 2) != val)
864          {
865             celt_warning("decode error");
866             RESTORE_STACK;
867             return CELT_CORRUPTED_DATA;
868          }
869          val = 1-val;
870       }
871    }
872
873    RESTORE_STACK;
874    return 0;
875    /*printf ("\n");*/
876 }
877