Adds SMALL_FOOTPRINT hack to the CELT PLC too
[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 #ifdef SMALL_FOOTPRINT
446       celt_norm *X;
447 #else
448       VARDECL(celt_norm, X);
449 #endif
450       opus_uint32 seed;
451       opus_val16 *plcLogE;
452       int end;
453       int effEnd;
454
455       end = st->end;
456       effEnd = IMAX(start, IMIN(end, mode->effEBands));
457
458 #ifdef SMALL_FOOTPRINT
459       /* This is an ugly hack that breaks aliasing rules and would be easily broken,
460          but it saves almost 4kB of stack. */
461       X = (celt_norm*)(out_syn[C-1]+overlap/2);
462 #else
463       ALLOC(X, C*N, celt_norm);   /**< Interleaved normalised MDCTs */
464 #endif
465
466       if (loss_count >= 5)
467          plcLogE = backgroundLogE;
468       else {
469          /* Energy decay */
470          opus_val16 decay = loss_count==0 ?
471                QCONST16(1.5f, DB_SHIFT) : QCONST16(.5f, DB_SHIFT);
472          c=0; do
473          {
474             for (i=start;i<end;i++)
475                oldBandE[c*nbEBands+i] -= decay;
476          } while (++c<C);
477          plcLogE = oldBandE;
478       }
479       seed = st->rng;
480       for (c=0;c<C;c++)
481       {
482          for (i=start;i<effEnd;i++)
483          {
484             int j;
485             int boffs;
486             int blen;
487             boffs = N*c+(eBands[i]<<LM);
488             blen = (eBands[i+1]-eBands[i])<<LM;
489             for (j=0;j<blen;j++)
490             {
491                seed = celt_lcg_rand(seed);
492                X[boffs+j] = (celt_norm)((opus_int32)seed>>20);
493             }
494             renormalise_vector(X+boffs, blen, Q15ONE);
495          }
496       }
497       st->rng = seed;
498
499       c=0; do {
500          OPUS_MOVE(decode_mem[c], decode_mem[c]+N,
501                DECODE_BUFFER_SIZE-N+(overlap>>1));
502       } while (++c<C);
503
504       celt_synthesis(mode, X, out_syn, plcLogE, start, effEnd, C, C, 0, LM, st->downsample, 0);
505    } else {
506       /* Pitch-based PLC */
507       const opus_val16 *window;
508       opus_val16 fade = Q15ONE;
509       int pitch_index;
510       VARDECL(opus_val32, etmp);
511       VARDECL(opus_val16, exc);
512
513       if (loss_count == 0)
514       {
515          VARDECL( opus_val16, lp_pitch_buf );
516          ALLOC( lp_pitch_buf, DECODE_BUFFER_SIZE>>1, opus_val16 );
517          pitch_downsample(decode_mem, lp_pitch_buf,
518                DECODE_BUFFER_SIZE, C, st->arch);
519          pitch_search(lp_pitch_buf+(PLC_PITCH_LAG_MAX>>1), lp_pitch_buf,
520                DECODE_BUFFER_SIZE-PLC_PITCH_LAG_MAX,
521                PLC_PITCH_LAG_MAX-PLC_PITCH_LAG_MIN, &pitch_index, st->arch);
522          pitch_index = PLC_PITCH_LAG_MAX-pitch_index;
523          st->last_pitch_index = pitch_index;
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 SMALL_FOOTPRINT
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 SMALL_FOOTPRINT
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 }