2774edd69c2054acf2ac386a86ca4b495432e1e9
[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
660    /*fprintf (stderr, "remaining bits after encode = %d\n", nbCompressedBytes*8-ec_enc_tell(&enc, 0));*/
661    /*if (ec_enc_tell(&enc, 0) < nbCompressedBytes*8 - 7)
662       celt_warning_int ("many unused bits: ", nbCompressedBytes*8-ec_enc_tell(&enc, 0));*/
663
664    /* Finishing the stream with a 0101... pattern so that the decoder can check is everything's right */
665    {
666       int val = 0;
667       while (ec_enc_tell(&enc, 0) < nbCompressedBytes*8)
668       {
669          ec_enc_uint(&enc, val, 2);
670          val = 1-val;
671       }
672    }
673    ec_enc_done(&enc);
674    {
675       /*unsigned char *data;*/
676       int nbBytes = ec_byte_bytes(&buf);
677       if (nbBytes > nbCompressedBytes)
678       {
679          celt_warning_int ("got too many bytes:", nbBytes);
680          RESTORE_STACK;
681          return CELT_INTERNAL_ERROR;
682       }
683    }
684
685    RESTORE_STACK;
686    return nbCompressedBytes;
687 }
688
689 #ifdef FIXED_POINT
690 #ifndef DISABLE_FLOAT_API
691 int celt_encode_float(CELTEncoder * restrict st, const float * pcm, float * optional_synthesis, unsigned char *compressed, int nbCompressedBytes)
692 {
693    int j, ret;
694    const int C = CHANNELS(st->mode);
695    const int N = st->block_size;
696    VARDECL(celt_int16_t, in);
697    SAVE_STACK;
698    ALLOC(in, C*N, celt_int16_t);
699
700    for (j=0;j<C*N;j++)
701      in[j] = FLOAT2INT16(pcm[j]);
702
703    if (optional_synthesis != NULL) {
704      ret=celt_encode(st,in,in,compressed,nbCompressedBytes);
705       for (j=0;j<C*N;j++)
706          optional_synthesis[j]=in[j]*(1/32768.);
707    } else {
708      ret=celt_encode(st,in,NULL,compressed,nbCompressedBytes);
709    }
710    RESTORE_STACK;
711    return ret;
712
713 }
714 #endif /*DISABLE_FLOAT_API*/
715 #else
716 int celt_encode(CELTEncoder * restrict st, const celt_int16_t * pcm, celt_int16_t * optional_synthesis, unsigned char *compressed, int nbCompressedBytes)
717 {
718    int j, ret;
719    VARDECL(celt_sig_t, in);
720    const int C = CHANNELS(st->mode);
721    const int N = st->block_size;
722    SAVE_STACK;
723    ALLOC(in, C*N, celt_sig_t);
724    for (j=0;j<C*N;j++) {
725      in[j] = SCALEOUT(pcm[j]);
726    }
727
728    if (optional_synthesis != NULL) {
729       ret = celt_encode_float(st,in,in,compressed,nbCompressedBytes);
730       for (j=0;j<C*N;j++)
731          optional_synthesis[j] = FLOAT2INT16(in[j]);
732    } else {
733       ret = celt_encode_float(st,in,NULL,compressed,nbCompressedBytes);
734    }
735    RESTORE_STACK;
736    return ret;
737 }
738 #endif
739
740 int celt_encoder_ctl(CELTEncoder * restrict st, int request, ...)
741 {
742    va_list ap;
743    va_start(ap, request);
744    switch (request)
745    {
746       case CELT_SET_COMPLEXITY_REQUEST:
747       {
748          int value = va_arg(ap, int);
749          if (value<0 || value>10)
750             goto bad_arg;
751          if (value<=2) {
752             st->pitch_enabled = 0; 
753             st->pitch_available = 0;
754          } else {
755               st->pitch_enabled = 1;
756               if (st->pitch_available<1)
757                 st->pitch_available = 1;
758          }   
759       }
760       break;
761       case CELT_SET_LTP_REQUEST:
762       {
763          int value = va_arg(ap, int);
764          if (value<0 || value>1 || (value==1 && st->pitch_available==0))
765             goto bad_arg;
766          if (value==0)
767             st->pitch_enabled = 0;
768          else
769             st->pitch_enabled = 1;
770       }
771       break;
772       default:
773          goto bad_request;
774    }
775    va_end(ap);
776    return CELT_OK;
777 bad_arg:
778    va_end(ap);
779    return CELT_BAD_ARG;
780 bad_request:
781    va_end(ap);
782    return CELT_UNIMPLEMENTED;
783 }
784
785 /****************************************************************************/
786 /*                                                                          */
787 /*                                DECODER                                   */
788 /*                                                                          */
789 /****************************************************************************/
790
791
792 /** Decoder state 
793  @brief Decoder state
794  */
795 struct CELTDecoder {
796    const CELTMode *mode;
797    int frame_size;
798    int block_size;
799    int overlap;
800
801    ec_byte_buffer buf;
802    ec_enc         enc;
803
804    celt_sig_t * restrict preemph_memD;
805
806    celt_sig_t *out_mem;
807
808    celt_word16_t *oldBandE;
809    
810    int last_pitch_index;
811 };
812
813 CELTDecoder *celt_decoder_create(const CELTMode *mode)
814 {
815    int N, C;
816    CELTDecoder *st;
817
818    if (check_mode(mode) != CELT_OK)
819       return NULL;
820
821    N = mode->mdctSize;
822    C = CHANNELS(mode);
823    st = celt_alloc(sizeof(CELTDecoder));
824    
825    st->mode = mode;
826    st->frame_size = N;
827    st->block_size = N;
828    st->overlap = mode->overlap;
829
830    st->out_mem = celt_alloc((MAX_PERIOD+st->overlap)*C*sizeof(celt_sig_t));
831    
832    st->oldBandE = (celt_word16_t*)celt_alloc(C*mode->nbEBands*sizeof(celt_word16_t));
833
834    st->preemph_memD = (celt_sig_t*)celt_alloc(C*sizeof(celt_sig_t));
835
836    st->last_pitch_index = 0;
837    return st;
838 }
839
840 void celt_decoder_destroy(CELTDecoder *st)
841 {
842    if (st == NULL)
843    {
844       celt_warning("NULL passed to celt_encoder_destroy");
845       return;
846    }
847    if (check_mode(st->mode) != CELT_OK)
848       return;
849
850
851    celt_free(st->out_mem);
852    
853    celt_free(st->oldBandE);
854    
855    celt_free(st->preemph_memD);
856
857    celt_free(st);
858 }
859
860 /** Handles lost packets by just copying past data with the same offset as the last
861     pitch period */
862 static void celt_decode_lost(CELTDecoder * restrict st, celt_word16_t * restrict pcm)
863 {
864    int c, N;
865    int pitch_index;
866    int i, len;
867    VARDECL(celt_sig_t, freq);
868    const int C = CHANNELS(st->mode);
869    int offset;
870    SAVE_STACK;
871    N = st->block_size;
872    ALLOC(freq,C*N, celt_sig_t);         /**< Interleaved signal MDCTs */
873    
874    len = N+st->mode->overlap;
875 #if 0
876    pitch_index = st->last_pitch_index;
877    
878    /* Use the pitch MDCT as the "guessed" signal */
879    compute_mdcts(st->mode, st->mode->window, st->out_mem+pitch_index*C, freq);
880
881 #else
882    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);
883    pitch_index = MAX_PERIOD-len-pitch_index;
884    offset = MAX_PERIOD-pitch_index;
885    while (offset+len >= MAX_PERIOD)
886       offset -= pitch_index;
887    compute_mdcts(st->mode, 0, st->out_mem+offset*C, freq);
888    for (i=0;i<N;i++)
889       freq[i] = ADD32(EPSILON, MULT16_32_Q15(QCONST16(.9f,15),freq[i]));
890 #endif
891    
892    
893    
894    CELT_MOVE(st->out_mem, st->out_mem+C*N, C*(MAX_PERIOD+st->mode->overlap-N));
895    /* Compute inverse MDCTs */
896    compute_inv_mdcts(st->mode, 0, freq, -1, 1, st->out_mem);
897
898    for (c=0;c<C;c++)
899    {
900       int j;
901       for (j=0;j<N;j++)
902       {
903          celt_sig_t tmp = MAC16_32_Q15(st->out_mem[C*(MAX_PERIOD-N)+C*j+c],
904                                 preemph,st->preemph_memD[c]);
905          st->preemph_memD[c] = tmp;
906          pcm[C*j+c] = SCALEOUT(SIG2WORD16(tmp));
907       }
908    }
909    RESTORE_STACK;
910 }
911
912 #ifdef FIXED_POINT
913 int celt_decode(CELTDecoder * restrict st, unsigned char *data, int len, celt_int16_t * restrict pcm)
914 {
915 #else
916 int celt_decode_float(CELTDecoder * restrict st, unsigned char *data, int len, celt_sig_t * restrict pcm)
917 {
918 #endif
919    int i, c, N, N4;
920    int has_pitch, has_fold;
921    int pitch_index;
922    int bits;
923    ec_dec dec;
924    ec_byte_buffer buf;
925    VARDECL(celt_sig_t, freq);
926    VARDECL(celt_norm_t, X);
927    VARDECL(celt_norm_t, P);
928    VARDECL(celt_ener_t, bandE);
929    VARDECL(celt_pgain_t, gains);
930    VARDECL(int, stereo_mode);
931    VARDECL(int, fine_quant);
932    VARDECL(int, pulses);
933    VARDECL(int, offsets);
934
935    int shortBlocks;
936    int transient_time;
937    int transient_shift;
938    const int C = CHANNELS(st->mode);
939    SAVE_STACK;
940
941    if (check_mode(st->mode) != CELT_OK)
942       return CELT_INVALID_MODE;
943
944    N = st->block_size;
945    N4 = (N-st->overlap)>>1;
946
947    ALLOC(freq, C*N, celt_sig_t); /**< Interleaved signal MDCTs */
948    ALLOC(X, C*N, celt_norm_t);         /**< Interleaved normalised MDCTs */
949    ALLOC(P, C*N, celt_norm_t);         /**< Interleaved normalised pitch MDCTs*/
950    ALLOC(bandE, st->mode->nbEBands*C, celt_ener_t);
951    ALLOC(gains, st->mode->nbPBands, celt_pgain_t);
952    
953    if (check_mode(st->mode) != CELT_OK)
954    {
955       RESTORE_STACK;
956       return CELT_INVALID_MODE;
957    }
958    if (data == NULL)
959    {
960       celt_decode_lost(st, pcm);
961       RESTORE_STACK;
962       return 0;
963    }
964    if (len<0) {
965      RESTORE_STACK;
966      return CELT_BAD_ARG;
967    }
968    
969    ec_byte_readinit(&buf,data,len);
970    ec_dec_init(&dec,&buf);
971    
972    has_pitch = ec_dec_bits(&dec, 1);
973    if (has_pitch)
974    {
975       has_fold = ec_dec_bits(&dec, 1);
976       shortBlocks = 0;
977    } else if (st->mode->nbShortMdcts > 1){
978       shortBlocks = ec_dec_bits(&dec, 1);
979       has_fold = 1;
980    } else {
981       shortBlocks = 0;
982       has_fold = 1;
983    }
984    if (shortBlocks)
985    {
986       transient_shift = ec_dec_bits(&dec, 2);
987       if (transient_shift)
988          transient_time = ec_dec_uint(&dec, N+st->mode->overlap);
989       else
990          transient_time = 0;
991    } else {
992       transient_time = -1;
993       transient_shift = 0;
994    }
995    
996    if (has_pitch)
997    {
998       pitch_index = ec_dec_uint(&dec, MAX_PERIOD-(2*N-2*N4));
999       st->last_pitch_index = pitch_index;
1000    } else {
1001       pitch_index = 0;
1002       for (i=0;i<st->mode->nbPBands;i++)
1003          gains[i] = 0;
1004    }
1005
1006    ALLOC(fine_quant, st->mode->nbEBands, int);
1007    /* Get band energies */
1008    unquant_coarse_energy(st->mode, bandE, st->oldBandE, len*8/3, st->mode->prob, &dec);
1009    
1010    ALLOC(pulses, st->mode->nbEBands, int);
1011    ALLOC(offsets, st->mode->nbEBands, int);
1012    ALLOC(stereo_mode, st->mode->nbEBands, int);
1013    stereo_decision(st->mode, X, stereo_mode, st->mode->nbEBands);
1014
1015    for (i=0;i<st->mode->nbEBands;i++)
1016       offsets[i] = 0;
1017
1018    bits = len*8 - ec_dec_tell(&dec, 0) - 1;
1019    if (has_pitch)
1020       bits -= st->mode->nbPBands;
1021    compute_allocation(st->mode, offsets, stereo_mode, bits, pulses, fine_quant);
1022    /*bits = ec_dec_tell(&dec, 0);
1023    compute_fine_allocation(st->mode, fine_quant, (20*C+len*8/5-(ec_dec_tell(&dec, 0)-bits))/C);*/
1024    
1025    unquant_fine_energy(st->mode, bandE, st->oldBandE, fine_quant, &dec);
1026
1027
1028    if (has_pitch) 
1029    {
1030       VARDECL(celt_ener_t, bandEp);
1031       
1032       /* Pitch MDCT */
1033       compute_mdcts(st->mode, 0, st->out_mem+pitch_index*C, freq);
1034       ALLOC(bandEp, st->mode->nbEBands*C, celt_ener_t);
1035       compute_band_energies(st->mode, freq, bandEp);
1036       normalise_bands(st->mode, freq, P, bandEp);
1037       /* Apply pitch gains */
1038    } else {
1039       for (i=0;i<C*N;i++)
1040          P[i] = 0;
1041    }
1042
1043    /* Decode fixed codebook and merge with pitch */
1044    unquant_bands(st->mode, X, P, has_pitch, gains, bandE, stereo_mode, pulses, shortBlocks, has_fold, len*8, &dec);
1045
1046    if (C==2)
1047    {
1048       renormalise_bands(st->mode, X);
1049    }
1050    /* Synthesis */
1051    denormalise_bands(st->mode, X, freq, bandE);
1052
1053
1054    CELT_MOVE(st->out_mem, st->out_mem+C*N, C*(MAX_PERIOD+st->overlap-N));
1055    /* Compute inverse MDCTs */
1056    compute_inv_mdcts(st->mode, shortBlocks, freq, transient_time, transient_shift, st->out_mem);
1057
1058    for (c=0;c<C;c++)
1059    {
1060       int j;
1061       for (j=0;j<N;j++)
1062       {
1063          celt_sig_t tmp = MAC16_32_Q15(st->out_mem[C*(MAX_PERIOD-N)+C*j+c],
1064                                 preemph,st->preemph_memD[c]);
1065          st->preemph_memD[c] = tmp;
1066          pcm[C*j+c] = SCALEOUT(SIG2WORD16(tmp));
1067       }
1068    }
1069
1070    {
1071       unsigned int val = 0;
1072       while (ec_dec_tell(&dec, 0) < len*8)
1073       {
1074          if (ec_dec_uint(&dec, 2) != val)
1075          {
1076             celt_warning("decode error");
1077             RESTORE_STACK;
1078             return CELT_CORRUPTED_DATA;
1079          }
1080          val = 1-val;
1081       }
1082    }
1083
1084    RESTORE_STACK;
1085    return 0;
1086    /*printf ("\n");*/
1087 }
1088
1089 #ifdef FIXED_POINT
1090 #ifndef DISABLE_FLOAT_API
1091 int celt_decode_float(CELTDecoder * restrict st, unsigned char *data, int len, float * restrict pcm)
1092 {
1093    int j, ret;
1094    const int C = CHANNELS(st->mode);
1095    const int N = st->block_size;
1096    VARDECL(celt_int16_t, out);
1097    SAVE_STACK;
1098    ALLOC(out, C*N, celt_int16_t);
1099
1100    ret=celt_decode(st, data, len, out);
1101
1102    for (j=0;j<C*N;j++)
1103      pcm[j]=out[j]*(1/32768.);
1104    RESTORE_STACK;
1105    return ret;
1106 }
1107 #endif /*DISABLE_FLOAT_API*/
1108 #else
1109 int celt_decode(CELTDecoder * restrict st, unsigned char *data, int len, celt_int16_t * restrict pcm)
1110 {
1111    int j, ret;
1112    VARDECL(celt_sig_t, out);
1113    const int C = CHANNELS(st->mode);
1114    const int N = st->block_size;
1115    SAVE_STACK;
1116    ALLOC(out, C*N, celt_sig_t);
1117
1118    ret=celt_decode_float(st, data, len, out);
1119
1120    for (j=0;j<C*N;j++)
1121      pcm[j] = FLOAT2INT16 (out[j]);
1122
1123    RESTORE_STACK;
1124    return ret;
1125 }
1126 #endif