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