Preparing for version 0.5.2
[opus.git] / libcelt / celt.c
1 /* (C) 2007-2008 Jean-Marc Valin, CSIRO
2    (C) 2008 Gregory Maxwell */
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_bands.h"
48 #include "psy.h"
49 #include "rate.h"
50 #include "stack_alloc.h"
51 #include "mathops.h"
52 #include "float_cast.h"
53 #include <stdarg.h>
54
55 static const celt_word16_t preemph = QCONST16(0.8f,15);
56
57 #ifdef FIXED_POINT
58 static const celt_word16_t transientWindow[16] = {
59      279,  1106,  2454,  4276,  6510,  9081, 11900, 14872,
60    17896, 20868, 23687, 26258, 28492, 30314, 31662, 32489};
61 #else
62 static const float transientWindow[16] = {
63    0.0085135, 0.0337639, 0.0748914, 0.1304955, 0.1986827, 0.2771308, 0.3631685, 0.4538658,
64    0.5461342, 0.6368315, 0.7228692, 0.8013173, 0.8695045, 0.9251086, 0.9662361, 0.9914865};
65 #endif
66
67    
68 /** Encoder state 
69  @brief Encoder state
70  */
71 struct CELTEncoder {
72    const CELTMode *mode;     /**< Mode used by the encoder */
73    int frame_size;
74    int block_size;
75    int overlap;
76    int channels;
77    
78    int pitch_enabled;
79    int pitch_available;
80
81    celt_word16_t * restrict preemph_memE; /* Input is 16-bit, so why bother with 32 */
82    celt_sig_t    * restrict preemph_memD;
83
84    celt_sig_t *in_mem;
85    celt_sig_t *out_mem;
86
87    celt_word16_t *oldBandE;
88 #ifdef EXP_PSY
89    celt_word16_t *psy_mem;
90    struct PsyDecay psy;
91 #endif
92 };
93
94 CELTEncoder *celt_encoder_create(const CELTMode *mode)
95 {
96    int N, C;
97    CELTEncoder *st;
98
99    if (check_mode(mode) != CELT_OK)
100       return NULL;
101
102    N = mode->mdctSize;
103    C = mode->nbChannels;
104    st = celt_alloc(sizeof(CELTEncoder));
105    
106    st->mode = mode;
107    st->frame_size = N;
108    st->block_size = N;
109    st->overlap = mode->overlap;
110
111    st->pitch_enabled = 1;
112    st->pitch_available = 1;
113
114    st->in_mem = celt_alloc(st->overlap*C*sizeof(celt_sig_t));
115    st->out_mem = celt_alloc((MAX_PERIOD+st->overlap)*C*sizeof(celt_sig_t));
116
117    st->oldBandE = (celt_word16_t*)celt_alloc(C*mode->nbEBands*sizeof(celt_word16_t));
118
119    st->preemph_memE = (celt_word16_t*)celt_alloc(C*sizeof(celt_word16_t));
120    st->preemph_memD = (celt_sig_t*)celt_alloc(C*sizeof(celt_sig_t));
121
122 #ifdef EXP_PSY
123    st->psy_mem = celt_alloc(MAX_PERIOD*sizeof(celt_word16_t));
124    psydecay_init(&st->psy, MAX_PERIOD/2, st->mode->Fs);
125 #endif
126
127    return st;
128 }
129
130 void celt_encoder_destroy(CELTEncoder *st)
131 {
132    if (st == NULL)
133    {
134       celt_warning("NULL passed to celt_encoder_destroy");
135       return;
136    }
137    if (check_mode(st->mode) != CELT_OK)
138       return;
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)float2int(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] = ABS32(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], ABS32(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->overlap;
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, N);
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, N2);
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, const celt_int16_t * pcm, celt_int16_t * optional_synthesis, unsigned char *compressed, int nbCompressedBytes)
369 {
370 #else
371 int celt_encode_float(CELTEncoder * restrict st, const celt_sig_t * pcm, celt_sig_t * optional_synthesis, 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    int has_fold=1;
379    ec_byte_buffer buf;
380    ec_enc         enc;
381    VARDECL(celt_sig_t, in);
382    VARDECL(celt_sig_t, freq);
383    VARDECL(celt_norm_t, X);
384    VARDECL(celt_norm_t, P);
385    VARDECL(celt_ener_t, bandE);
386    VARDECL(celt_pgain_t, gains);
387    VARDECL(int, stereo_mode);
388    VARDECL(int, fine_quant);
389    VARDECL(celt_word16_t, error);
390    VARDECL(int, pulses);
391    VARDECL(int, offsets);
392 #ifdef EXP_PSY
393    VARDECL(celt_word32_t, mask);
394    VARDECL(celt_word32_t, tonality);
395    VARDECL(celt_word32_t, bandM);
396    VARDECL(celt_ener_t, bandN);
397 #endif
398    int shortBlocks=0;
399    int transient_time;
400    int transient_shift;
401    const int C = CHANNELS(st->mode);
402    SAVE_STACK;
403
404    if (check_mode(st->mode) != CELT_OK)
405       return CELT_INVALID_MODE;
406
407    if (nbCompressedBytes<0)
408      return CELT_BAD_ARG; 
409
410    /* The memset is important for now in case the encoder doesn't fill up all the bytes */
411    CELT_MEMSET(compressed, 0, nbCompressedBytes);
412    ec_byte_writeinit_buffer(&buf, compressed, nbCompressedBytes);
413    ec_enc_init(&enc,&buf);
414
415    N = st->block_size;
416    N4 = (N-st->overlap)>>1;
417    ALLOC(in, 2*C*N-2*C*N4, celt_sig_t);
418
419    CELT_COPY(in, st->in_mem, C*st->overlap);
420    for (c=0;c<C;c++)
421    {
422       const celt_word16_t * restrict pcmp = pcm+c;
423       celt_sig_t * restrict inp = in+C*st->overlap+c;
424       for (i=0;i<N;i++)
425       {
426          /* Apply pre-emphasis */
427          celt_sig_t tmp = SCALEIN(SHL32(EXTEND32(*pcmp), SIG_SHIFT));
428          *inp = SUB32(tmp, SHR32(MULT16_16(preemph,st->preemph_memE[c]),3));
429          st->preemph_memE[c] = SCALEIN(*pcmp);
430          inp += C;
431          pcmp += C;
432       }
433    }
434    CELT_COPY(st->in_mem, in+C*(2*N-2*N4-st->overlap), C*st->overlap);
435    
436    /* Transient handling */
437    if (st->mode->nbShortMdcts > 1)
438    {
439       if (transient_analysis(in, N+st->overlap, C, &transient_time, &transient_shift))
440       {
441 #ifndef FIXED_POINT
442          float gain_1;
443 #endif
444          ec_enc_bits(&enc, 0, 1); /*Pitch off */
445          ec_enc_bits(&enc, 1, 1); /*Transient on */
446          ec_enc_bits(&enc, transient_shift, 2);
447          if (transient_shift)
448             ec_enc_uint(&enc, transient_time, N+st->overlap);
449          /* Apply the inverse shaping window */
450          if (transient_shift)
451          {
452 #ifdef FIXED_POINT
453             for (c=0;c<C;c++)
454                for (i=0;i<16;i++)
455                   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]);
456             for (c=0;c<C;c++)
457                for (i=transient_time;i<N+st->overlap;i++)
458                   in[C*i+c] = SHR32(in[C*i+c], transient_shift);
459 #else
460             for (c=0;c<C;c++)
461                for (i=0;i<16;i++)
462                   in[C*(transient_time+i-16)+c] /= 1+transientWindow[i]*((1<<transient_shift)-1);
463             gain_1 = 1./(1<<transient_shift);
464             for (c=0;c<C;c++)
465                for (i=transient_time;i<N+st->overlap;i++)
466                   in[C*i+c] *= gain_1;
467 #endif
468          }
469          shortBlocks = 1;
470       } else {
471          transient_time = -1;
472          transient_shift = 0;
473          shortBlocks = 0;
474       }
475    } else {
476       transient_time = -1;
477       transient_shift = 0;
478       shortBlocks = 0;
479    }
480
481    /* Pitch analysis: we do it early to save on the peak stack space */
482    /* Don't use pitch if there isn't enough data available yet, or if we're using shortBlocks */
483    has_pitch = st->pitch_enabled && (st->pitch_available >= MAX_PERIOD) && (!shortBlocks);
484 #ifdef EXP_PSY
485    ALLOC(tonality, MAX_PERIOD/4, celt_word16_t);
486    {
487       VARDECL(celt_word16_t, X);
488       ALLOC(X, MAX_PERIOD/2, celt_word16_t);
489       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);
490       compute_tonality(st->mode, X, st->psy_mem, MAX_PERIOD, tonality, MAX_PERIOD/4);
491    }
492 #else
493    if (has_pitch)
494    {
495       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);
496    }
497 #endif
498    ALLOC(freq, C*N, celt_sig_t); /**< Interleaved signal MDCTs */
499    
500    /* Compute MDCTs */
501    compute_mdcts(st->mode, shortBlocks, in, freq);
502
503 #ifdef EXP_PSY
504    ALLOC(mask, N, celt_sig_t);
505    compute_mdct_masking(&st->psy, freq, tonality, st->psy_mem, mask, C*N);
506    /*for (i=0;i<256;i++)
507       printf ("%f %f %f ", freq[i], tonality[i], mask[i]);
508    printf ("\n");*/
509 #endif
510
511    /* Deferred allocation after find_spectral_pitch() to reduce the peak memory usage */
512    ALLOC(X, C*N, celt_norm_t);         /**< Interleaved normalised MDCTs */
513    ALLOC(P, C*N, celt_norm_t);         /**< Interleaved normalised pitch MDCTs*/
514    ALLOC(bandE,st->mode->nbEBands*C, celt_ener_t);
515    ALLOC(gains,st->mode->nbPBands, celt_pgain_t);
516
517
518    /* Band normalisation */
519    compute_band_energies(st->mode, freq, bandE);
520    normalise_bands(st->mode, freq, X, bandE);
521
522 #ifdef EXP_PSY
523    ALLOC(bandN,C*st->mode->nbEBands, celt_ener_t);
524    ALLOC(bandM,st->mode->nbEBands, celt_ener_t);
525    compute_noise_energies(st->mode, freq, tonality, bandN);
526
527    /*for (i=0;i<st->mode->nbEBands;i++)
528       printf ("%f ", (.1+bandN[i])/(.1+bandE[i]));
529    printf ("\n");*/
530    has_fold = 0;
531    for (i=st->mode->nbPBands;i<st->mode->nbEBands;i++)
532       if (bandN[i] < .4*bandE[i])
533          has_fold++;
534    /*printf ("%d\n", has_fold);*/
535    if (has_fold>=2)
536       has_fold = 0;
537    else
538       has_fold = 1;
539    for (i=0;i<N;i++)
540       mask[i] = sqrt(mask[i]);
541    compute_band_energies(st->mode, mask, bandM);
542    /*for (i=0;i<st->mode->nbEBands;i++)
543       printf ("%f %f ", bandE[i], bandM[i]);
544    printf ("\n");*/
545 #endif
546
547    /* Compute MDCTs of the pitch part */
548    if (has_pitch)
549    {
550       celt_word32_t curr_power, pitch_power=0;
551       /* Normalise the pitch vector as well (discard the energies) */
552       VARDECL(celt_ener_t, bandEp);
553       
554       compute_mdcts(st->mode, 0, st->out_mem+pitch_index*C, freq);
555       ALLOC(bandEp, st->mode->nbEBands*st->mode->nbChannels, celt_ener_t);
556       compute_band_energies(st->mode, freq, bandEp);
557       normalise_bands(st->mode, freq, P, bandEp);
558       pitch_power = bandEp[0]+bandEp[1]+bandEp[2];
559       /* Check if we can safely use the pitch (i.e. effective gain isn't too high) */
560       curr_power = bandE[0]+bandE[1]+bandE[2];
561       if ((MULT16_32_Q15(QCONST16(.1f, 15),curr_power) + QCONST32(10.f,ENER_SHIFT) < pitch_power))
562       {
563          /* Pitch prediction */
564          has_pitch = compute_pitch_gain(st->mode, X, P, gains);
565       } else {
566          has_pitch = 0;
567       }
568    }
569    
570    if (has_pitch) 
571    {  
572       ec_enc_bits(&enc, has_pitch, 1); /* Pitch flag */
573       ec_enc_bits(&enc, has_fold, 1); /* Folding flag */
574       ec_enc_uint(&enc, pitch_index, MAX_PERIOD-(2*N-2*N4));
575    } else {
576       if (!shortBlocks)
577       {
578          ec_enc_bits(&enc, 0, 1); /* Pitch off */
579          if (st->mode->nbShortMdcts > 1)
580            ec_enc_bits(&enc, 0, 1); /* Transient off */
581       }
582       has_fold = 1;
583       /* No pitch, so we just pretend we found a gain of zero */
584       for (i=0;i<st->mode->nbPBands;i++)
585          gains[i] = 0;
586       for (i=0;i<C*N;i++)
587          P[i] = 0;
588    }
589
590 #ifdef STDIN_TUNING2
591    static int fine_quant[30];
592    static int pulses[30];
593    static int init=0;
594    if (!init)
595    {
596       for (i=0;i<st->mode->nbEBands;i++)
597          scanf("%d ", &fine_quant[i]);
598       for (i=0;i<st->mode->nbEBands;i++)
599          scanf("%d ", &pulses[i]);
600       init = 1;
601    }
602 #else
603    ALLOC(fine_quant, st->mode->nbEBands, int);
604    ALLOC(pulses, st->mode->nbEBands, int);
605 #endif
606
607    /* Bit allocation */
608    ALLOC(error, C*st->mode->nbEBands, celt_word16_t);
609    quant_coarse_energy(st->mode, bandE, st->oldBandE, nbCompressedBytes*8/3, st->mode->prob, error, &enc);
610    
611    ALLOC(offsets, st->mode->nbEBands, int);
612    ALLOC(stereo_mode, st->mode->nbEBands, int);
613    stereo_decision(st->mode, X, stereo_mode, st->mode->nbEBands);
614
615    for (i=0;i<st->mode->nbEBands;i++)
616       offsets[i] = 0;
617    bits = nbCompressedBytes*8 - ec_enc_tell(&enc, 0) - 1;
618    if (has_pitch)
619       bits -= st->mode->nbPBands;
620 #ifndef STDIN_TUNING
621    compute_allocation(st->mode, offsets, stereo_mode, bits, pulses, fine_quant);
622 #endif
623
624    quant_fine_energy(st->mode, bandE, st->oldBandE, error, fine_quant, &enc);
625
626    /* Residual quantisation */
627    quant_bands(st->mode, X, P, NULL, has_pitch, gains, bandE, stereo_mode, pulses, shortBlocks, has_fold, nbCompressedBytes*8, &enc);
628
629    /* Re-synthesis of the coded audio if required */
630    if (st->pitch_available>0 || optional_synthesis!=NULL)
631    {
632       if (st->pitch_available>0 && st->pitch_available<MAX_PERIOD)
633         st->pitch_available+=st->frame_size;
634
635       if (C==2)
636          renormalise_bands(st->mode, X);
637       /* Synthesis */
638       denormalise_bands(st->mode, X, freq, bandE);
639       
640       
641       CELT_MOVE(st->out_mem, st->out_mem+C*N, C*(MAX_PERIOD+st->overlap-N));
642       
643       compute_inv_mdcts(st->mode, shortBlocks, freq, transient_time, transient_shift, st->out_mem);
644       /* De-emphasis and put everything back at the right place in the synthesis history */
645       if (optional_synthesis != NULL) {
646          for (c=0;c<C;c++)
647          {
648             int j;
649             for (j=0;j<N;j++)
650             {
651                celt_sig_t tmp = MAC16_32_Q15(st->out_mem[C*(MAX_PERIOD-N)+C*j+c],
652                                    preemph,st->preemph_memD[c]);
653                st->preemph_memD[c] = tmp;
654                optional_synthesis[C*j+c] = SCALEOUT(SIG2WORD16(tmp));
655             }
656          }
657       }
658    }
659    /*fprintf (stderr, "remaining bits after encode = %d\n", nbCompressedBytes*8-ec_enc_tell(&st->enc, 0));*/
660    /*if (ec_enc_tell(&st->enc, 0) < nbCompressedBytes*8 - 7)
661       celt_warning_int ("many unused bits: ", nbCompressedBytes*8-ec_enc_tell(&st->enc, 0));*/
662    /*printf ("%d\n", ec_enc_tell(&st->enc, 0)-8*nbCompressedBytes);*/
663    /* Finishing the stream with a 0101... pattern so that the decoder can check is everything's right */
664    {
665       int val = 0;
666       while (ec_enc_tell(&enc, 0) < nbCompressedBytes*8)
667       {
668          ec_enc_uint(&enc, val, 2);
669          val = 1-val;
670       }
671    }
672    ec_enc_done(&enc);
673    {
674       /*unsigned char *data;*/
675       int nbBytes = ec_byte_bytes(&buf);
676       if (nbBytes > nbCompressedBytes)
677       {
678          celt_warning_int ("got too many bytes:", nbBytes);
679          RESTORE_STACK;
680          return CELT_INTERNAL_ERROR;
681       }
682    }
683
684    RESTORE_STACK;
685    return nbCompressedBytes;
686 }
687
688 #ifdef FIXED_POINT
689 #ifndef DISABLE_FLOAT_API
690 int celt_encode_float(CELTEncoder * restrict st, const float * pcm, float * optional_synthesis, unsigned char *compressed, int nbCompressedBytes)
691 {
692    int j, ret;
693    const int C = CHANNELS(st->mode);
694    const int N = st->block_size;
695    VARDECL(celt_int16_t, in);
696    SAVE_STACK;
697    ALLOC(in, C*N, celt_int16_t);
698
699    for (j=0;j<C*N;j++)
700      in[j] = FLOAT2INT16(pcm[j]);
701
702    if (optional_synthesis != NULL) {
703      ret=celt_encode(st,in,in,compressed,nbCompressedBytes);
704       for (j=0;j<C*N;j++)
705          optional_synthesis[j]=in[j]*(1/32768.);
706    } else {
707      ret=celt_encode(st,in,NULL,compressed,nbCompressedBytes);
708    }
709    RESTORE_STACK;
710    return ret;
711
712 }
713 #endif /*DISABLE_FLOAT_API*/
714 #else
715 int celt_encode(CELTEncoder * restrict st, const celt_int16_t * pcm, celt_int16_t * optional_synthesis, unsigned char *compressed, int nbCompressedBytes)
716 {
717    int j, ret;
718    VARDECL(celt_sig_t, in);
719    const int C = CHANNELS(st->mode);
720    const int N = st->block_size;
721    SAVE_STACK;
722    ALLOC(in, C*N, celt_sig_t);
723    for (j=0;j<C*N;j++) {
724      in[j] = SCALEOUT(pcm[j]);
725    }
726
727    if (optional_synthesis != NULL) {
728       ret = celt_encode_float(st,in,in,compressed,nbCompressedBytes);
729       for (j=0;j<C*N;j++)
730          optional_synthesis[j] = FLOAT2INT16(in[j]);
731    } else {
732       ret = celt_encode_float(st,in,NULL,compressed,nbCompressedBytes);
733    }
734    RESTORE_STACK;
735    return ret;
736 }
737 #endif
738
739 int celt_encoder_ctl(CELTEncoder * restrict st, int request, ...)
740 {
741    va_list ap;
742    va_start(ap, request);
743    switch (request)
744    {
745       case CELT_SET_COMPLEXITY_REQUEST:
746       {
747          int value = va_arg(ap, int);
748          if (value<0 || value>10)
749             goto bad_arg;
750          if (value<=2) {
751             st->pitch_enabled = 0; 
752             st->pitch_available = 0;
753          } else {
754               st->pitch_enabled = 1;
755               if (st->pitch_available<1)
756                 st->pitch_available = 1;
757          }   
758       }
759       break;
760       case CELT_SET_LTP_REQUEST:
761       {
762          int value = va_arg(ap, int);
763          if (value<0 || value>1 || (value==1 && st->pitch_available==0))
764             goto bad_arg;
765          if (value==0)
766             st->pitch_enabled = 0;
767          else
768             st->pitch_enabled = 1;
769       }
770       break;
771       default:
772          goto bad_request;
773    }
774    va_end(ap);
775    return CELT_OK;
776 bad_arg:
777    va_end(ap);
778    return CELT_BAD_ARG;
779 bad_request:
780    va_end(ap);
781    return CELT_UNIMPLEMENTED;
782 }
783
784 /****************************************************************************/
785 /*                                                                          */
786 /*                                DECODER                                   */
787 /*                                                                          */
788 /****************************************************************************/
789
790
791 /** Decoder state 
792  @brief Decoder state
793  */
794 struct CELTDecoder {
795    const CELTMode *mode;
796    int frame_size;
797    int block_size;
798    int overlap;
799
800    ec_byte_buffer buf;
801    ec_enc         enc;
802
803    celt_sig_t * restrict preemph_memD;
804
805    celt_sig_t *out_mem;
806
807    celt_word16_t *oldBandE;
808    
809    int last_pitch_index;
810 };
811
812 CELTDecoder *celt_decoder_create(const CELTMode *mode)
813 {
814    int N, C;
815    CELTDecoder *st;
816
817    if (check_mode(mode) != CELT_OK)
818       return NULL;
819
820    N = mode->mdctSize;
821    C = CHANNELS(mode);
822    st = celt_alloc(sizeof(CELTDecoder));
823    
824    st->mode = mode;
825    st->frame_size = N;
826    st->block_size = N;
827    st->overlap = mode->overlap;
828
829    st->out_mem = celt_alloc((MAX_PERIOD+st->overlap)*C*sizeof(celt_sig_t));
830    
831    st->oldBandE = (celt_word16_t*)celt_alloc(C*mode->nbEBands*sizeof(celt_word16_t));
832
833    st->preemph_memD = (celt_sig_t*)celt_alloc(C*sizeof(celt_sig_t));
834
835    st->last_pitch_index = 0;
836    return st;
837 }
838
839 void celt_decoder_destroy(CELTDecoder *st)
840 {
841    if (st == NULL)
842    {
843       celt_warning("NULL passed to celt_encoder_destroy");
844       return;
845    }
846    if (check_mode(st->mode) != CELT_OK)
847       return;
848
849
850    celt_free(st->out_mem);
851    
852    celt_free(st->oldBandE);
853    
854    celt_free(st->preemph_memD);
855
856    celt_free(st);
857 }
858
859 /** Handles lost packets by just copying past data with the same offset as the last
860     pitch period */
861 static void celt_decode_lost(CELTDecoder * restrict st, celt_word16_t * restrict pcm)
862 {
863    int c, N;
864    int pitch_index;
865    int i, len;
866    VARDECL(celt_sig_t, freq);
867    const int C = CHANNELS(st->mode);
868    int offset;
869    SAVE_STACK;
870    N = st->block_size;
871    ALLOC(freq,C*N, celt_sig_t);         /**< Interleaved signal MDCTs */
872    
873    len = N+st->mode->overlap;
874 #if 0
875    pitch_index = st->last_pitch_index;
876    
877    /* Use the pitch MDCT as the "guessed" signal */
878    compute_mdcts(st->mode, st->mode->window, st->out_mem+pitch_index*C, freq);
879
880 #else
881    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);
882    pitch_index = MAX_PERIOD-len-pitch_index;
883    offset = MAX_PERIOD-pitch_index;
884    while (offset+len >= MAX_PERIOD)
885       offset -= pitch_index;
886    compute_mdcts(st->mode, 0, st->out_mem+offset*C, freq);
887    for (i=0;i<N;i++)
888       freq[i] = ADD32(EPSILON, MULT16_32_Q15(QCONST16(.9f,15),freq[i]));
889 #endif
890    
891    
892    
893    CELT_MOVE(st->out_mem, st->out_mem+C*N, C*(MAX_PERIOD+st->mode->overlap-N));
894    /* Compute inverse MDCTs */
895    compute_inv_mdcts(st->mode, 0, freq, -1, 1, st->out_mem);
896
897    for (c=0;c<C;c++)
898    {
899       int j;
900       for (j=0;j<N;j++)
901       {
902          celt_sig_t tmp = MAC16_32_Q15(st->out_mem[C*(MAX_PERIOD-N)+C*j+c],
903                                 preemph,st->preemph_memD[c]);
904          st->preemph_memD[c] = tmp;
905          pcm[C*j+c] = SCALEOUT(SIG2WORD16(tmp));
906       }
907    }
908    RESTORE_STACK;
909 }
910
911 #ifdef FIXED_POINT
912 int celt_decode(CELTDecoder * restrict st, unsigned char *data, int len, celt_int16_t * restrict pcm)
913 {
914 #else
915 int celt_decode_float(CELTDecoder * restrict st, unsigned char *data, int len, celt_sig_t * restrict pcm)
916 {
917 #endif
918    int i, c, N, N4;
919    int has_pitch, has_fold;
920    int pitch_index;
921    int bits;
922    ec_dec dec;
923    ec_byte_buffer buf;
924    VARDECL(celt_sig_t, freq);
925    VARDECL(celt_norm_t, X);
926    VARDECL(celt_norm_t, P);
927    VARDECL(celt_ener_t, bandE);
928    VARDECL(celt_pgain_t, gains);
929    VARDECL(int, stereo_mode);
930    VARDECL(int, fine_quant);
931    VARDECL(int, pulses);
932    VARDECL(int, offsets);
933
934    int shortBlocks;
935    int transient_time;
936    int transient_shift;
937    const int C = CHANNELS(st->mode);
938    SAVE_STACK;
939
940    if (check_mode(st->mode) != CELT_OK)
941       return CELT_INVALID_MODE;
942
943    N = st->block_size;
944    N4 = (N-st->overlap)>>1;
945
946    ALLOC(freq, C*N, celt_sig_t); /**< Interleaved signal MDCTs */
947    ALLOC(X, C*N, celt_norm_t);         /**< Interleaved normalised MDCTs */
948    ALLOC(P, C*N, celt_norm_t);         /**< Interleaved normalised pitch MDCTs*/
949    ALLOC(bandE, st->mode->nbEBands*C, celt_ener_t);
950    ALLOC(gains, st->mode->nbPBands, celt_pgain_t);
951    
952    if (check_mode(st->mode) != CELT_OK)
953    {
954       RESTORE_STACK;
955       return CELT_INVALID_MODE;
956    }
957    if (data == NULL)
958    {
959       celt_decode_lost(st, pcm);
960       RESTORE_STACK;
961       return 0;
962    }
963    if (len<0) {
964      RESTORE_STACK;
965      return CELT_BAD_ARG;
966    }
967    
968    ec_byte_readinit(&buf,data,len);
969    ec_dec_init(&dec,&buf);
970    
971    has_pitch = ec_dec_bits(&dec, 1);
972    if (has_pitch)
973    {
974       has_fold = ec_dec_bits(&dec, 1);
975       shortBlocks = 0;
976    } else if (st->mode->nbShortMdcts > 1){
977       shortBlocks = ec_dec_bits(&dec, 1);
978       has_fold = 1;
979    } else {
980       shortBlocks = 0;
981       has_fold = 1;
982    }
983    if (shortBlocks)
984    {
985       transient_shift = ec_dec_bits(&dec, 2);
986       if (transient_shift)
987          transient_time = ec_dec_uint(&dec, N+st->mode->overlap);
988       else
989          transient_time = 0;
990    } else {
991       transient_time = -1;
992       transient_shift = 0;
993    }
994    
995    if (has_pitch)
996    {
997       pitch_index = ec_dec_uint(&dec, MAX_PERIOD-(2*N-2*N4));
998       st->last_pitch_index = pitch_index;
999    } else {
1000       pitch_index = 0;
1001       for (i=0;i<st->mode->nbPBands;i++)
1002          gains[i] = 0;
1003    }
1004
1005    ALLOC(fine_quant, st->mode->nbEBands, int);
1006    /* Get band energies */
1007    unquant_coarse_energy(st->mode, bandE, st->oldBandE, len*8/3, st->mode->prob, &dec);
1008    
1009    ALLOC(pulses, st->mode->nbEBands, int);
1010    ALLOC(offsets, st->mode->nbEBands, int);
1011    ALLOC(stereo_mode, st->mode->nbEBands, int);
1012    stereo_decision(st->mode, X, stereo_mode, st->mode->nbEBands);
1013
1014    for (i=0;i<st->mode->nbEBands;i++)
1015       offsets[i] = 0;
1016
1017    bits = len*8 - ec_dec_tell(&dec, 0) - 1;
1018    if (has_pitch)
1019       bits -= st->mode->nbPBands;
1020    compute_allocation(st->mode, offsets, stereo_mode, bits, pulses, fine_quant);
1021    /*bits = ec_dec_tell(&dec, 0);
1022    compute_fine_allocation(st->mode, fine_quant, (20*C+len*8/5-(ec_dec_tell(&dec, 0)-bits))/C);*/
1023    
1024    unquant_fine_energy(st->mode, bandE, st->oldBandE, fine_quant, &dec);
1025
1026
1027    if (has_pitch) 
1028    {
1029       VARDECL(celt_ener_t, bandEp);
1030       
1031       /* Pitch MDCT */
1032       compute_mdcts(st->mode, 0, st->out_mem+pitch_index*C, freq);
1033       ALLOC(bandEp, st->mode->nbEBands*C, celt_ener_t);
1034       compute_band_energies(st->mode, freq, bandEp);
1035       normalise_bands(st->mode, freq, P, bandEp);
1036       /* Apply pitch gains */
1037    } else {
1038       for (i=0;i<C*N;i++)
1039          P[i] = 0;
1040    }
1041
1042    /* Decode fixed codebook and merge with pitch */
1043    unquant_bands(st->mode, X, P, has_pitch, gains, bandE, stereo_mode, pulses, shortBlocks, has_fold, len*8, &dec);
1044
1045    if (C==2)
1046    {
1047       renormalise_bands(st->mode, X);
1048    }
1049    /* Synthesis */
1050    denormalise_bands(st->mode, X, freq, bandE);
1051
1052
1053    CELT_MOVE(st->out_mem, st->out_mem+C*N, C*(MAX_PERIOD+st->overlap-N));
1054    /* Compute inverse MDCTs */
1055    compute_inv_mdcts(st->mode, shortBlocks, freq, transient_time, transient_shift, st->out_mem);
1056
1057    for (c=0;c<C;c++)
1058    {
1059       int j;
1060       for (j=0;j<N;j++)
1061       {
1062          celt_sig_t tmp = MAC16_32_Q15(st->out_mem[C*(MAX_PERIOD-N)+C*j+c],
1063                                 preemph,st->preemph_memD[c]);
1064          st->preemph_memD[c] = tmp;
1065          pcm[C*j+c] = SCALEOUT(SIG2WORD16(tmp));
1066       }
1067    }
1068
1069    {
1070       unsigned int val = 0;
1071       while (ec_dec_tell(&dec, 0) < len*8)
1072       {
1073          if (ec_dec_uint(&dec, 2) != val)
1074          {
1075             celt_warning("decode error");
1076             RESTORE_STACK;
1077             return CELT_CORRUPTED_DATA;
1078          }
1079          val = 1-val;
1080       }
1081    }
1082
1083    RESTORE_STACK;
1084    return 0;
1085    /*printf ("\n");*/
1086 }
1087
1088 #ifdef FIXED_POINT
1089 #ifndef DISABLE_FLOAT_API
1090 int celt_decode_float(CELTDecoder * restrict st, unsigned char *data, int len, float * restrict pcm)
1091 {
1092    int j, ret;
1093    const int C = CHANNELS(st->mode);
1094    const int N = st->block_size;
1095    VARDECL(celt_int16_t, out);
1096    SAVE_STACK;
1097    ALLOC(out, C*N, celt_int16_t);
1098
1099    ret=celt_decode(st, data, len, out);
1100
1101    for (j=0;j<C*N;j++)
1102      pcm[j]=out[j]*(1/32768.);
1103    RESTORE_STACK;
1104    return ret;
1105 }
1106 #endif /*DISABLE_FLOAT_API*/
1107 #else
1108 int celt_decode(CELTDecoder * restrict st, unsigned char *data, int len, celt_int16_t * restrict pcm)
1109 {
1110    int j, ret;
1111    VARDECL(celt_sig_t, out);
1112    const int C = CHANNELS(st->mode);
1113    const int N = st->block_size;
1114    SAVE_STACK;
1115    ALLOC(out, C*N, celt_sig_t);
1116
1117    ret=celt_decode_float(st, data, len, out);
1118
1119    for (j=0;j<C*N;j++)
1120      pcm[j] = FLOAT2INT16 (out[j]);
1121
1122    RESTORE_STACK;
1123    return ret;
1124 }
1125 #endif