Making it easier to Torben to develop his new PLC code
[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 #ifdef NEW_PLC
792 #define DECODE_BUFFER_SIZE 2048
793 #else
794 #define DECODE_BUFFER_SIZE MAX_PERIOD
795 #endif
796
797 /** Decoder state 
798  @brief Decoder state
799  */
800 struct CELTDecoder {
801    const CELTMode *mode;
802    int frame_size;
803    int block_size;
804    int overlap;
805
806    ec_byte_buffer buf;
807    ec_enc         enc;
808
809    celt_sig_t * restrict preemph_memD;
810
811    celt_sig_t *out_mem;
812    celt_sig_t *decode_mem;
813
814    celt_word16_t *oldBandE;
815    
816    int last_pitch_index;
817 };
818
819 CELTDecoder *celt_decoder_create(const CELTMode *mode)
820 {
821    int N, C;
822    CELTDecoder *st;
823
824    if (check_mode(mode) != CELT_OK)
825       return NULL;
826
827    N = mode->mdctSize;
828    C = CHANNELS(mode);
829    st = celt_alloc(sizeof(CELTDecoder));
830    
831    st->mode = mode;
832    st->frame_size = N;
833    st->block_size = N;
834    st->overlap = mode->overlap;
835
836    st->decode_mem = celt_alloc((DECODE_BUFFER_SIZE+st->overlap)*C*sizeof(celt_sig_t));
837    st->out_mem = st->decode_mem+DECODE_BUFFER_SIZE-MAX_PERIOD;
838    
839    st->oldBandE = (celt_word16_t*)celt_alloc(C*mode->nbEBands*sizeof(celt_word16_t));
840
841    st->preemph_memD = (celt_sig_t*)celt_alloc(C*sizeof(celt_sig_t));
842
843    st->last_pitch_index = 0;
844    return st;
845 }
846
847 void celt_decoder_destroy(CELTDecoder *st)
848 {
849    if (st == NULL)
850    {
851       celt_warning("NULL passed to celt_encoder_destroy");
852       return;
853    }
854    if (check_mode(st->mode) != CELT_OK)
855       return;
856
857
858    celt_free(st->decode_mem);
859    
860    celt_free(st->oldBandE);
861    
862    celt_free(st->preemph_memD);
863
864    celt_free(st);
865 }
866
867 /** Handles lost packets by just copying past data with the same offset as the last
868     pitch period */
869 #ifdef NEW_PLC
870 #include "plc.c"
871 #else
872 static void celt_decode_lost(CELTDecoder * restrict st, celt_word16_t * restrict pcm)
873 {
874    int c, N;
875    int pitch_index;
876    int i, len;
877    VARDECL(celt_sig_t, freq);
878    const int C = CHANNELS(st->mode);
879    int offset;
880    SAVE_STACK;
881    N = st->block_size;
882    ALLOC(freq,C*N, celt_sig_t);         /**< Interleaved signal MDCTs */
883    
884    len = N+st->mode->overlap;
885 #if 0
886    pitch_index = st->last_pitch_index;
887    
888    /* Use the pitch MDCT as the "guessed" signal */
889    compute_mdcts(st->mode, st->mode->window, st->out_mem+pitch_index*C, freq);
890
891 #else
892    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);
893    pitch_index = MAX_PERIOD-len-pitch_index;
894    offset = MAX_PERIOD-pitch_index;
895    while (offset+len >= MAX_PERIOD)
896       offset -= pitch_index;
897    compute_mdcts(st->mode, 0, st->out_mem+offset*C, freq);
898    for (i=0;i<N;i++)
899       freq[i] = ADD32(EPSILON, MULT16_32_Q15(QCONST16(.9f,15),freq[i]));
900 #endif
901    
902    
903    
904    CELT_MOVE(st->out_mem, st->out_mem+C*N, C*(MAX_PERIOD+st->mode->overlap-N));
905    /* Compute inverse MDCTs */
906    compute_inv_mdcts(st->mode, 0, freq, -1, 1, st->out_mem);
907
908    for (c=0;c<C;c++)
909    {
910       int j;
911       for (j=0;j<N;j++)
912       {
913          celt_sig_t tmp = MAC16_32_Q15(st->out_mem[C*(MAX_PERIOD-N)+C*j+c],
914                                 preemph,st->preemph_memD[c]);
915          st->preemph_memD[c] = tmp;
916          pcm[C*j+c] = SCALEOUT(SIG2WORD16(tmp));
917       }
918    }
919    RESTORE_STACK;
920 }
921 #endif
922
923 #ifdef FIXED_POINT
924 int celt_decode(CELTDecoder * restrict st, const unsigned char *data, int len, celt_int16_t * restrict pcm)
925 {
926 #else
927 int celt_decode_float(CELTDecoder * restrict st, const unsigned char *data, int len, celt_sig_t * restrict pcm)
928 {
929 #endif
930    int i, c, N, N4;
931    int has_pitch, has_fold;
932    int pitch_index;
933    int bits;
934    ec_dec dec;
935    ec_byte_buffer buf;
936    VARDECL(celt_sig_t, freq);
937    VARDECL(celt_norm_t, X);
938    VARDECL(celt_norm_t, P);
939    VARDECL(celt_ener_t, bandE);
940    VARDECL(celt_pgain_t, gains);
941    VARDECL(int, stereo_mode);
942    VARDECL(int, fine_quant);
943    VARDECL(int, pulses);
944    VARDECL(int, offsets);
945
946    int shortBlocks;
947    int transient_time;
948    int transient_shift;
949    const int C = CHANNELS(st->mode);
950    SAVE_STACK;
951
952    if (check_mode(st->mode) != CELT_OK)
953       return CELT_INVALID_MODE;
954
955    N = st->block_size;
956    N4 = (N-st->overlap)>>1;
957
958    ALLOC(freq, C*N, celt_sig_t); /**< Interleaved signal MDCTs */
959    ALLOC(X, C*N, celt_norm_t);         /**< Interleaved normalised MDCTs */
960    ALLOC(P, C*N, celt_norm_t);         /**< Interleaved normalised pitch MDCTs*/
961    ALLOC(bandE, st->mode->nbEBands*C, celt_ener_t);
962    ALLOC(gains, st->mode->nbPBands, celt_pgain_t);
963    
964    if (check_mode(st->mode) != CELT_OK)
965    {
966       RESTORE_STACK;
967       return CELT_INVALID_MODE;
968    }
969    if (data == NULL)
970    {
971       celt_decode_lost(st, pcm);
972       RESTORE_STACK;
973       return 0;
974    }
975    if (len<0) {
976      RESTORE_STACK;
977      return CELT_BAD_ARG;
978    }
979    
980    ec_byte_readinit(&buf,data,len);
981    ec_dec_init(&dec,&buf);
982    
983    has_pitch = ec_dec_bits(&dec, 1);
984    if (has_pitch)
985    {
986       has_fold = ec_dec_bits(&dec, 1);
987       shortBlocks = 0;
988    } else if (st->mode->nbShortMdcts > 1){
989       shortBlocks = ec_dec_bits(&dec, 1);
990       has_fold = 1;
991    } else {
992       shortBlocks = 0;
993       has_fold = 1;
994    }
995    if (shortBlocks)
996    {
997       transient_shift = ec_dec_bits(&dec, 2);
998       if (transient_shift)
999          transient_time = ec_dec_uint(&dec, N+st->mode->overlap);
1000       else
1001          transient_time = 0;
1002    } else {
1003       transient_time = -1;
1004       transient_shift = 0;
1005    }
1006    
1007    if (has_pitch)
1008    {
1009       pitch_index = ec_dec_uint(&dec, MAX_PERIOD-(2*N-2*N4));
1010       st->last_pitch_index = pitch_index;
1011    } else {
1012       pitch_index = 0;
1013       for (i=0;i<st->mode->nbPBands;i++)
1014          gains[i] = 0;
1015    }
1016
1017    ALLOC(fine_quant, st->mode->nbEBands, int);
1018    /* Get band energies */
1019    unquant_coarse_energy(st->mode, bandE, st->oldBandE, len*8/3, st->mode->prob, &dec);
1020    
1021    ALLOC(pulses, st->mode->nbEBands, int);
1022    ALLOC(offsets, st->mode->nbEBands, int);
1023    ALLOC(stereo_mode, st->mode->nbEBands, int);
1024    stereo_decision(st->mode, X, stereo_mode, st->mode->nbEBands);
1025
1026    for (i=0;i<st->mode->nbEBands;i++)
1027       offsets[i] = 0;
1028
1029    bits = len*8 - ec_dec_tell(&dec, 0) - 1;
1030    if (has_pitch)
1031       bits -= st->mode->nbPBands;
1032    compute_allocation(st->mode, offsets, stereo_mode, bits, pulses, fine_quant);
1033    /*bits = ec_dec_tell(&dec, 0);
1034    compute_fine_allocation(st->mode, fine_quant, (20*C+len*8/5-(ec_dec_tell(&dec, 0)-bits))/C);*/
1035    
1036    unquant_fine_energy(st->mode, bandE, st->oldBandE, fine_quant, &dec);
1037
1038
1039    if (has_pitch) 
1040    {
1041       VARDECL(celt_ener_t, bandEp);
1042       
1043       /* Pitch MDCT */
1044       compute_mdcts(st->mode, 0, st->out_mem+pitch_index*C, freq);
1045       ALLOC(bandEp, st->mode->nbEBands*C, celt_ener_t);
1046       compute_band_energies(st->mode, freq, bandEp);
1047       normalise_bands(st->mode, freq, P, bandEp);
1048       /* Apply pitch gains */
1049    } else {
1050       for (i=0;i<C*N;i++)
1051          P[i] = 0;
1052    }
1053
1054    /* Decode fixed codebook and merge with pitch */
1055    if (C==1)
1056       unquant_bands(st->mode, X, P, has_pitch, gains, bandE, pulses, shortBlocks, has_fold, len*8, &dec);
1057    else
1058       unquant_bands_stereo(st->mode, X, P, has_pitch, gains, bandE, pulses, shortBlocks, has_fold, len*8, &dec);
1059
1060    /* Synthesis */
1061    denormalise_bands(st->mode, X, freq, bandE);
1062
1063
1064    CELT_MOVE(st->decode_mem, st->decode_mem+C*N, C*(DECODE_BUFFER_SIZE+st->overlap-N));
1065    /* Compute inverse MDCTs */
1066    compute_inv_mdcts(st->mode, shortBlocks, freq, transient_time, transient_shift, st->out_mem);
1067
1068    for (c=0;c<C;c++)
1069    {
1070       int j;
1071       for (j=0;j<N;j++)
1072       {
1073          celt_sig_t tmp = MAC16_32_Q15(st->out_mem[C*(MAX_PERIOD-N)+C*j+c],
1074                                 preemph,st->preemph_memD[c]);
1075          st->preemph_memD[c] = tmp;
1076          pcm[C*j+c] = SCALEOUT(SIG2WORD16(tmp));
1077       }
1078    }
1079
1080    {
1081       unsigned int val = 0;
1082       while (ec_dec_tell(&dec, 0) < len*8)
1083       {
1084          if (ec_dec_uint(&dec, 2) != val)
1085          {
1086             celt_warning("decode error");
1087             RESTORE_STACK;
1088             return CELT_CORRUPTED_DATA;
1089          }
1090          val = 1-val;
1091       }
1092    }
1093
1094    RESTORE_STACK;
1095    return 0;
1096    /*printf ("\n");*/
1097 }
1098
1099 #ifdef FIXED_POINT
1100 #ifndef DISABLE_FLOAT_API
1101 int celt_decode_float(CELTDecoder * restrict st, const unsigned char *data, int len, float * restrict pcm)
1102 {
1103    int j, ret;
1104    const int C = CHANNELS(st->mode);
1105    const int N = st->block_size;
1106    VARDECL(celt_int16_t, out);
1107    SAVE_STACK;
1108    ALLOC(out, C*N, celt_int16_t);
1109
1110    ret=celt_decode(st, data, len, out);
1111
1112    for (j=0;j<C*N;j++)
1113      pcm[j]=out[j]*(1/32768.);
1114    RESTORE_STACK;
1115    return ret;
1116 }
1117 #endif /*DISABLE_FLOAT_API*/
1118 #else
1119 int celt_decode(CELTDecoder * restrict st, const unsigned char *data, int len, celt_int16_t * restrict pcm)
1120 {
1121    int j, ret;
1122    VARDECL(celt_sig_t, out);
1123    const int C = CHANNELS(st->mode);
1124    const int N = st->block_size;
1125    SAVE_STACK;
1126    ALLOC(out, C*N, celt_sig_t);
1127
1128    ret=celt_decode_float(st, data, len, out);
1129
1130    for (j=0;j<C*N;j++)
1131      pcm[j] = FLOAT2INT16 (out[j]);
1132
1133    RESTORE_STACK;
1134    return ret;
1135 }
1136 #endif