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