Moves deemphasis() call out of celt_decode_lost() to reduce peak stack
[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 void celt_decode_lost(CELTDecoder * OPUS_RESTRICT st, int N, int LM)
407 {
408    int c;
409    int i;
410    const int C = st->channels;
411    celt_sig *decode_mem[2];
412    celt_sig *out_syn[2];
413    opus_val16 *lpc;
414    opus_val16 *oldBandE, *oldLogE, *oldLogE2, *backgroundLogE;
415    const OpusCustomMode *mode;
416    int nbEBands;
417    int overlap;
418    int start;
419    int loss_count;
420    int noise_based;
421    const opus_int16 *eBands;
422    SAVE_STACK;
423
424    mode = st->mode;
425    nbEBands = mode->nbEBands;
426    overlap = mode->overlap;
427    eBands = mode->eBands;
428
429    c=0; do {
430       decode_mem[c] = st->_decode_mem + c*(DECODE_BUFFER_SIZE+overlap);
431       out_syn[c] = decode_mem[c]+DECODE_BUFFER_SIZE-N;
432    } while (++c<C);
433    lpc = (opus_val16*)(st->_decode_mem+(DECODE_BUFFER_SIZE+overlap)*C);
434    oldBandE = lpc+C*LPC_ORDER;
435    oldLogE = oldBandE + 2*nbEBands;
436    oldLogE2 = oldLogE + 2*nbEBands;
437    backgroundLogE = oldLogE2  + 2*nbEBands;
438
439    loss_count = st->loss_count;
440    start = st->start;
441    noise_based = loss_count >= 5 || start != 0;
442    if (noise_based)
443    {
444       /* Noise-based PLC/CNG */
445       VARDECL(celt_norm, X);
446       opus_uint32 seed;
447       opus_val16 *plcLogE;
448       int end;
449       int effEnd;
450
451       end = st->end;
452       effEnd = IMAX(start, IMIN(end, mode->effEBands));
453
454       /* Share the interleaved signal MDCT coefficient buffer with the
455          deemphasis scratch buffer. */
456       ALLOC(X, C*N, celt_norm);   /**< Interleaved normalised MDCTs */
457
458       if (loss_count >= 5)
459          plcLogE = backgroundLogE;
460       else {
461          /* Energy decay */
462          opus_val16 decay = loss_count==0 ?
463                QCONST16(1.5f, DB_SHIFT) : QCONST16(.5f, DB_SHIFT);
464          c=0; do
465          {
466             for (i=start;i<end;i++)
467                oldBandE[c*nbEBands+i] -= decay;
468          } while (++c<C);
469          plcLogE = oldBandE;
470       }
471       seed = st->rng;
472       for (c=0;c<C;c++)
473       {
474          for (i=start;i<effEnd;i++)
475          {
476             int j;
477             int boffs;
478             int blen;
479             boffs = N*c+(eBands[i]<<LM);
480             blen = (eBands[i+1]-eBands[i])<<LM;
481             for (j=0;j<blen;j++)
482             {
483                seed = celt_lcg_rand(seed);
484                X[boffs+j] = (celt_norm)((opus_int32)seed>>20);
485             }
486             renormalise_vector(X+boffs, blen, Q15ONE);
487          }
488       }
489       st->rng = seed;
490
491       c=0; do {
492          OPUS_MOVE(decode_mem[c], decode_mem[c]+N,
493                DECODE_BUFFER_SIZE-N+(overlap>>1));
494       } while (++c<C);
495
496       celt_synthesis(mode, X, out_syn, plcLogE, start, effEnd, C, C, 0, LM, st->downsample, 0);
497    } else {
498       /* Pitch-based PLC */
499       const opus_val16 *window;
500       opus_val16 fade = Q15ONE;
501       int pitch_index;
502       VARDECL(opus_val32, etmp);
503       VARDECL(opus_val16, exc);
504
505       if (loss_count == 0)
506       {
507          VARDECL( opus_val16, lp_pitch_buf );
508          ALLOC( lp_pitch_buf, DECODE_BUFFER_SIZE>>1, opus_val16 );
509          pitch_downsample(decode_mem, lp_pitch_buf,
510                DECODE_BUFFER_SIZE, C, st->arch);
511          pitch_search(lp_pitch_buf+(PLC_PITCH_LAG_MAX>>1), lp_pitch_buf,
512                DECODE_BUFFER_SIZE-PLC_PITCH_LAG_MAX,
513                PLC_PITCH_LAG_MAX-PLC_PITCH_LAG_MIN, &pitch_index, st->arch);
514          pitch_index = PLC_PITCH_LAG_MAX-pitch_index;
515          st->last_pitch_index = pitch_index;
516       } else {
517          pitch_index = st->last_pitch_index;
518          fade = QCONST16(.8f,15);
519       }
520
521       ALLOC(etmp, overlap, opus_val32);
522       ALLOC(exc, MAX_PERIOD, opus_val16);
523       window = mode->window;
524       c=0; do {
525          opus_val16 decay;
526          opus_val16 attenuation;
527          opus_val32 S1=0;
528          celt_sig *buf;
529          int extrapolation_offset;
530          int extrapolation_len;
531          int exc_length;
532          int j;
533
534          buf = decode_mem[c];
535          for (i=0;i<MAX_PERIOD;i++) {
536             exc[i] = ROUND16(buf[DECODE_BUFFER_SIZE-MAX_PERIOD+i], SIG_SHIFT);
537          }
538
539          if (loss_count == 0)
540          {
541             opus_val32 ac[LPC_ORDER+1];
542             /* Compute LPC coefficients for the last MAX_PERIOD samples before
543                the first loss so we can work in the excitation-filter domain. */
544             _celt_autocorr(exc, ac, window, overlap,
545                    LPC_ORDER, MAX_PERIOD, st->arch);
546             /* Add a noise floor of -40 dB. */
547 #ifdef FIXED_POINT
548             ac[0] += SHR32(ac[0],13);
549 #else
550             ac[0] *= 1.0001f;
551 #endif
552             /* Use lag windowing to stabilize the Levinson-Durbin recursion. */
553             for (i=1;i<=LPC_ORDER;i++)
554             {
555                /*ac[i] *= exp(-.5*(2*M_PI*.002*i)*(2*M_PI*.002*i));*/
556 #ifdef FIXED_POINT
557                ac[i] -= MULT16_32_Q15(2*i*i, ac[i]);
558 #else
559                ac[i] -= ac[i]*(0.008f*0.008f)*i*i;
560 #endif
561             }
562             _celt_lpc(lpc+c*LPC_ORDER, ac, LPC_ORDER);
563          }
564          /* We want the excitation for 2 pitch periods in order to look for a
565             decaying signal, but we can't get more than MAX_PERIOD. */
566          exc_length = IMIN(2*pitch_index, MAX_PERIOD);
567          /* Initialize the LPC history with the samples just before the start
568             of the region for which we're computing the excitation. */
569          {
570             opus_val16 lpc_mem[LPC_ORDER];
571             for (i=0;i<LPC_ORDER;i++)
572             {
573                lpc_mem[i] =
574                      ROUND16(buf[DECODE_BUFFER_SIZE-exc_length-1-i], SIG_SHIFT);
575             }
576             /* Compute the excitation for exc_length samples before the loss. */
577             celt_fir(exc+MAX_PERIOD-exc_length, lpc+c*LPC_ORDER,
578                   exc+MAX_PERIOD-exc_length, exc_length, LPC_ORDER, lpc_mem);
579          }
580
581          /* Check if the waveform is decaying, and if so how fast.
582             We do this to avoid adding energy when concealing in a segment
583             with decaying energy. */
584          {
585             opus_val32 E1=1, E2=1;
586             int decay_length;
587 #ifdef FIXED_POINT
588             int shift = IMAX(0,2*celt_zlog2(celt_maxabs16(&exc[MAX_PERIOD-exc_length], exc_length))-20);
589 #endif
590             decay_length = exc_length>>1;
591             for (i=0;i<decay_length;i++)
592             {
593                opus_val16 e;
594                e = exc[MAX_PERIOD-decay_length+i];
595                E1 += SHR32(MULT16_16(e, e), shift);
596                e = exc[MAX_PERIOD-2*decay_length+i];
597                E2 += SHR32(MULT16_16(e, e), shift);
598             }
599             E1 = MIN32(E1, E2);
600             decay = celt_sqrt(frac_div32(SHR32(E1, 1), E2));
601          }
602
603          /* Move the decoder memory one frame to the left to give us room to
604             add the data for the new frame. We ignore the overlap that extends
605             past the end of the buffer, because we aren't going to use it. */
606          OPUS_MOVE(buf, buf+N, DECODE_BUFFER_SIZE-N);
607
608          /* Extrapolate from the end of the excitation with a period of
609             "pitch_index", scaling down each period by an additional factor of
610             "decay". */
611          extrapolation_offset = MAX_PERIOD-pitch_index;
612          /* We need to extrapolate enough samples to cover a complete MDCT
613             window (including overlap/2 samples on both sides). */
614          extrapolation_len = N+overlap;
615          /* We also apply fading if this is not the first loss. */
616          attenuation = MULT16_16_Q15(fade, decay);
617          for (i=j=0;i<extrapolation_len;i++,j++)
618          {
619             opus_val16 tmp;
620             if (j >= pitch_index) {
621                j -= pitch_index;
622                attenuation = MULT16_16_Q15(attenuation, decay);
623             }
624             buf[DECODE_BUFFER_SIZE-N+i] =
625                   SHL32(EXTEND32(MULT16_16_Q15(attenuation,
626                         exc[extrapolation_offset+j])), SIG_SHIFT);
627             /* Compute the energy of the previously decoded signal whose
628                excitation we're copying. */
629             tmp = ROUND16(
630                   buf[DECODE_BUFFER_SIZE-MAX_PERIOD-N+extrapolation_offset+j],
631                   SIG_SHIFT);
632             S1 += SHR32(MULT16_16(tmp, tmp), 8);
633          }
634
635          {
636             opus_val16 lpc_mem[LPC_ORDER];
637             /* Copy the last decoded samples (prior to the overlap region) to
638                synthesis filter memory so we can have a continuous signal. */
639             for (i=0;i<LPC_ORDER;i++)
640                lpc_mem[i] = ROUND16(buf[DECODE_BUFFER_SIZE-N-1-i], SIG_SHIFT);
641             /* Apply the synthesis filter to convert the excitation back into
642                the signal domain. */
643             celt_iir(buf+DECODE_BUFFER_SIZE-N, lpc+c*LPC_ORDER,
644                   buf+DECODE_BUFFER_SIZE-N, extrapolation_len, LPC_ORDER,
645                   lpc_mem);
646          }
647
648          /* Check if the synthesis energy is higher than expected, which can
649             happen with the signal changes during our window. If so,
650             attenuate. */
651          {
652             opus_val32 S2=0;
653             for (i=0;i<extrapolation_len;i++)
654             {
655                opus_val16 tmp = ROUND16(buf[DECODE_BUFFER_SIZE-N+i], SIG_SHIFT);
656                S2 += SHR32(MULT16_16(tmp, tmp), 8);
657             }
658             /* This checks for an "explosion" in the synthesis. */
659 #ifdef FIXED_POINT
660             if (!(S1 > SHR32(S2,2)))
661 #else
662             /* The float test is written this way to catch NaNs in the output
663                of the IIR filter at the same time. */
664             if (!(S1 > 0.2f*S2))
665 #endif
666             {
667                for (i=0;i<extrapolation_len;i++)
668                   buf[DECODE_BUFFER_SIZE-N+i] = 0;
669             } else if (S1 < S2)
670             {
671                opus_val16 ratio = celt_sqrt(frac_div32(SHR32(S1,1)+1,S2+1));
672                for (i=0;i<overlap;i++)
673                {
674                   opus_val16 tmp_g = Q15ONE
675                         - MULT16_16_Q15(window[i], Q15ONE-ratio);
676                   buf[DECODE_BUFFER_SIZE-N+i] =
677                         MULT16_32_Q15(tmp_g, buf[DECODE_BUFFER_SIZE-N+i]);
678                }
679                for (i=overlap;i<extrapolation_len;i++)
680                {
681                   buf[DECODE_BUFFER_SIZE-N+i] =
682                         MULT16_32_Q15(ratio, buf[DECODE_BUFFER_SIZE-N+i]);
683                }
684             }
685          }
686
687          /* Apply the pre-filter to the MDCT overlap for the next frame because
688             the post-filter will be re-applied in the decoder after the MDCT
689             overlap. */
690          comb_filter(etmp, buf+DECODE_BUFFER_SIZE,
691               st->postfilter_period, st->postfilter_period, overlap,
692               -st->postfilter_gain, -st->postfilter_gain,
693               st->postfilter_tapset, st->postfilter_tapset, NULL, 0);
694
695          /* Simulate TDAC on the concealed audio so that it blends with the
696             MDCT of the next frame. */
697          for (i=0;i<overlap/2;i++)
698          {
699             buf[DECODE_BUFFER_SIZE+i] =
700                MULT16_32_Q15(window[i], etmp[overlap-1-i])
701                + MULT16_32_Q15(window[overlap-i-1], etmp[i]);
702          }
703       } while (++c<C);
704    }
705
706    st->loss_count = loss_count+1;
707
708    RESTORE_STACK;
709 }
710
711 int celt_decode_with_ec(CELTDecoder * OPUS_RESTRICT st, const unsigned char *data,
712       int len, opus_val16 * OPUS_RESTRICT pcm, int frame_size, ec_dec *dec, int accum)
713 {
714    int c, i, N;
715    int spread_decision;
716    opus_int32 bits;
717    ec_dec _dec;
718 #ifdef SMALL_FOOTPRINT
719    celt_norm *X;
720 #else
721    VARDECL(celt_norm, X);
722 #endif
723    VARDECL(int, fine_quant);
724    VARDECL(int, pulses);
725    VARDECL(int, cap);
726    VARDECL(int, offsets);
727    VARDECL(int, fine_priority);
728    VARDECL(int, tf_res);
729    VARDECL(unsigned char, collapse_masks);
730    celt_sig *decode_mem[2];
731    celt_sig *out_syn[2];
732    opus_val16 *lpc;
733    opus_val16 *oldBandE, *oldLogE, *oldLogE2, *backgroundLogE;
734
735    int shortBlocks;
736    int isTransient;
737    int intra_ener;
738    const int CC = st->channels;
739    int LM, M;
740    int start;
741    int end;
742    int effEnd;
743    int codedBands;
744    int alloc_trim;
745    int postfilter_pitch;
746    opus_val16 postfilter_gain;
747    int intensity=0;
748    int dual_stereo=0;
749    opus_int32 total_bits;
750    opus_int32 balance;
751    opus_int32 tell;
752    int dynalloc_logp;
753    int postfilter_tapset;
754    int anti_collapse_rsv;
755    int anti_collapse_on=0;
756    int silence;
757    int C = st->stream_channels;
758    const OpusCustomMode *mode;
759    int nbEBands;
760    int overlap;
761    const opus_int16 *eBands;
762    ALLOC_STACK;
763
764    mode = st->mode;
765    nbEBands = mode->nbEBands;
766    overlap = mode->overlap;
767    eBands = mode->eBands;
768    start = st->start;
769    end = st->end;
770    frame_size *= st->downsample;
771
772    lpc = (opus_val16*)(st->_decode_mem+(DECODE_BUFFER_SIZE+overlap)*CC);
773    oldBandE = lpc+CC*LPC_ORDER;
774    oldLogE = oldBandE + 2*nbEBands;
775    oldLogE2 = oldLogE + 2*nbEBands;
776    backgroundLogE = oldLogE2  + 2*nbEBands;
777
778 #ifdef CUSTOM_MODES
779    if (st->signalling && data!=NULL)
780    {
781       int data0=data[0];
782       /* Convert "standard mode" to Opus header */
783       if (mode->Fs==48000 && mode->shortMdctSize==120)
784       {
785          data0 = fromOpus(data0);
786          if (data0<0)
787             return OPUS_INVALID_PACKET;
788       }
789       st->end = end = IMAX(1, mode->effEBands-2*(data0>>5));
790       LM = (data0>>3)&0x3;
791       C = 1 + ((data0>>2)&0x1);
792       data++;
793       len--;
794       if (LM>mode->maxLM)
795          return OPUS_INVALID_PACKET;
796       if (frame_size < mode->shortMdctSize<<LM)
797          return OPUS_BUFFER_TOO_SMALL;
798       else
799          frame_size = mode->shortMdctSize<<LM;
800    } else {
801 #else
802    {
803 #endif
804       for (LM=0;LM<=mode->maxLM;LM++)
805          if (mode->shortMdctSize<<LM==frame_size)
806             break;
807       if (LM>mode->maxLM)
808          return OPUS_BAD_ARG;
809    }
810    M=1<<LM;
811
812    if (len<0 || len>1275 || pcm==NULL)
813       return OPUS_BAD_ARG;
814
815    N = M*mode->shortMdctSize;
816    c=0; do {
817       decode_mem[c] = st->_decode_mem + c*(DECODE_BUFFER_SIZE+overlap);
818       out_syn[c] = decode_mem[c]+DECODE_BUFFER_SIZE-N;
819    } while (++c<CC);
820
821    effEnd = end;
822    if (effEnd > mode->effEBands)
823       effEnd = mode->effEBands;
824
825    if (data == NULL || len<=1)
826    {
827       celt_decode_lost(st, N, LM);
828       deemphasis(out_syn, pcm, N, CC, st->downsample, mode->preemph, st->preemph_memD, accum);
829       RESTORE_STACK;
830       return frame_size/st->downsample;
831    }
832
833    if (dec == NULL)
834    {
835       ec_dec_init(&_dec,(unsigned char*)data,len);
836       dec = &_dec;
837    }
838
839    if (C==1)
840    {
841       for (i=0;i<nbEBands;i++)
842          oldBandE[i]=MAX16(oldBandE[i],oldBandE[nbEBands+i]);
843    }
844
845    total_bits = len*8;
846    tell = ec_tell(dec);
847
848    if (tell >= total_bits)
849       silence = 1;
850    else if (tell==1)
851       silence = ec_dec_bit_logp(dec, 15);
852    else
853       silence = 0;
854    if (silence)
855    {
856       /* Pretend we've read all the remaining bits */
857       tell = len*8;
858       dec->nbits_total+=tell-ec_tell(dec);
859    }
860
861    postfilter_gain = 0;
862    postfilter_pitch = 0;
863    postfilter_tapset = 0;
864    if (start==0 && tell+16 <= total_bits)
865    {
866       if(ec_dec_bit_logp(dec, 1))
867       {
868          int qg, octave;
869          octave = ec_dec_uint(dec, 6);
870          postfilter_pitch = (16<<octave)+ec_dec_bits(dec, 4+octave)-1;
871          qg = ec_dec_bits(dec, 3);
872          if (ec_tell(dec)+2<=total_bits)
873             postfilter_tapset = ec_dec_icdf(dec, tapset_icdf, 2);
874          postfilter_gain = QCONST16(.09375f,15)*(qg+1);
875       }
876       tell = ec_tell(dec);
877    }
878
879    if (LM > 0 && tell+3 <= total_bits)
880    {
881       isTransient = ec_dec_bit_logp(dec, 3);
882       tell = ec_tell(dec);
883    }
884    else
885       isTransient = 0;
886
887    if (isTransient)
888       shortBlocks = M;
889    else
890       shortBlocks = 0;
891
892    /* Decode the global flags (first symbols in the stream) */
893    intra_ener = tell+3<=total_bits ? ec_dec_bit_logp(dec, 3) : 0;
894    /* Get band energies */
895    unquant_coarse_energy(mode, start, end, oldBandE,
896          intra_ener, dec, C, LM);
897
898    ALLOC(tf_res, nbEBands, int);
899    tf_decode(start, end, isTransient, tf_res, LM, dec);
900
901    tell = ec_tell(dec);
902    spread_decision = SPREAD_NORMAL;
903    if (tell+4 <= total_bits)
904       spread_decision = ec_dec_icdf(dec, spread_icdf, 5);
905
906    ALLOC(cap, nbEBands, int);
907
908    init_caps(mode,cap,LM,C);
909
910    ALLOC(offsets, nbEBands, int);
911
912    dynalloc_logp = 6;
913    total_bits<<=BITRES;
914    tell = ec_tell_frac(dec);
915    for (i=start;i<end;i++)
916    {
917       int width, quanta;
918       int dynalloc_loop_logp;
919       int boost;
920       width = C*(eBands[i+1]-eBands[i])<<LM;
921       /* quanta is 6 bits, but no more than 1 bit/sample
922          and no less than 1/8 bit/sample */
923       quanta = IMIN(width<<BITRES, IMAX(6<<BITRES, width));
924       dynalloc_loop_logp = dynalloc_logp;
925       boost = 0;
926       while (tell+(dynalloc_loop_logp<<BITRES) < total_bits && boost < cap[i])
927       {
928          int flag;
929          flag = ec_dec_bit_logp(dec, dynalloc_loop_logp);
930          tell = ec_tell_frac(dec);
931          if (!flag)
932             break;
933          boost += quanta;
934          total_bits -= quanta;
935          dynalloc_loop_logp = 1;
936       }
937       offsets[i] = boost;
938       /* Making dynalloc more likely */
939       if (boost>0)
940          dynalloc_logp = IMAX(2, dynalloc_logp-1);
941    }
942
943    ALLOC(fine_quant, nbEBands, int);
944    alloc_trim = tell+(6<<BITRES) <= total_bits ?
945          ec_dec_icdf(dec, trim_icdf, 7) : 5;
946
947    bits = (((opus_int32)len*8)<<BITRES) - ec_tell_frac(dec) - 1;
948    anti_collapse_rsv = isTransient&&LM>=2&&bits>=((LM+2)<<BITRES) ? (1<<BITRES) : 0;
949    bits -= anti_collapse_rsv;
950
951    ALLOC(pulses, nbEBands, int);
952    ALLOC(fine_priority, nbEBands, int);
953
954    codedBands = compute_allocation(mode, start, end, offsets, cap,
955          alloc_trim, &intensity, &dual_stereo, bits, &balance, pulses,
956          fine_quant, fine_priority, C, LM, dec, 0, 0, 0);
957
958    unquant_fine_energy(mode, start, end, oldBandE, fine_quant, dec, C);
959
960    c=0; do {
961       OPUS_MOVE(decode_mem[c], decode_mem[c]+N, DECODE_BUFFER_SIZE-N+overlap/2);
962    } while (++c<CC);
963
964    /* Decode fixed codebook */
965    ALLOC(collapse_masks, C*nbEBands, unsigned char);
966
967 #ifdef SMALL_FOOTPRINT
968    /* This is an ugly hack that breaks aliasing rules and would be easily broken,
969       but it saves almost 4kB of stack. */
970    X = (celt_norm*)(out_syn[CC-1]+overlap/2);
971 #else
972    ALLOC(X, C*N, celt_norm);   /**< Interleaved normalised MDCTs */
973 #endif
974
975    quant_all_bands(0, mode, start, end, X, C==2 ? X+N : NULL, collapse_masks,
976          NULL, pulses, shortBlocks, spread_decision, dual_stereo, intensity, tf_res,
977          len*(8<<BITRES)-anti_collapse_rsv, balance, dec, LM, codedBands, &st->rng);
978
979    if (anti_collapse_rsv > 0)
980    {
981       anti_collapse_on = ec_dec_bits(dec, 1);
982    }
983
984    unquant_energy_finalise(mode, start, end, oldBandE,
985          fine_quant, fine_priority, len*8-ec_tell(dec), dec, C);
986
987    if (anti_collapse_on)
988       anti_collapse(mode, X, collapse_masks, LM, C, N,
989             start, end, oldBandE, oldLogE, oldLogE2, pulses, st->rng);
990
991    if (silence)
992    {
993       for (i=0;i<C*nbEBands;i++)
994          oldBandE[i] = -QCONST16(28.f,DB_SHIFT);
995    }
996
997    celt_synthesis(mode, X, out_syn, oldBandE, start, effEnd, C, CC, isTransient, LM, st->downsample, silence);
998
999    c=0; do {
1000       st->postfilter_period=IMAX(st->postfilter_period, COMBFILTER_MINPERIOD);
1001       st->postfilter_period_old=IMAX(st->postfilter_period_old, COMBFILTER_MINPERIOD);
1002       comb_filter(out_syn[c], out_syn[c], st->postfilter_period_old, st->postfilter_period, mode->shortMdctSize,
1003             st->postfilter_gain_old, st->postfilter_gain, st->postfilter_tapset_old, st->postfilter_tapset,
1004             mode->window, overlap);
1005       if (LM!=0)
1006          comb_filter(out_syn[c]+mode->shortMdctSize, out_syn[c]+mode->shortMdctSize, st->postfilter_period, postfilter_pitch, N-mode->shortMdctSize,
1007                st->postfilter_gain, postfilter_gain, st->postfilter_tapset, postfilter_tapset,
1008                mode->window, overlap);
1009
1010    } while (++c<CC);
1011    st->postfilter_period_old = st->postfilter_period;
1012    st->postfilter_gain_old = st->postfilter_gain;
1013    st->postfilter_tapset_old = st->postfilter_tapset;
1014    st->postfilter_period = postfilter_pitch;
1015    st->postfilter_gain = postfilter_gain;
1016    st->postfilter_tapset = postfilter_tapset;
1017    if (LM!=0)
1018    {
1019       st->postfilter_period_old = st->postfilter_period;
1020       st->postfilter_gain_old = st->postfilter_gain;
1021       st->postfilter_tapset_old = st->postfilter_tapset;
1022    }
1023
1024    if (C==1)
1025       OPUS_COPY(&oldBandE[nbEBands], oldBandE, nbEBands);
1026
1027    /* In case start or end were to change */
1028    if (!isTransient)
1029    {
1030       OPUS_COPY(oldLogE2, oldLogE, 2*nbEBands);
1031       OPUS_COPY(oldLogE, oldBandE, 2*nbEBands);
1032       for (i=0;i<2*nbEBands;i++)
1033          backgroundLogE[i] = MIN16(backgroundLogE[i] + M*QCONST16(0.001f,DB_SHIFT), oldBandE[i]);
1034    } else {
1035       for (i=0;i<2*nbEBands;i++)
1036          oldLogE[i] = MIN16(oldLogE[i], oldBandE[i]);
1037    }
1038    c=0; do
1039    {
1040       for (i=0;i<start;i++)
1041       {
1042          oldBandE[c*nbEBands+i]=0;
1043          oldLogE[c*nbEBands+i]=oldLogE2[c*nbEBands+i]=-QCONST16(28.f,DB_SHIFT);
1044       }
1045       for (i=end;i<nbEBands;i++)
1046       {
1047          oldBandE[c*nbEBands+i]=0;
1048          oldLogE[c*nbEBands+i]=oldLogE2[c*nbEBands+i]=-QCONST16(28.f,DB_SHIFT);
1049       }
1050    } while (++c<2);
1051    st->rng = dec->rng;
1052
1053    deemphasis(out_syn, pcm, N, CC, st->downsample, mode->preemph, st->preemph_memD, accum);
1054    st->loss_count = 0;
1055    RESTORE_STACK;
1056    if (ec_tell(dec) > 8*len)
1057       return OPUS_INTERNAL_ERROR;
1058    if(ec_get_error(dec))
1059       st->error = 1;
1060    return frame_size/st->downsample;
1061 }
1062
1063
1064 #ifdef CUSTOM_MODES
1065
1066 #ifdef FIXED_POINT
1067 int opus_custom_decode(CELTDecoder * OPUS_RESTRICT st, const unsigned char *data, int len, opus_int16 * OPUS_RESTRICT pcm, int frame_size)
1068 {
1069    return celt_decode_with_ec(st, data, len, pcm, frame_size, NULL, 0);
1070 }
1071
1072 #ifndef DISABLE_FLOAT_API
1073 int opus_custom_decode_float(CELTDecoder * OPUS_RESTRICT st, const unsigned char *data, int len, float * OPUS_RESTRICT pcm, int frame_size)
1074 {
1075    int j, ret, C, N;
1076    VARDECL(opus_int16, out);
1077    ALLOC_STACK;
1078
1079    if (pcm==NULL)
1080       return OPUS_BAD_ARG;
1081
1082    C = st->channels;
1083    N = frame_size;
1084
1085    ALLOC(out, C*N, opus_int16);
1086    ret=celt_decode_with_ec(st, data, len, out, frame_size, NULL, 0);
1087    if (ret>0)
1088       for (j=0;j<C*ret;j++)
1089          pcm[j]=out[j]*(1.f/32768.f);
1090
1091    RESTORE_STACK;
1092    return ret;
1093 }
1094 #endif /* DISABLE_FLOAT_API */
1095
1096 #else
1097
1098 int opus_custom_decode_float(CELTDecoder * OPUS_RESTRICT st, const unsigned char *data, int len, float * OPUS_RESTRICT pcm, int frame_size)
1099 {
1100    return celt_decode_with_ec(st, data, len, pcm, frame_size, NULL, 0);
1101 }
1102
1103 int opus_custom_decode(CELTDecoder * OPUS_RESTRICT st, const unsigned char *data, int len, opus_int16 * OPUS_RESTRICT pcm, int frame_size)
1104 {
1105    int j, ret, C, N;
1106    VARDECL(celt_sig, out);
1107    ALLOC_STACK;
1108
1109    if (pcm==NULL)
1110       return OPUS_BAD_ARG;
1111
1112    C = st->channels;
1113    N = frame_size;
1114    ALLOC(out, C*N, celt_sig);
1115
1116    ret=celt_decode_with_ec(st, data, len, out, frame_size, NULL, 0);
1117
1118    if (ret>0)
1119       for (j=0;j<C*ret;j++)
1120          pcm[j] = FLOAT2INT16 (out[j]);
1121
1122    RESTORE_STACK;
1123    return ret;
1124 }
1125
1126 #endif
1127 #endif /* CUSTOM_MODES */
1128
1129 int opus_custom_decoder_ctl(CELTDecoder * OPUS_RESTRICT st, int request, ...)
1130 {
1131    va_list ap;
1132
1133    va_start(ap, request);
1134    switch (request)
1135    {
1136       case CELT_SET_START_BAND_REQUEST:
1137       {
1138          opus_int32 value = va_arg(ap, opus_int32);
1139          if (value<0 || value>=st->mode->nbEBands)
1140             goto bad_arg;
1141          st->start = value;
1142       }
1143       break;
1144       case CELT_SET_END_BAND_REQUEST:
1145       {
1146          opus_int32 value = va_arg(ap, opus_int32);
1147          if (value<1 || value>st->mode->nbEBands)
1148             goto bad_arg;
1149          st->end = value;
1150       }
1151       break;
1152       case CELT_SET_CHANNELS_REQUEST:
1153       {
1154          opus_int32 value = va_arg(ap, opus_int32);
1155          if (value<1 || value>2)
1156             goto bad_arg;
1157          st->stream_channels = value;
1158       }
1159       break;
1160       case CELT_GET_AND_CLEAR_ERROR_REQUEST:
1161       {
1162          opus_int32 *value = va_arg(ap, opus_int32*);
1163          if (value==NULL)
1164             goto bad_arg;
1165          *value=st->error;
1166          st->error = 0;
1167       }
1168       break;
1169       case OPUS_GET_LOOKAHEAD_REQUEST:
1170       {
1171          opus_int32 *value = va_arg(ap, opus_int32*);
1172          if (value==NULL)
1173             goto bad_arg;
1174          *value = st->overlap/st->downsample;
1175       }
1176       break;
1177       case OPUS_RESET_STATE:
1178       {
1179          int i;
1180          opus_val16 *lpc, *oldBandE, *oldLogE, *oldLogE2;
1181          lpc = (opus_val16*)(st->_decode_mem+(DECODE_BUFFER_SIZE+st->overlap)*st->channels);
1182          oldBandE = lpc+st->channels*LPC_ORDER;
1183          oldLogE = oldBandE + 2*st->mode->nbEBands;
1184          oldLogE2 = oldLogE + 2*st->mode->nbEBands;
1185          OPUS_CLEAR((char*)&st->DECODER_RESET_START,
1186                opus_custom_decoder_get_size(st->mode, st->channels)-
1187                ((char*)&st->DECODER_RESET_START - (char*)st));
1188          for (i=0;i<2*st->mode->nbEBands;i++)
1189             oldLogE[i]=oldLogE2[i]=-QCONST16(28.f,DB_SHIFT);
1190       }
1191       break;
1192       case OPUS_GET_PITCH_REQUEST:
1193       {
1194          opus_int32 *value = va_arg(ap, opus_int32*);
1195          if (value==NULL)
1196             goto bad_arg;
1197          *value = st->postfilter_period;
1198       }
1199       break;
1200       case CELT_GET_MODE_REQUEST:
1201       {
1202          const CELTMode ** value = va_arg(ap, const CELTMode**);
1203          if (value==0)
1204             goto bad_arg;
1205          *value=st->mode;
1206       }
1207       break;
1208       case CELT_SET_SIGNALLING_REQUEST:
1209       {
1210          opus_int32 value = va_arg(ap, opus_int32);
1211          st->signalling = value;
1212       }
1213       break;
1214       case OPUS_GET_FINAL_RANGE_REQUEST:
1215       {
1216          opus_uint32 * value = va_arg(ap, opus_uint32 *);
1217          if (value==0)
1218             goto bad_arg;
1219          *value=st->rng;
1220       }
1221       break;
1222       default:
1223          goto bad_request;
1224    }
1225    va_end(ap);
1226    return OPUS_OK;
1227 bad_arg:
1228    va_end(ap);
1229    return OPUS_BAD_ARG;
1230 bad_request:
1231       va_end(ap);
1232   return OPUS_UNIMPLEMENTED;
1233 }