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