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