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