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