887574b3a9dcf4a97df764416f0f9ec4c237d022
[opus.git] / celt / celt_decoder.c
1 /* Copyright (c) 2007-2008 CSIRO
2    Copyright (c) 2007-2010 Xiph.Org Foundation
3    Copyright (c) 2008 Gregory Maxwell
4    Written by Jean-Marc Valin and Gregory Maxwell */
5 /*
6    Redistribution and use in source and binary forms, with or without
7    modification, are permitted provided that the following conditions
8    are met:
9
10    - Redistributions of source code must retain the above copyright
11    notice, this list of conditions and the following disclaimer.
12
13    - Redistributions in binary form must reproduce the above copyright
14    notice, this list of conditions and the following disclaimer in the
15    documentation and/or other materials provided with the distribution.
16
17    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
18    ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
19    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
20    A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
21    OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
22    EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
23    PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
24    PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
25    LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
26    NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
27    SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
28 */
29
30 #ifdef HAVE_CONFIG_H
31 #include "config.h"
32 #endif
33
34 #define CELT_DECODER_C
35
36 #include "cpu_support.h"
37 #include "os_support.h"
38 #include "mdct.h"
39 #include <math.h>
40 #include "celt.h"
41 #include "pitch.h"
42 #include "bands.h"
43 #include "modes.h"
44 #include "entcode.h"
45 #include "quant_bands.h"
46 #include "rate.h"
47 #include "stack_alloc.h"
48 #include "mathops.h"
49 #include "float_cast.h"
50 #include <stdarg.h>
51 #include "celt_lpc.h"
52 #include "vq.h"
53
54 /**********************************************************************/
55 /*                                                                    */
56 /*                             DECODER                                */
57 /*                                                                    */
58 /**********************************************************************/
59 #define DECODE_BUFFER_SIZE 2048
60
61 /** Decoder state
62  @brief Decoder state
63  */
64 struct OpusCustomDecoder {
65    const OpusCustomMode *mode;
66    int overlap;
67    int channels;
68    int stream_channels;
69
70    int downsample;
71    int start, end;
72    int signalling;
73    int arch;
74
75    /* Everything beyond this point gets cleared on a reset */
76 #define DECODER_RESET_START rng
77
78    opus_uint32 rng;
79    int error;
80    int last_pitch_index;
81    int loss_count;
82    int postfilter_period;
83    int postfilter_period_old;
84    opus_val16 postfilter_gain;
85    opus_val16 postfilter_gain_old;
86    int postfilter_tapset;
87    int postfilter_tapset_old;
88
89    celt_sig preemph_memD[2];
90
91    celt_sig _decode_mem[1]; /* Size = channels*(DECODE_BUFFER_SIZE+mode->overlap) */
92    /* opus_val16 lpc[],  Size = channels*LPC_ORDER */
93    /* opus_val16 oldEBands[], Size = 2*mode->nbEBands */
94    /* opus_val16 oldLogE[], Size = 2*mode->nbEBands */
95    /* opus_val16 oldLogE2[], Size = 2*mode->nbEBands */
96    /* opus_val16 backgroundLogE[], Size = 2*mode->nbEBands */
97 };
98
99 int celt_decoder_get_size(int channels)
100 {
101    const CELTMode *mode = opus_custom_mode_create(48000, 960, NULL);
102    return opus_custom_decoder_get_size(mode, channels);
103 }
104
105 OPUS_CUSTOM_NOSTATIC int opus_custom_decoder_get_size(const CELTMode *mode, int channels)
106 {
107    int size = sizeof(struct CELTDecoder)
108             + (channels*(DECODE_BUFFER_SIZE+mode->overlap)-1)*sizeof(celt_sig)
109             + channels*LPC_ORDER*sizeof(opus_val16)
110             + 4*2*mode->nbEBands*sizeof(opus_val16);
111    return size;
112 }
113
114 #ifdef CUSTOM_MODES
115 CELTDecoder *opus_custom_decoder_create(const CELTMode *mode, int channels, int *error)
116 {
117    int ret;
118    CELTDecoder *st = (CELTDecoder *)opus_alloc(opus_custom_decoder_get_size(mode, channels));
119    ret = opus_custom_decoder_init(st, mode, channels);
120    if (ret != OPUS_OK)
121    {
122       opus_custom_decoder_destroy(st);
123       st = NULL;
124    }
125    if (error)
126       *error = ret;
127    return st;
128 }
129 #endif /* CUSTOM_MODES */
130
131 int celt_decoder_init(CELTDecoder *st, opus_int32 sampling_rate, int channels)
132 {
133    int ret;
134    ret = opus_custom_decoder_init(st, opus_custom_mode_create(48000, 960, NULL), channels);
135    if (ret != OPUS_OK)
136       return ret;
137    st->downsample = resampling_factor(sampling_rate);
138    if (st->downsample==0)
139       return OPUS_BAD_ARG;
140    else
141       return OPUS_OK;
142 }
143
144 OPUS_CUSTOM_NOSTATIC int opus_custom_decoder_init(CELTDecoder *st, const CELTMode *mode, int channels)
145 {
146    if (channels < 0 || channels > 2)
147       return OPUS_BAD_ARG;
148
149    if (st==NULL)
150       return OPUS_ALLOC_FAIL;
151
152    OPUS_CLEAR((char*)st, opus_custom_decoder_get_size(mode, channels));
153
154    st->mode = mode;
155    st->overlap = mode->overlap;
156    st->stream_channels = st->channels = channels;
157
158    st->downsample = 1;
159    st->start = 0;
160    st->end = st->mode->effEBands;
161    st->signalling = 1;
162    st->arch = opus_select_arch();
163
164    st->loss_count = 0;
165
166    opus_custom_decoder_ctl(st, OPUS_RESET_STATE);
167
168    return OPUS_OK;
169 }
170
171 #ifdef CUSTOM_MODES
172 void opus_custom_decoder_destroy(CELTDecoder *st)
173 {
174    opus_free(st);
175 }
176 #endif /* CUSTOM_MODES */
177
178 static OPUS_INLINE opus_val16 SIG2WORD16(celt_sig x)
179 {
180 #ifdef FIXED_POINT
181    x = PSHR32(x, SIG_SHIFT);
182    x = MAX32(x, -32768);
183    x = MIN32(x, 32767);
184    return EXTRACT16(x);
185 #else
186    return (opus_val16)x;
187 #endif
188 }
189
190 #ifndef RESYNTH
191 static
192 #endif
193 void deemphasis(celt_sig *in[], opus_val16 *pcm, int N, int C, int downsample, const opus_val16 *coef,
194       celt_sig *mem, int accum)
195 {
196    int c;
197    int Nd;
198    int apply_downsampling=0;
199    opus_val16 coef0;
200    VARDECL(celt_sig, scratch);
201    SAVE_STACK;
202 #ifndef FIXED_POINT
203    (void)accum;
204    celt_assert(accum==0);
205 #endif
206    ALLOC(scratch, N, celt_sig);
207    coef0 = coef[0];
208    Nd = N/downsample;
209    c=0; do {
210       int j;
211       celt_sig * OPUS_RESTRICT x;
212       opus_val16  * OPUS_RESTRICT y;
213       celt_sig m = mem[c];
214       x =in[c];
215       y = pcm+c;
216 #ifdef CUSTOM_MODES
217       if (coef[1] != 0)
218       {
219          opus_val16 coef1 = coef[1];
220          opus_val16 coef3 = coef[3];
221          for (j=0;j<N;j++)
222          {
223             celt_sig tmp = x[j] + m + VERY_SMALL;
224             m = MULT16_32_Q15(coef0, tmp)
225                           - MULT16_32_Q15(coef1, x[j]);
226             tmp = SHL32(MULT16_32_Q15(coef3, tmp), 2);
227             scratch[j] = tmp;
228          }
229          apply_downsampling=1;
230       } else
231 #endif
232       if (downsample>1)
233       {
234          /* Shortcut for the standard (non-custom modes) case */
235          for (j=0;j<N;j++)
236          {
237             celt_sig tmp = x[j] + m + VERY_SMALL;
238             m = MULT16_32_Q15(coef0, tmp);
239             scratch[j] = tmp;
240          }
241          apply_downsampling=1;
242       } else {
243          /* Shortcut for the standard (non-custom modes) case */
244 #ifdef FIXED_POINT
245          if (accum)
246          {
247             for (j=0;j<N;j++)
248             {
249                celt_sig tmp = x[j] + m + VERY_SMALL;
250                m = MULT16_32_Q15(coef0, tmp);
251                y[j*C] = SAT16(ADD32(y[j*C], SCALEOUT(SIG2WORD16(tmp))));
252             }
253          } else
254 #endif
255          {
256             for (j=0;j<N;j++)
257             {
258                celt_sig tmp = x[j] + m + VERY_SMALL;
259                m = MULT16_32_Q15(coef0, tmp);
260                y[j*C] = SCALEOUT(SIG2WORD16(tmp));
261             }
262          }
263       }
264       mem[c] = m;
265
266       if (apply_downsampling)
267       {
268          /* Perform down-sampling */
269 #ifdef FIXED_POINT
270          if (accum)
271          {
272             for (j=0;j<Nd;j++)
273                y[j*C] = SAT16(ADD32(y[j*C], SCALEOUT(SIG2WORD16(scratch[j*downsample]))));
274          } else
275 #endif
276          {
277             for (j=0;j<Nd;j++)
278                y[j*C] = SCALEOUT(SIG2WORD16(scratch[j*downsample]));
279          }
280       }
281    } while (++c<C);
282    RESTORE_STACK;
283 }
284
285 #ifndef RESYNTH
286 static
287 #endif
288 void celt_synthesis(const CELTMode *mode, celt_norm *X, celt_sig * out_syn[],
289       opus_val16 *oldBandE, int start, int effEnd, int C, int CC, int isTransient,
290       int LM, int downsample, int silence)
291 {
292    int c, i;
293    int M;
294    int b;
295    int B;
296    int N, NB;
297    int shift;
298    int nbEBands;
299    int overlap;
300    VARDECL(celt_sig, freq);
301    SAVE_STACK;
302
303    overlap = mode->overlap;
304    nbEBands = mode->nbEBands;
305    N = mode->shortMdctSize<<LM;
306    ALLOC(freq, N, celt_sig); /**< Interleaved signal MDCTs */
307    M = 1<<LM;
308
309    if (isTransient)
310    {
311       B = M;
312       NB = mode->shortMdctSize;
313       shift = mode->maxLM;
314    } else {
315       B = 1;
316       NB = mode->shortMdctSize<<LM;
317       shift = mode->maxLM-LM;
318    }
319
320    if (CC==2&&C==1)
321    {
322       /* Copying a mono streams to two channels */
323       celt_sig *freq2;
324       denormalise_bands(mode, X, freq, oldBandE, start, effEnd, M,
325             downsample, silence);
326       /* Store a temporary copy in the output buffer because the IMDCT destroys its input. */
327       freq2 = out_syn[1]+overlap/2;
328       OPUS_COPY(freq2, freq, N);
329       for (b=0;b<B;b++)
330          clt_mdct_backward(&mode->mdct, &freq2[b], out_syn[0]+NB*b, mode->window, overlap, shift, B);
331       for (b=0;b<B;b++)
332          clt_mdct_backward(&mode->mdct, &freq[b], out_syn[1]+NB*b, mode->window, overlap, shift, B);
333    } else if (CC==1&&C==2)
334    {
335       /* Downmixing a stereo stream to mono */
336       celt_sig *freq2;
337       freq2 = out_syn[0]+overlap/2;
338       denormalise_bands(mode, X, freq, oldBandE, start, effEnd, M,
339             downsample, silence);
340       /* Use the output buffer as temp array before downmixing. */
341       denormalise_bands(mode, X+N, freq2, oldBandE+nbEBands, start, effEnd, M,
342             downsample, silence);
343       for (i=0;i<N;i++)
344          freq[i] = HALF32(ADD32(freq[i],freq2[i]));
345       for (b=0;b<B;b++)
346          clt_mdct_backward(&mode->mdct, &freq[b], out_syn[0]+NB*b, mode->window, overlap, shift, B);
347    } else {
348       /* Normal case (mono or stereo) */
349       c=0; do {
350          denormalise_bands(mode, X+c*N, freq, oldBandE+c*nbEBands, start, effEnd, M,
351                downsample, silence);
352          for (b=0;b<B;b++)
353             clt_mdct_backward(&mode->mdct, &freq[b], out_syn[c]+NB*b, mode->window, overlap, shift, B);
354       } while (++c<CC);
355    }
356    RESTORE_STACK;
357 }
358
359 static void tf_decode(int start, int end, int isTransient, int *tf_res, int LM, ec_dec *dec)
360 {
361    int i, curr, tf_select;
362    int tf_select_rsv;
363    int tf_changed;
364    int logp;
365    opus_uint32 budget;
366    opus_uint32 tell;
367
368    budget = dec->storage*8;
369    tell = ec_tell(dec);
370    logp = isTransient ? 2 : 4;
371    tf_select_rsv = LM>0 && tell+logp+1<=budget;
372    budget -= tf_select_rsv;
373    tf_changed = curr = 0;
374    for (i=start;i<end;i++)
375    {
376       if (tell+logp<=budget)
377       {
378          curr ^= ec_dec_bit_logp(dec, logp);
379          tell = ec_tell(dec);
380          tf_changed |= curr;
381       }
382       tf_res[i] = curr;
383       logp = isTransient ? 4 : 5;
384    }
385    tf_select = 0;
386    if (tf_select_rsv &&
387      tf_select_table[LM][4*isTransient+0+tf_changed] !=
388      tf_select_table[LM][4*isTransient+2+tf_changed])
389    {
390       tf_select = ec_dec_bit_logp(dec, 1);
391    }
392    for (i=start;i<end;i++)
393    {
394       tf_res[i] = tf_select_table[LM][4*isTransient+2*tf_select+tf_res[i]];
395    }
396 }
397
398 /* The maximum pitch lag to allow in the pitch-based PLC. It's possible to save
399    CPU time in the PLC pitch search by making this smaller than MAX_PERIOD. The
400    current value corresponds to a pitch of 66.67 Hz. */
401 #define PLC_PITCH_LAG_MAX (720)
402 /* The minimum pitch lag to allow in the pitch-based PLC. This corresponds to a
403    pitch of 480 Hz. */
404 #define PLC_PITCH_LAG_MIN (100)
405
406 static int celt_plc_pitch_search(celt_sig *decode_mem[2], int C, int arch)
407 {
408    int pitch_index;
409    SAVE_STACK;
410    VARDECL( opus_val16, lp_pitch_buf );
411    ALLOC( lp_pitch_buf, DECODE_BUFFER_SIZE>>1, opus_val16 );
412    pitch_downsample(decode_mem, lp_pitch_buf,
413          DECODE_BUFFER_SIZE, C, arch);
414    pitch_search(lp_pitch_buf+(PLC_PITCH_LAG_MAX>>1), lp_pitch_buf,
415          DECODE_BUFFER_SIZE-PLC_PITCH_LAG_MAX,
416          PLC_PITCH_LAG_MAX-PLC_PITCH_LAG_MIN, &pitch_index, arch);
417    pitch_index = PLC_PITCH_LAG_MAX-pitch_index;
418    RESTORE_STACK;
419    return pitch_index;
420 }
421
422 static void celt_decode_lost(CELTDecoder * OPUS_RESTRICT st, int N, int LM)
423 {
424    int c;
425    int i;
426    const int C = st->channels;
427    celt_sig *decode_mem[2];
428    celt_sig *out_syn[2];
429    opus_val16 *lpc;
430    opus_val16 *oldBandE, *oldLogE, *oldLogE2, *backgroundLogE;
431    const OpusCustomMode *mode;
432    int nbEBands;
433    int overlap;
434    int start;
435    int loss_count;
436    int noise_based;
437    const opus_int16 *eBands;
438    SAVE_STACK;
439
440    mode = st->mode;
441    nbEBands = mode->nbEBands;
442    overlap = mode->overlap;
443    eBands = mode->eBands;
444
445    c=0; do {
446       decode_mem[c] = st->_decode_mem + c*(DECODE_BUFFER_SIZE+overlap);
447       out_syn[c] = decode_mem[c]+DECODE_BUFFER_SIZE-N;
448    } while (++c<C);
449    lpc = (opus_val16*)(st->_decode_mem+(DECODE_BUFFER_SIZE+overlap)*C);
450    oldBandE = lpc+C*LPC_ORDER;
451    oldLogE = oldBandE + 2*nbEBands;
452    oldLogE2 = oldLogE + 2*nbEBands;
453    backgroundLogE = oldLogE2  + 2*nbEBands;
454
455    loss_count = st->loss_count;
456    start = st->start;
457    noise_based = loss_count >= 5 || start != 0;
458    if (noise_based)
459    {
460       /* Noise-based PLC/CNG */
461 #ifdef SMALL_FOOTPRINT
462       celt_norm *X;
463 #else
464       VARDECL(celt_norm, X);
465 #endif
466       opus_uint32 seed;
467       opus_val16 *plcLogE;
468       int end;
469       int effEnd;
470
471       end = st->end;
472       effEnd = IMAX(start, IMIN(end, mode->effEBands));
473
474 #ifdef SMALL_FOOTPRINT
475       /* This is an ugly hack that breaks aliasing rules and would be easily broken,
476          but it saves almost 4kB of stack. */
477       X = (celt_norm*)(out_syn[C-1]+overlap/2);
478 #else
479       ALLOC(X, C*N, celt_norm);   /**< Interleaved normalised MDCTs */
480 #endif
481
482       if (loss_count >= 5)
483          plcLogE = backgroundLogE;
484       else {
485          /* Energy decay */
486          opus_val16 decay = loss_count==0 ?
487                QCONST16(1.5f, DB_SHIFT) : QCONST16(.5f, DB_SHIFT);
488          c=0; do
489          {
490             for (i=start;i<end;i++)
491                oldBandE[c*nbEBands+i] -= decay;
492          } while (++c<C);
493          plcLogE = oldBandE;
494       }
495       seed = st->rng;
496       for (c=0;c<C;c++)
497       {
498          for (i=start;i<effEnd;i++)
499          {
500             int j;
501             int boffs;
502             int blen;
503             boffs = N*c+(eBands[i]<<LM);
504             blen = (eBands[i+1]-eBands[i])<<LM;
505             for (j=0;j<blen;j++)
506             {
507                seed = celt_lcg_rand(seed);
508                X[boffs+j] = (celt_norm)((opus_int32)seed>>20);
509             }
510             renormalise_vector(X+boffs, blen, Q15ONE);
511          }
512       }
513       st->rng = seed;
514
515       c=0; do {
516          OPUS_MOVE(decode_mem[c], decode_mem[c]+N,
517                DECODE_BUFFER_SIZE-N+(overlap>>1));
518       } while (++c<C);
519
520       celt_synthesis(mode, X, out_syn, plcLogE, start, effEnd, C, C, 0, LM, st->downsample, 0);
521    } else {
522       /* Pitch-based PLC */
523       const opus_val16 *window;
524       opus_val16 fade = Q15ONE;
525       int pitch_index;
526       VARDECL(opus_val32, etmp);
527       VARDECL(opus_val16, exc);
528
529       if (loss_count == 0)
530       {
531          st->last_pitch_index = pitch_index = celt_plc_pitch_search(decode_mem, C, st->arch);
532       } else {
533          pitch_index = st->last_pitch_index;
534          fade = QCONST16(.8f,15);
535       }
536
537       ALLOC(etmp, overlap, opus_val32);
538       ALLOC(exc, MAX_PERIOD, opus_val16);
539       window = mode->window;
540       c=0; do {
541          opus_val16 decay;
542          opus_val16 attenuation;
543          opus_val32 S1=0;
544          celt_sig *buf;
545          int extrapolation_offset;
546          int extrapolation_len;
547          int exc_length;
548          int j;
549
550          buf = decode_mem[c];
551          for (i=0;i<MAX_PERIOD;i++) {
552             exc[i] = ROUND16(buf[DECODE_BUFFER_SIZE-MAX_PERIOD+i], SIG_SHIFT);
553          }
554
555          if (loss_count == 0)
556          {
557             opus_val32 ac[LPC_ORDER+1];
558             /* Compute LPC coefficients for the last MAX_PERIOD samples before
559                the first loss so we can work in the excitation-filter domain. */
560             _celt_autocorr(exc, ac, window, overlap,
561                    LPC_ORDER, MAX_PERIOD, st->arch);
562             /* Add a noise floor of -40 dB. */
563 #ifdef FIXED_POINT
564             ac[0] += SHR32(ac[0],13);
565 #else
566             ac[0] *= 1.0001f;
567 #endif
568             /* Use lag windowing to stabilize the Levinson-Durbin recursion. */
569             for (i=1;i<=LPC_ORDER;i++)
570             {
571                /*ac[i] *= exp(-.5*(2*M_PI*.002*i)*(2*M_PI*.002*i));*/
572 #ifdef FIXED_POINT
573                ac[i] -= MULT16_32_Q15(2*i*i, ac[i]);
574 #else
575                ac[i] -= ac[i]*(0.008f*0.008f)*i*i;
576 #endif
577             }
578             _celt_lpc(lpc+c*LPC_ORDER, ac, LPC_ORDER);
579          }
580          /* We want the excitation for 2 pitch periods in order to look for a
581             decaying signal, but we can't get more than MAX_PERIOD. */
582          exc_length = IMIN(2*pitch_index, MAX_PERIOD);
583          /* Initialize the LPC history with the samples just before the start
584             of the region for which we're computing the excitation. */
585          {
586             opus_val16 lpc_mem[LPC_ORDER];
587             for (i=0;i<LPC_ORDER;i++)
588             {
589                lpc_mem[i] =
590                      ROUND16(buf[DECODE_BUFFER_SIZE-exc_length-1-i], SIG_SHIFT);
591             }
592             /* Compute the excitation for exc_length samples before the loss. */
593             celt_fir(exc+MAX_PERIOD-exc_length, lpc+c*LPC_ORDER,
594                   exc+MAX_PERIOD-exc_length, exc_length, LPC_ORDER, lpc_mem);
595          }
596
597          /* Check if the waveform is decaying, and if so how fast.
598             We do this to avoid adding energy when concealing in a segment
599             with decaying energy. */
600          {
601             opus_val32 E1=1, E2=1;
602             int decay_length;
603 #ifdef FIXED_POINT
604             int shift = IMAX(0,2*celt_zlog2(celt_maxabs16(&exc[MAX_PERIOD-exc_length], exc_length))-20);
605 #endif
606             decay_length = exc_length>>1;
607             for (i=0;i<decay_length;i++)
608             {
609                opus_val16 e;
610                e = exc[MAX_PERIOD-decay_length+i];
611                E1 += SHR32(MULT16_16(e, e), shift);
612                e = exc[MAX_PERIOD-2*decay_length+i];
613                E2 += SHR32(MULT16_16(e, e), shift);
614             }
615             E1 = MIN32(E1, E2);
616             decay = celt_sqrt(frac_div32(SHR32(E1, 1), E2));
617          }
618
619          /* Move the decoder memory one frame to the left to give us room to
620             add the data for the new frame. We ignore the overlap that extends
621             past the end of the buffer, because we aren't going to use it. */
622          OPUS_MOVE(buf, buf+N, DECODE_BUFFER_SIZE-N);
623
624          /* Extrapolate from the end of the excitation with a period of
625             "pitch_index", scaling down each period by an additional factor of
626             "decay". */
627          extrapolation_offset = MAX_PERIOD-pitch_index;
628          /* We need to extrapolate enough samples to cover a complete MDCT
629             window (including overlap/2 samples on both sides). */
630          extrapolation_len = N+overlap;
631          /* We also apply fading if this is not the first loss. */
632          attenuation = MULT16_16_Q15(fade, decay);
633          for (i=j=0;i<extrapolation_len;i++,j++)
634          {
635             opus_val16 tmp;
636             if (j >= pitch_index) {
637                j -= pitch_index;
638                attenuation = MULT16_16_Q15(attenuation, decay);
639             }
640             buf[DECODE_BUFFER_SIZE-N+i] =
641                   SHL32(EXTEND32(MULT16_16_Q15(attenuation,
642                         exc[extrapolation_offset+j])), SIG_SHIFT);
643             /* Compute the energy of the previously decoded signal whose
644                excitation we're copying. */
645             tmp = ROUND16(
646                   buf[DECODE_BUFFER_SIZE-MAX_PERIOD-N+extrapolation_offset+j],
647                   SIG_SHIFT);
648             S1 += SHR32(MULT16_16(tmp, tmp), 8);
649          }
650
651          {
652             opus_val16 lpc_mem[LPC_ORDER];
653             /* Copy the last decoded samples (prior to the overlap region) to
654                synthesis filter memory so we can have a continuous signal. */
655             for (i=0;i<LPC_ORDER;i++)
656                lpc_mem[i] = ROUND16(buf[DECODE_BUFFER_SIZE-N-1-i], SIG_SHIFT);
657             /* Apply the synthesis filter to convert the excitation back into
658                the signal domain. */
659             celt_iir(buf+DECODE_BUFFER_SIZE-N, lpc+c*LPC_ORDER,
660                   buf+DECODE_BUFFER_SIZE-N, extrapolation_len, LPC_ORDER,
661                   lpc_mem);
662          }
663
664          /* Check if the synthesis energy is higher than expected, which can
665             happen with the signal changes during our window. If so,
666             attenuate. */
667          {
668             opus_val32 S2=0;
669             for (i=0;i<extrapolation_len;i++)
670             {
671                opus_val16 tmp = ROUND16(buf[DECODE_BUFFER_SIZE-N+i], SIG_SHIFT);
672                S2 += SHR32(MULT16_16(tmp, tmp), 8);
673             }
674             /* This checks for an "explosion" in the synthesis. */
675 #ifdef FIXED_POINT
676             if (!(S1 > SHR32(S2,2)))
677 #else
678             /* The float test is written this way to catch NaNs in the output
679                of the IIR filter at the same time. */
680             if (!(S1 > 0.2f*S2))
681 #endif
682             {
683                for (i=0;i<extrapolation_len;i++)
684                   buf[DECODE_BUFFER_SIZE-N+i] = 0;
685             } else if (S1 < S2)
686             {
687                opus_val16 ratio = celt_sqrt(frac_div32(SHR32(S1,1)+1,S2+1));
688                for (i=0;i<overlap;i++)
689                {
690                   opus_val16 tmp_g = Q15ONE
691                         - MULT16_16_Q15(window[i], Q15ONE-ratio);
692                   buf[DECODE_BUFFER_SIZE-N+i] =
693                         MULT16_32_Q15(tmp_g, buf[DECODE_BUFFER_SIZE-N+i]);
694                }
695                for (i=overlap;i<extrapolation_len;i++)
696                {
697                   buf[DECODE_BUFFER_SIZE-N+i] =
698                         MULT16_32_Q15(ratio, buf[DECODE_BUFFER_SIZE-N+i]);
699                }
700             }
701          }
702
703          /* Apply the pre-filter to the MDCT overlap for the next frame because
704             the post-filter will be re-applied in the decoder after the MDCT
705             overlap. */
706          comb_filter(etmp, buf+DECODE_BUFFER_SIZE,
707               st->postfilter_period, st->postfilter_period, overlap,
708               -st->postfilter_gain, -st->postfilter_gain,
709               st->postfilter_tapset, st->postfilter_tapset, NULL, 0);
710
711          /* Simulate TDAC on the concealed audio so that it blends with the
712             MDCT of the next frame. */
713          for (i=0;i<overlap/2;i++)
714          {
715             buf[DECODE_BUFFER_SIZE+i] =
716                MULT16_32_Q15(window[i], etmp[overlap-1-i])
717                + MULT16_32_Q15(window[overlap-i-1], etmp[i]);
718          }
719       } while (++c<C);
720    }
721
722    st->loss_count = loss_count+1;
723
724    RESTORE_STACK;
725 }
726
727 int celt_decode_with_ec(CELTDecoder * OPUS_RESTRICT st, const unsigned char *data,
728       int len, opus_val16 * OPUS_RESTRICT pcm, int frame_size, ec_dec *dec, int accum)
729 {
730    int c, i, N;
731    int spread_decision;
732    opus_int32 bits;
733    ec_dec _dec;
734 #ifdef SMALL_FOOTPRINT
735    celt_norm *X;
736 #else
737    VARDECL(celt_norm, X);
738 #endif
739    VARDECL(int, fine_quant);
740    VARDECL(int, pulses);
741    VARDECL(int, cap);
742    VARDECL(int, offsets);
743    VARDECL(int, fine_priority);
744    VARDECL(int, tf_res);
745    VARDECL(unsigned char, collapse_masks);
746    celt_sig *decode_mem[2];
747    celt_sig *out_syn[2];
748    opus_val16 *lpc;
749    opus_val16 *oldBandE, *oldLogE, *oldLogE2, *backgroundLogE;
750
751    int shortBlocks;
752    int isTransient;
753    int intra_ener;
754    const int CC = st->channels;
755    int LM, M;
756    int start;
757    int end;
758    int effEnd;
759    int codedBands;
760    int alloc_trim;
761    int postfilter_pitch;
762    opus_val16 postfilter_gain;
763    int intensity=0;
764    int dual_stereo=0;
765    opus_int32 total_bits;
766    opus_int32 balance;
767    opus_int32 tell;
768    int dynalloc_logp;
769    int postfilter_tapset;
770    int anti_collapse_rsv;
771    int anti_collapse_on=0;
772    int silence;
773    int C = st->stream_channels;
774    const OpusCustomMode *mode;
775    int nbEBands;
776    int overlap;
777    const opus_int16 *eBands;
778    ALLOC_STACK;
779
780    mode = st->mode;
781    nbEBands = mode->nbEBands;
782    overlap = mode->overlap;
783    eBands = mode->eBands;
784    start = st->start;
785    end = st->end;
786    frame_size *= st->downsample;
787
788    lpc = (opus_val16*)(st->_decode_mem+(DECODE_BUFFER_SIZE+overlap)*CC);
789    oldBandE = lpc+CC*LPC_ORDER;
790    oldLogE = oldBandE + 2*nbEBands;
791    oldLogE2 = oldLogE + 2*nbEBands;
792    backgroundLogE = oldLogE2  + 2*nbEBands;
793
794 #ifdef CUSTOM_MODES
795    if (st->signalling && data!=NULL)
796    {
797       int data0=data[0];
798       /* Convert "standard mode" to Opus header */
799       if (mode->Fs==48000 && mode->shortMdctSize==120)
800       {
801          data0 = fromOpus(data0);
802          if (data0<0)
803             return OPUS_INVALID_PACKET;
804       }
805       st->end = end = IMAX(1, mode->effEBands-2*(data0>>5));
806       LM = (data0>>3)&0x3;
807       C = 1 + ((data0>>2)&0x1);
808       data++;
809       len--;
810       if (LM>mode->maxLM)
811          return OPUS_INVALID_PACKET;
812       if (frame_size < mode->shortMdctSize<<LM)
813          return OPUS_BUFFER_TOO_SMALL;
814       else
815          frame_size = mode->shortMdctSize<<LM;
816    } else {
817 #else
818    {
819 #endif
820       for (LM=0;LM<=mode->maxLM;LM++)
821          if (mode->shortMdctSize<<LM==frame_size)
822             break;
823       if (LM>mode->maxLM)
824          return OPUS_BAD_ARG;
825    }
826    M=1<<LM;
827
828    if (len<0 || len>1275 || pcm==NULL)
829       return OPUS_BAD_ARG;
830
831    N = M*mode->shortMdctSize;
832    c=0; do {
833       decode_mem[c] = st->_decode_mem + c*(DECODE_BUFFER_SIZE+overlap);
834       out_syn[c] = decode_mem[c]+DECODE_BUFFER_SIZE-N;
835    } while (++c<CC);
836
837    effEnd = end;
838    if (effEnd > mode->effEBands)
839       effEnd = mode->effEBands;
840
841    if (data == NULL || len<=1)
842    {
843       celt_decode_lost(st, N, LM);
844       deemphasis(out_syn, pcm, N, CC, st->downsample, mode->preemph, st->preemph_memD, accum);
845       RESTORE_STACK;
846       return frame_size/st->downsample;
847    }
848
849    if (dec == NULL)
850    {
851       ec_dec_init(&_dec,(unsigned char*)data,len);
852       dec = &_dec;
853    }
854
855    if (C==1)
856    {
857       for (i=0;i<nbEBands;i++)
858          oldBandE[i]=MAX16(oldBandE[i],oldBandE[nbEBands+i]);
859    }
860
861    total_bits = len*8;
862    tell = ec_tell(dec);
863
864    if (tell >= total_bits)
865       silence = 1;
866    else if (tell==1)
867       silence = ec_dec_bit_logp(dec, 15);
868    else
869       silence = 0;
870    if (silence)
871    {
872       /* Pretend we've read all the remaining bits */
873       tell = len*8;
874       dec->nbits_total+=tell-ec_tell(dec);
875    }
876
877    postfilter_gain = 0;
878    postfilter_pitch = 0;
879    postfilter_tapset = 0;
880    if (start==0 && tell+16 <= total_bits)
881    {
882       if(ec_dec_bit_logp(dec, 1))
883       {
884          int qg, octave;
885          octave = ec_dec_uint(dec, 6);
886          postfilter_pitch = (16<<octave)+ec_dec_bits(dec, 4+octave)-1;
887          qg = ec_dec_bits(dec, 3);
888          if (ec_tell(dec)+2<=total_bits)
889             postfilter_tapset = ec_dec_icdf(dec, tapset_icdf, 2);
890          postfilter_gain = QCONST16(.09375f,15)*(qg+1);
891       }
892       tell = ec_tell(dec);
893    }
894
895    if (LM > 0 && tell+3 <= total_bits)
896    {
897       isTransient = ec_dec_bit_logp(dec, 3);
898       tell = ec_tell(dec);
899    }
900    else
901       isTransient = 0;
902
903    if (isTransient)
904       shortBlocks = M;
905    else
906       shortBlocks = 0;
907
908    /* Decode the global flags (first symbols in the stream) */
909    intra_ener = tell+3<=total_bits ? ec_dec_bit_logp(dec, 3) : 0;
910    /* Get band energies */
911    unquant_coarse_energy(mode, start, end, oldBandE,
912          intra_ener, dec, C, LM);
913
914    ALLOC(tf_res, nbEBands, int);
915    tf_decode(start, end, isTransient, tf_res, LM, dec);
916
917    tell = ec_tell(dec);
918    spread_decision = SPREAD_NORMAL;
919    if (tell+4 <= total_bits)
920       spread_decision = ec_dec_icdf(dec, spread_icdf, 5);
921
922    ALLOC(cap, nbEBands, int);
923
924    init_caps(mode,cap,LM,C);
925
926    ALLOC(offsets, nbEBands, int);
927
928    dynalloc_logp = 6;
929    total_bits<<=BITRES;
930    tell = ec_tell_frac(dec);
931    for (i=start;i<end;i++)
932    {
933       int width, quanta;
934       int dynalloc_loop_logp;
935       int boost;
936       width = C*(eBands[i+1]-eBands[i])<<LM;
937       /* quanta is 6 bits, but no more than 1 bit/sample
938          and no less than 1/8 bit/sample */
939       quanta = IMIN(width<<BITRES, IMAX(6<<BITRES, width));
940       dynalloc_loop_logp = dynalloc_logp;
941       boost = 0;
942       while (tell+(dynalloc_loop_logp<<BITRES) < total_bits && boost < cap[i])
943       {
944          int flag;
945          flag = ec_dec_bit_logp(dec, dynalloc_loop_logp);
946          tell = ec_tell_frac(dec);
947          if (!flag)
948             break;
949          boost += quanta;
950          total_bits -= quanta;
951          dynalloc_loop_logp = 1;
952       }
953       offsets[i] = boost;
954       /* Making dynalloc more likely */
955       if (boost>0)
956          dynalloc_logp = IMAX(2, dynalloc_logp-1);
957    }
958
959    ALLOC(fine_quant, nbEBands, int);
960    alloc_trim = tell+(6<<BITRES) <= total_bits ?
961          ec_dec_icdf(dec, trim_icdf, 7) : 5;
962
963    bits = (((opus_int32)len*8)<<BITRES) - ec_tell_frac(dec) - 1;
964    anti_collapse_rsv = isTransient&&LM>=2&&bits>=((LM+2)<<BITRES) ? (1<<BITRES) : 0;
965    bits -= anti_collapse_rsv;
966
967    ALLOC(pulses, nbEBands, int);
968    ALLOC(fine_priority, nbEBands, int);
969
970    codedBands = compute_allocation(mode, start, end, offsets, cap,
971          alloc_trim, &intensity, &dual_stereo, bits, &balance, pulses,
972          fine_quant, fine_priority, C, LM, dec, 0, 0, 0);
973
974    unquant_fine_energy(mode, start, end, oldBandE, fine_quant, dec, C);
975
976    c=0; do {
977       OPUS_MOVE(decode_mem[c], decode_mem[c]+N, DECODE_BUFFER_SIZE-N+overlap/2);
978    } while (++c<CC);
979
980    /* Decode fixed codebook */
981    ALLOC(collapse_masks, C*nbEBands, unsigned char);
982
983 #ifdef SMALL_FOOTPRINT
984    /* This is an ugly hack that breaks aliasing rules and would be easily broken,
985       but it saves almost 4kB of stack. */
986    X = (celt_norm*)(out_syn[CC-1]+overlap/2);
987 #else
988    ALLOC(X, C*N, celt_norm);   /**< Interleaved normalised MDCTs */
989 #endif
990
991    quant_all_bands(0, mode, start, end, X, C==2 ? X+N : NULL, collapse_masks,
992          NULL, pulses, shortBlocks, spread_decision, dual_stereo, intensity, tf_res,
993          len*(8<<BITRES)-anti_collapse_rsv, balance, dec, LM, codedBands, &st->rng);
994
995    if (anti_collapse_rsv > 0)
996    {
997       anti_collapse_on = ec_dec_bits(dec, 1);
998    }
999
1000    unquant_energy_finalise(mode, start, end, oldBandE,
1001          fine_quant, fine_priority, len*8-ec_tell(dec), dec, C);
1002
1003    if (anti_collapse_on)
1004       anti_collapse(mode, X, collapse_masks, LM, C, N,
1005             start, end, oldBandE, oldLogE, oldLogE2, pulses, st->rng);
1006
1007    if (silence)
1008    {
1009       for (i=0;i<C*nbEBands;i++)
1010          oldBandE[i] = -QCONST16(28.f,DB_SHIFT);
1011    }
1012
1013    celt_synthesis(mode, X, out_syn, oldBandE, start, effEnd, C, CC, isTransient, LM, st->downsample, silence);
1014
1015    c=0; do {
1016       st->postfilter_period=IMAX(st->postfilter_period, COMBFILTER_MINPERIOD);
1017       st->postfilter_period_old=IMAX(st->postfilter_period_old, COMBFILTER_MINPERIOD);
1018       comb_filter(out_syn[c], out_syn[c], st->postfilter_period_old, st->postfilter_period, mode->shortMdctSize,
1019             st->postfilter_gain_old, st->postfilter_gain, st->postfilter_tapset_old, st->postfilter_tapset,
1020             mode->window, overlap);
1021       if (LM!=0)
1022          comb_filter(out_syn[c]+mode->shortMdctSize, out_syn[c]+mode->shortMdctSize, st->postfilter_period, postfilter_pitch, N-mode->shortMdctSize,
1023                st->postfilter_gain, postfilter_gain, st->postfilter_tapset, postfilter_tapset,
1024                mode->window, overlap);
1025
1026    } while (++c<CC);
1027    st->postfilter_period_old = st->postfilter_period;
1028    st->postfilter_gain_old = st->postfilter_gain;
1029    st->postfilter_tapset_old = st->postfilter_tapset;
1030    st->postfilter_period = postfilter_pitch;
1031    st->postfilter_gain = postfilter_gain;
1032    st->postfilter_tapset = postfilter_tapset;
1033    if (LM!=0)
1034    {
1035       st->postfilter_period_old = st->postfilter_period;
1036       st->postfilter_gain_old = st->postfilter_gain;
1037       st->postfilter_tapset_old = st->postfilter_tapset;
1038    }
1039
1040    if (C==1)
1041       OPUS_COPY(&oldBandE[nbEBands], oldBandE, nbEBands);
1042
1043    /* In case start or end were to change */
1044    if (!isTransient)
1045    {
1046       OPUS_COPY(oldLogE2, oldLogE, 2*nbEBands);
1047       OPUS_COPY(oldLogE, oldBandE, 2*nbEBands);
1048       for (i=0;i<2*nbEBands;i++)
1049          backgroundLogE[i] = MIN16(backgroundLogE[i] + M*QCONST16(0.001f,DB_SHIFT), oldBandE[i]);
1050    } else {
1051       for (i=0;i<2*nbEBands;i++)
1052          oldLogE[i] = MIN16(oldLogE[i], oldBandE[i]);
1053    }
1054    c=0; do
1055    {
1056       for (i=0;i<start;i++)
1057       {
1058          oldBandE[c*nbEBands+i]=0;
1059          oldLogE[c*nbEBands+i]=oldLogE2[c*nbEBands+i]=-QCONST16(28.f,DB_SHIFT);
1060       }
1061       for (i=end;i<nbEBands;i++)
1062       {
1063          oldBandE[c*nbEBands+i]=0;
1064          oldLogE[c*nbEBands+i]=oldLogE2[c*nbEBands+i]=-QCONST16(28.f,DB_SHIFT);
1065       }
1066    } while (++c<2);
1067    st->rng = dec->rng;
1068
1069    deemphasis(out_syn, pcm, N, CC, st->downsample, mode->preemph, st->preemph_memD, accum);
1070    st->loss_count = 0;
1071    RESTORE_STACK;
1072    if (ec_tell(dec) > 8*len)
1073       return OPUS_INTERNAL_ERROR;
1074    if(ec_get_error(dec))
1075       st->error = 1;
1076    return frame_size/st->downsample;
1077 }
1078
1079
1080 #ifdef CUSTOM_MODES
1081
1082 #ifdef FIXED_POINT
1083 int opus_custom_decode(CELTDecoder * OPUS_RESTRICT st, const unsigned char *data, int len, opus_int16 * OPUS_RESTRICT pcm, int frame_size)
1084 {
1085    return celt_decode_with_ec(st, data, len, pcm, frame_size, NULL, 0);
1086 }
1087
1088 #ifndef DISABLE_FLOAT_API
1089 int opus_custom_decode_float(CELTDecoder * OPUS_RESTRICT st, const unsigned char *data, int len, float * OPUS_RESTRICT pcm, int frame_size)
1090 {
1091    int j, ret, C, N;
1092    VARDECL(opus_int16, out);
1093    ALLOC_STACK;
1094
1095    if (pcm==NULL)
1096       return OPUS_BAD_ARG;
1097
1098    C = st->channels;
1099    N = frame_size;
1100
1101    ALLOC(out, C*N, opus_int16);
1102    ret=celt_decode_with_ec(st, data, len, out, frame_size, NULL, 0);
1103    if (ret>0)
1104       for (j=0;j<C*ret;j++)
1105          pcm[j]=out[j]*(1.f/32768.f);
1106
1107    RESTORE_STACK;
1108    return ret;
1109 }
1110 #endif /* DISABLE_FLOAT_API */
1111
1112 #else
1113
1114 int opus_custom_decode_float(CELTDecoder * OPUS_RESTRICT st, const unsigned char *data, int len, float * OPUS_RESTRICT pcm, int frame_size)
1115 {
1116    return celt_decode_with_ec(st, data, len, pcm, frame_size, NULL, 0);
1117 }
1118
1119 int opus_custom_decode(CELTDecoder * OPUS_RESTRICT st, const unsigned char *data, int len, opus_int16 * OPUS_RESTRICT pcm, int frame_size)
1120 {
1121    int j, ret, C, N;
1122    VARDECL(celt_sig, out);
1123    ALLOC_STACK;
1124
1125    if (pcm==NULL)
1126       return OPUS_BAD_ARG;
1127
1128    C = st->channels;
1129    N = frame_size;
1130    ALLOC(out, C*N, celt_sig);
1131
1132    ret=celt_decode_with_ec(st, data, len, out, frame_size, NULL, 0);
1133
1134    if (ret>0)
1135       for (j=0;j<C*ret;j++)
1136          pcm[j] = FLOAT2INT16 (out[j]);
1137
1138    RESTORE_STACK;
1139    return ret;
1140 }
1141
1142 #endif
1143 #endif /* CUSTOM_MODES */
1144
1145 int opus_custom_decoder_ctl(CELTDecoder * OPUS_RESTRICT st, int request, ...)
1146 {
1147    va_list ap;
1148
1149    va_start(ap, request);
1150    switch (request)
1151    {
1152       case CELT_SET_START_BAND_REQUEST:
1153       {
1154          opus_int32 value = va_arg(ap, opus_int32);
1155          if (value<0 || value>=st->mode->nbEBands)
1156             goto bad_arg;
1157          st->start = value;
1158       }
1159       break;
1160       case CELT_SET_END_BAND_REQUEST:
1161       {
1162          opus_int32 value = va_arg(ap, opus_int32);
1163          if (value<1 || value>st->mode->nbEBands)
1164             goto bad_arg;
1165          st->end = value;
1166       }
1167       break;
1168       case CELT_SET_CHANNELS_REQUEST:
1169       {
1170          opus_int32 value = va_arg(ap, opus_int32);
1171          if (value<1 || value>2)
1172             goto bad_arg;
1173          st->stream_channels = value;
1174       }
1175       break;
1176       case CELT_GET_AND_CLEAR_ERROR_REQUEST:
1177       {
1178          opus_int32 *value = va_arg(ap, opus_int32*);
1179          if (value==NULL)
1180             goto bad_arg;
1181          *value=st->error;
1182          st->error = 0;
1183       }
1184       break;
1185       case OPUS_GET_LOOKAHEAD_REQUEST:
1186       {
1187          opus_int32 *value = va_arg(ap, opus_int32*);
1188          if (value==NULL)
1189             goto bad_arg;
1190          *value = st->overlap/st->downsample;
1191       }
1192       break;
1193       case OPUS_RESET_STATE:
1194       {
1195          int i;
1196          opus_val16 *lpc, *oldBandE, *oldLogE, *oldLogE2;
1197          lpc = (opus_val16*)(st->_decode_mem+(DECODE_BUFFER_SIZE+st->overlap)*st->channels);
1198          oldBandE = lpc+st->channels*LPC_ORDER;
1199          oldLogE = oldBandE + 2*st->mode->nbEBands;
1200          oldLogE2 = oldLogE + 2*st->mode->nbEBands;
1201          OPUS_CLEAR((char*)&st->DECODER_RESET_START,
1202                opus_custom_decoder_get_size(st->mode, st->channels)-
1203                ((char*)&st->DECODER_RESET_START - (char*)st));
1204          for (i=0;i<2*st->mode->nbEBands;i++)
1205             oldLogE[i]=oldLogE2[i]=-QCONST16(28.f,DB_SHIFT);
1206       }
1207       break;
1208       case OPUS_GET_PITCH_REQUEST:
1209       {
1210          opus_int32 *value = va_arg(ap, opus_int32*);
1211          if (value==NULL)
1212             goto bad_arg;
1213          *value = st->postfilter_period;
1214       }
1215       break;
1216       case CELT_GET_MODE_REQUEST:
1217       {
1218          const CELTMode ** value = va_arg(ap, const CELTMode**);
1219          if (value==0)
1220             goto bad_arg;
1221          *value=st->mode;
1222       }
1223       break;
1224       case CELT_SET_SIGNALLING_REQUEST:
1225       {
1226          opus_int32 value = va_arg(ap, opus_int32);
1227          st->signalling = value;
1228       }
1229       break;
1230       case OPUS_GET_FINAL_RANGE_REQUEST:
1231       {
1232          opus_uint32 * value = va_arg(ap, opus_uint32 *);
1233          if (value==0)
1234             goto bad_arg;
1235          *value=st->rng;
1236       }
1237       break;
1238       default:
1239          goto bad_request;
1240    }
1241    va_end(ap);
1242    return OPUS_OK;
1243 bad_arg:
1244    va_end(ap);
1245    return OPUS_BAD_ARG;
1246 bad_request:
1247       va_end(ap);
1248   return OPUS_UNIMPLEMENTED;
1249 }