Moves VBR calculations to a separate function.
[opus.git] / celt / celt_encoder.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_ENCODER_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 /** Encoder state
56  @brief Encoder state
57  */
58 struct OpusCustomEncoder {
59    const OpusCustomMode *mode;     /**< Mode used by the encoder */
60    int overlap;
61    int channels;
62    int stream_channels;
63
64    int force_intra;
65    int clip;
66    int disable_pf;
67    int complexity;
68    int upsample;
69    int start, end;
70
71    opus_int32 bitrate;
72    int vbr;
73    int signalling;
74    int constrained_vbr;      /* If zero, VBR can do whatever it likes with the rate */
75    int loss_rate;
76    int lsb_depth;
77    int variable_duration;
78    int lfe;
79    int arch;
80
81    /* Everything beyond this point gets cleared on a reset */
82 #define ENCODER_RESET_START rng
83
84    opus_uint32 rng;
85    int spread_decision;
86    opus_val32 delayedIntra;
87    int tonal_average;
88    int lastCodedBands;
89    int hf_average;
90    int tapset_decision;
91
92    int prefilter_period;
93    opus_val16 prefilter_gain;
94    int prefilter_tapset;
95 #ifdef RESYNTH
96    int prefilter_period_old;
97    opus_val16 prefilter_gain_old;
98    int prefilter_tapset_old;
99 #endif
100    int consec_transient;
101    AnalysisInfo analysis;
102
103    opus_val32 preemph_memE[2];
104    opus_val32 preemph_memD[2];
105
106    /* VBR-related parameters */
107    opus_int32 vbr_reservoir;
108    opus_int32 vbr_drift;
109    opus_int32 vbr_offset;
110    opus_int32 vbr_count;
111    opus_val32 overlap_max;
112    opus_val16 stereo_saving;
113    int intensity;
114    opus_val16 *energy_save;
115    opus_val16 *energy_mask;
116
117 #ifdef RESYNTH
118    /* +MAX_PERIOD/2 to make space for overlap */
119    celt_sig syn_mem[2][2*MAX_PERIOD+MAX_PERIOD/2];
120 #endif
121
122    celt_sig in_mem[1]; /* Size = channels*mode->overlap */
123    /* celt_sig prefilter_mem[],  Size = channels*COMBFILTER_MAXPERIOD */
124    /* opus_val16 oldBandE[],     Size = channels*mode->nbEBands */
125    /* opus_val16 oldLogE[],      Size = channels*mode->nbEBands */
126    /* opus_val16 oldLogE2[],     Size = channels*mode->nbEBands */
127 };
128
129 int celt_encoder_get_size(int channels)
130 {
131    CELTMode *mode = opus_custom_mode_create(48000, 960, NULL);
132    return opus_custom_encoder_get_size(mode, channels);
133 }
134
135 OPUS_CUSTOM_NOSTATIC int opus_custom_encoder_get_size(const CELTMode *mode, int channels)
136 {
137    int size = sizeof(struct CELTEncoder)
138          + (channels*mode->overlap-1)*sizeof(celt_sig)    /* celt_sig in_mem[channels*mode->overlap]; */
139          + channels*COMBFILTER_MAXPERIOD*sizeof(celt_sig) /* celt_sig prefilter_mem[channels*COMBFILTER_MAXPERIOD]; */
140          + 3*channels*mode->nbEBands*sizeof(opus_val16);  /* opus_val16 oldBandE[channels*mode->nbEBands]; */
141                                                           /* opus_val16 oldLogE[channels*mode->nbEBands]; */
142                                                           /* opus_val16 oldLogE2[channels*mode->nbEBands]; */
143    return size;
144 }
145
146 #ifdef CUSTOM_MODES
147 CELTEncoder *opus_custom_encoder_create(const CELTMode *mode, int channels, int *error)
148 {
149    int ret;
150    CELTEncoder *st = (CELTEncoder *)opus_alloc(opus_custom_encoder_get_size(mode, channels));
151    /* init will handle the NULL case */
152    ret = opus_custom_encoder_init(st, mode, channels);
153    if (ret != OPUS_OK)
154    {
155       opus_custom_encoder_destroy(st);
156       st = NULL;
157    }
158    if (error)
159       *error = ret;
160    return st;
161 }
162 #endif /* CUSTOM_MODES */
163
164 int celt_encoder_init(CELTEncoder *st, opus_int32 sampling_rate, int channels)
165 {
166    int ret;
167    ret = opus_custom_encoder_init(st, opus_custom_mode_create(48000, 960, NULL), channels);
168    if (ret != OPUS_OK)
169       return ret;
170    st->upsample = resampling_factor(sampling_rate);
171    return OPUS_OK;
172 }
173
174 OPUS_CUSTOM_NOSTATIC int opus_custom_encoder_init(CELTEncoder *st, const CELTMode *mode, int channels)
175 {
176    if (channels < 0 || channels > 2)
177       return OPUS_BAD_ARG;
178
179    if (st==NULL || mode==NULL)
180       return OPUS_ALLOC_FAIL;
181
182    OPUS_CLEAR((char*)st, opus_custom_encoder_get_size(mode, channels));
183
184    st->mode = mode;
185    st->overlap = mode->overlap;
186    st->stream_channels = st->channels = channels;
187
188    st->upsample = 1;
189    st->start = 0;
190    st->end = st->mode->effEBands;
191    st->signalling = 1;
192
193    st->arch = opus_select_arch();
194
195    st->constrained_vbr = 1;
196    st->clip = 1;
197
198    st->bitrate = OPUS_BITRATE_MAX;
199    st->vbr = 0;
200    st->force_intra  = 0;
201    st->complexity = 5;
202    st->lsb_depth=24;
203
204    opus_custom_encoder_ctl(st, OPUS_RESET_STATE);
205
206    return OPUS_OK;
207 }
208
209 #ifdef CUSTOM_MODES
210 void opus_custom_encoder_destroy(CELTEncoder *st)
211 {
212    opus_free(st);
213 }
214 #endif /* CUSTOM_MODES */
215
216
217 static int transient_analysis(const opus_val32 * OPUS_RESTRICT in, int len, int C,
218                               opus_val16 *tf_estimate, int *tf_chan)
219 {
220    int i;
221    VARDECL(opus_val16, tmp);
222    opus_val32 mem0,mem1;
223    int is_transient = 0;
224    opus_int32 mask_metric = 0;
225    int c;
226    opus_val16 tf_max;
227    int len2;
228    /* Table of 6*64/x, trained on real data to minimize the average error */
229    static const unsigned char inv_table[128] = {
230          255,255,156,110, 86, 70, 59, 51, 45, 40, 37, 33, 31, 28, 26, 25,
231           23, 22, 21, 20, 19, 18, 17, 16, 16, 15, 15, 14, 13, 13, 12, 12,
232           12, 12, 11, 11, 11, 10, 10, 10,  9,  9,  9,  9,  9,  9,  8,  8,
233            8,  8,  8,  7,  7,  7,  7,  7,  7,  6,  6,  6,  6,  6,  6,  6,
234            6,  6,  6,  6,  6,  6,  6,  6,  6,  5,  5,  5,  5,  5,  5,  5,
235            5,  5,  5,  5,  5,  4,  4,  4,  4,  4,  4,  4,  4,  4,  4,  4,
236            4,  4,  4,  4,  4,  4,  4,  4,  4,  4,  4,  4,  4,  4,  3,  3,
237            3,  3,  3,  3,  3,  3,  3,  3,  3,  3,  3,  3,  3,  3,  3,  2,
238    };
239    SAVE_STACK;
240    ALLOC(tmp, len, opus_val16);
241
242    len2=len/2;
243    tf_max = 0;
244    for (c=0;c<C;c++)
245    {
246       opus_val32 mean;
247       opus_int32 unmask=0;
248       opus_val32 norm;
249       opus_val16 maxE;
250       mem0=0;
251       mem1=0;
252       /* High-pass filter: (1 - 2*z^-1 + z^-2) / (1 - z^-1 + .5*z^-2) */
253       for (i=0;i<len;i++)
254       {
255          opus_val32 x,y;
256          x = SHR32(in[i+c*len],SIG_SHIFT);
257          y = ADD32(mem0, x);
258 #ifdef FIXED_POINT
259          mem0 = mem1 + y - SHL32(x,1);
260          mem1 = x - SHR32(y,1);
261 #else
262          mem0 = mem1 + y - 2*x;
263          mem1 = x - .5f*y;
264 #endif
265          tmp[i] = EXTRACT16(SHR32(y,2));
266          /*printf("%f ", tmp[i]);*/
267       }
268       /*printf("\n");*/
269       /* First few samples are bad because we don't propagate the memory */
270       for (i=0;i<12;i++)
271          tmp[i] = 0;
272
273 #ifdef FIXED_POINT
274       /* Normalize tmp to max range */
275       {
276          int shift=0;
277          shift = 14-celt_ilog2(1+celt_maxabs16(tmp, len));
278          if (shift!=0)
279          {
280             for (i=0;i<len;i++)
281                tmp[i] = SHL16(tmp[i], shift);
282          }
283       }
284 #endif
285
286       mean=0;
287       mem0=0;
288       /* Grouping by two to reduce complexity */
289       /* Forward pass to compute the post-echo threshold*/
290       for (i=0;i<len2;i++)
291       {
292          opus_val16 x2 = PSHR32(MULT16_16(tmp[2*i],tmp[2*i]) + MULT16_16(tmp[2*i+1],tmp[2*i+1]),16);
293          mean += x2;
294 #ifdef FIXED_POINT
295          /* FIXME: Use PSHR16() instead */
296          tmp[i] = mem0 + PSHR32(x2-mem0,4);
297 #else
298          tmp[i] = mem0 + MULT16_16_P15(QCONST16(.0625f,15),x2-mem0);
299 #endif
300          mem0 = tmp[i];
301       }
302
303       mem0=0;
304       maxE=0;
305       /* Backward pass to compute the pre-echo threshold */
306       for (i=len2-1;i>=0;i--)
307       {
308 #ifdef FIXED_POINT
309          /* FIXME: Use PSHR16() instead */
310          tmp[i] = mem0 + PSHR32(tmp[i]-mem0,3);
311 #else
312          tmp[i] = mem0 + MULT16_16_P15(QCONST16(0.125f,15),tmp[i]-mem0);
313 #endif
314          mem0 = tmp[i];
315          maxE = MAX16(maxE, mem0);
316       }
317       /*for (i=0;i<len2;i++)printf("%f ", tmp[i]/mean);printf("\n");*/
318
319       /* Compute the ratio of the "frame energy" over the harmonic mean of the energy.
320          This essentially corresponds to a bitrate-normalized temporal noise-to-mask
321          ratio */
322
323       /* As a compromise with the old transient detector, frame energy is the
324          geometric mean of the energy and half the max */
325 #ifdef FIXED_POINT
326       /* Costs two sqrt() to avoid overflows */
327       mean = MULT16_16(celt_sqrt(mean), celt_sqrt(MULT16_16(maxE,len2>>1)));
328 #else
329       mean = celt_sqrt(mean * maxE*.5*len2);
330 #endif
331       /* Inverse of the mean energy in Q15+6 */
332       norm = SHL32(EXTEND32(len2),6+14)/ADD32(EPSILON,SHR32(mean,1));
333       /* Compute harmonic mean discarding the unreliable boundaries
334          The data is smooth, so we only take 1/4th of the samples */
335       unmask=0;
336       for (i=12;i<len2-5;i+=4)
337       {
338          int id;
339 #ifdef FIXED_POINT
340          id = IMAX(0,IMIN(127,MULT16_32_Q15(tmp[i],norm))); /* Do not round to nearest */
341 #else
342          id = IMAX(0,IMIN(127,(int)floor(64*norm*tmp[i]))); /* Do not round to nearest */
343 #endif
344          unmask += inv_table[id];
345       }
346       /*printf("%d\n", unmask);*/
347       /* Normalize, compensate for the 1/4th of the sample and the factor of 6 in the inverse table */
348       unmask = 64*unmask*4/(6*(len2-17));
349       if (unmask>mask_metric)
350       {
351          *tf_chan = c;
352          mask_metric = unmask;
353       }
354    }
355    is_transient = mask_metric>200;
356
357    /* Arbitrary metric for VBR boost */
358    tf_max = MAX16(0,celt_sqrt(27*mask_metric)-42);
359    /* *tf_estimate = 1 + MIN16(1, sqrt(MAX16(0, tf_max-30))/20); */
360    *tf_estimate = celt_sqrt(MAX16(0, SHL32(MULT16_16(QCONST16(0.0069,14),MIN16(163,tf_max)),14)-QCONST32(0.139,28)));
361    /*printf("%d %f\n", tf_max, mask_metric);*/
362    RESTORE_STACK;
363 #ifdef FUZZING
364    is_transient = rand()&0x1;
365 #endif
366    /*printf("%d %f %d\n", is_transient, (float)*tf_estimate, tf_max);*/
367    return is_transient;
368 }
369
370 /* Looks for sudden increases of energy to decide whether we need to patch
371    the transient decision */
372 int patch_transient_decision(opus_val16 *new, opus_val16 *old, int nbEBands,
373       int end, int C)
374 {
375    int i, c;
376    opus_val32 mean_diff=0;
377    opus_val16 spread_old[26];
378    /* Apply an aggressive (-6 dB/Bark) spreading function to the old frame to
379       avoid false detection caused by irrelevant bands */
380    if (C==1)
381    {
382       spread_old[0] = old[0];
383       for (i=1;i<end;i++)
384          spread_old[i] = MAX16(spread_old[i-1]-QCONST16(1.0f, DB_SHIFT), old[i]);
385    } else {
386       spread_old[0] = MAX16(old[0],old[nbEBands]);
387       for (i=1;i<end;i++)
388          spread_old[i] = MAX16(spread_old[i-1]-QCONST16(1.0f, DB_SHIFT),
389                                MAX16(old[i],old[i+nbEBands]));
390    }
391    for (i=end-2;i>=0;i--)
392       spread_old[i] = MAX16(spread_old[i], spread_old[i+1]-QCONST16(1.0f, DB_SHIFT));
393    /* Compute mean increase */
394    c=0; do {
395       for (i=2;i<end-1;i++)
396       {
397          opus_val16 x1, x2;
398          x1 = MAX16(0, new[i]);
399          x2 = MAX16(0, spread_old[i]);
400          mean_diff = ADD32(mean_diff, EXTEND32(MAX16(0, SUB16(x1, x2))));
401       }
402    } while (++c<C);
403    mean_diff = DIV32(mean_diff, C*(end-3));
404    /*printf("%f %f %d\n", mean_diff, max_diff, count);*/
405    return mean_diff > QCONST16(1.f, DB_SHIFT);
406 }
407
408 /** Apply window and compute the MDCT for all sub-frames and
409     all channels in a frame */
410 static void compute_mdcts(const CELTMode *mode, int shortBlocks, celt_sig * OPUS_RESTRICT in,
411                           celt_sig * OPUS_RESTRICT out, int C, int CC, int LM, int upsample)
412 {
413    const int overlap = OVERLAP(mode);
414    int N;
415    int B;
416    int shift;
417    int i, b, c;
418    if (shortBlocks)
419    {
420       B = shortBlocks;
421       N = mode->shortMdctSize;
422       shift = mode->maxLM;
423    } else {
424       B = 1;
425       N = mode->shortMdctSize<<LM;
426       shift = mode->maxLM-LM;
427    }
428    c=0; do {
429       for (b=0;b<B;b++)
430       {
431          /* Interleaving the sub-frames while doing the MDCTs */
432          clt_mdct_forward(&mode->mdct, in+c*(B*N+overlap)+b*N, &out[b+c*N*B], mode->window, overlap, shift, B);
433       }
434    } while (++c<CC);
435    if (CC==2&&C==1)
436    {
437       for (i=0;i<B*N;i++)
438          out[i] = ADD32(HALF32(out[i]), HALF32(out[B*N+i]));
439    }
440    if (upsample != 1)
441    {
442       c=0; do
443       {
444          int bound = B*N/upsample;
445          for (i=0;i<bound;i++)
446             out[c*B*N+i] *= upsample;
447          for (;i<B*N;i++)
448             out[c*B*N+i] = 0;
449       } while (++c<C);
450    }
451 }
452
453
454 static void preemphasis(const opus_val16 * OPUS_RESTRICT pcmp, celt_sig * OPUS_RESTRICT inp,
455                         int N, int CC, int upsample, const opus_val16 *coef, celt_sig *mem, int clip)
456 {
457    int i;
458    opus_val16 coef0;
459    celt_sig m;
460    int Nu;
461
462    coef0 = coef[0];
463
464
465    Nu = N/upsample;
466    if (upsample!=1)
467    {
468       for (i=0;i<N;i++)
469          inp[i] = 0;
470    }
471    for (i=0;i<Nu;i++)
472    {
473       celt_sig x;
474
475       x = SCALEIN(pcmp[CC*i]);
476 #ifndef FIXED_POINT
477       /* Replace NaNs with zeros */
478       if (!(x==x))
479          x = 0;
480 #endif
481       inp[i*upsample] = x;
482    }
483
484 #ifndef FIXED_POINT
485    if (clip)
486    {
487       /* Clip input to avoid encoding non-portable files */
488       for (i=0;i<Nu;i++)
489          inp[i*upsample] = MAX32(-65536.f, MIN32(65536.f,inp[i*upsample]));
490    }
491 #endif
492    m = *mem;
493 #ifdef CUSTOM_MODES
494    if (coef[1] != 0)
495    {
496       opus_val16 coef1 = coef[1];
497       opus_val16 coef2 = coef[2];
498       for (i=0;i<N;i++)
499       {
500          opus_val16 x, tmp;
501          x = inp[i];
502          /* Apply pre-emphasis */
503          tmp = MULT16_16(coef2, x);
504          inp[i] = tmp + m;
505          m = MULT16_32_Q15(coef1, inp[i]) - MULT16_32_Q15(coef0, tmp);
506       }
507    } else
508 #endif
509    {
510       for (i=0;i<N;i++)
511       {
512          celt_sig x;
513          x = SHL32(inp[i], SIG_SHIFT);
514          /* Apply pre-emphasis */
515          inp[i] = x + m;
516          m = - MULT16_32_Q15(coef0, x);
517       }
518    }
519    *mem = m;
520 }
521
522
523
524 static opus_val32 l1_metric(const celt_norm *tmp, int N, int LM, opus_val16 bias)
525 {
526    int i;
527    opus_val32 L1;
528    L1 = 0;
529    for (i=0;i<N;i++)
530       L1 += EXTEND32(ABS16(tmp[i]));
531    /* When in doubt, prefer good freq resolution */
532    L1 = MAC16_32_Q15(L1, LM*bias, L1);
533    return L1;
534
535 }
536
537 static int tf_analysis(const CELTMode *m, int len, int isTransient,
538       int *tf_res, int lambda, celt_norm *X, int N0, int LM,
539       int *tf_sum, opus_val16 tf_estimate, int tf_chan)
540 {
541    int i;
542    VARDECL(int, metric);
543    int cost0;
544    int cost1;
545    VARDECL(int, path0);
546    VARDECL(int, path1);
547    VARDECL(celt_norm, tmp);
548    VARDECL(celt_norm, tmp_1);
549    int sel;
550    int selcost[2];
551    int tf_select=0;
552    opus_val16 bias;
553
554    SAVE_STACK;
555    bias = MULT16_16_Q14(QCONST16(.04f,15), MAX16(-QCONST16(.25f,14), QCONST16(.5f,14)-tf_estimate));
556    /*printf("%f ", bias);*/
557
558    ALLOC(metric, len, int);
559    ALLOC(tmp, (m->eBands[len]-m->eBands[len-1])<<LM, celt_norm);
560    ALLOC(tmp_1, (m->eBands[len]-m->eBands[len-1])<<LM, celt_norm);
561    ALLOC(path0, len, int);
562    ALLOC(path1, len, int);
563
564    *tf_sum = 0;
565    for (i=0;i<len;i++)
566    {
567       int j, k, N;
568       int narrow;
569       opus_val32 L1, best_L1;
570       int best_level=0;
571       N = (m->eBands[i+1]-m->eBands[i])<<LM;
572       /* band is too narrow to be split down to LM=-1 */
573       narrow = (m->eBands[i+1]-m->eBands[i])==1;
574       for (j=0;j<N;j++)
575          tmp[j] = X[tf_chan*N0 + j+(m->eBands[i]<<LM)];
576       /* Just add the right channel if we're in stereo */
577       /*if (C==2)
578          for (j=0;j<N;j++)
579             tmp[j] = ADD16(SHR16(tmp[j], 1),SHR16(X[N0+j+(m->eBands[i]<<LM)], 1));*/
580       L1 = l1_metric(tmp, N, isTransient ? LM : 0, bias);
581       best_L1 = L1;
582       /* Check the -1 case for transients */
583       if (isTransient && !narrow)
584       {
585          for (j=0;j<N;j++)
586             tmp_1[j] = tmp[j];
587          haar1(tmp_1, N>>LM, 1<<LM);
588          L1 = l1_metric(tmp_1, N, LM+1, bias);
589          if (L1<best_L1)
590          {
591             best_L1 = L1;
592             best_level = -1;
593          }
594       }
595       /*printf ("%f ", L1);*/
596       for (k=0;k<LM+!(isTransient||narrow);k++)
597       {
598          int B;
599
600          if (isTransient)
601             B = (LM-k-1);
602          else
603             B = k+1;
604
605          haar1(tmp, N>>k, 1<<k);
606
607          L1 = l1_metric(tmp, N, B, bias);
608
609          if (L1 < best_L1)
610          {
611             best_L1 = L1;
612             best_level = k+1;
613          }
614       }
615       /*printf ("%d ", isTransient ? LM-best_level : best_level);*/
616       /* metric is in Q1 to be able to select the mid-point (-0.5) for narrower bands */
617       if (isTransient)
618          metric[i] = 2*best_level;
619       else
620          metric[i] = -2*best_level;
621       *tf_sum += (isTransient ? LM : 0) - metric[i]/2;
622       /* For bands that can't be split to -1, set the metric to the half-way point to avoid
623          biasing the decision */
624       if (narrow && (metric[i]==0 || metric[i]==-2*LM))
625          metric[i]-=1;
626       /*printf("%d ", metric[i]);*/
627    }
628    /*printf("\n");*/
629    /* Search for the optimal tf resolution, including tf_select */
630    tf_select = 0;
631    for (sel=0;sel<2;sel++)
632    {
633       cost0 = 0;
634       cost1 = isTransient ? 0 : lambda;
635       for (i=1;i<len;i++)
636       {
637          int curr0, curr1;
638          curr0 = IMIN(cost0, cost1 + lambda);
639          curr1 = IMIN(cost0 + lambda, cost1);
640          cost0 = curr0 + abs(metric[i]-2*tf_select_table[LM][4*isTransient+2*sel+0]);
641          cost1 = curr1 + abs(metric[i]-2*tf_select_table[LM][4*isTransient+2*sel+1]);
642       }
643       cost0 = IMIN(cost0, cost1);
644       selcost[sel]=cost0;
645    }
646    /* For now, we're conservative and only allow tf_select=1 for transients.
647     * If tests confirm it's useful for non-transients, we could allow it. */
648    if (selcost[1]<selcost[0] && isTransient)
649       tf_select=1;
650    cost0 = 0;
651    cost1 = isTransient ? 0 : lambda;
652    /* Viterbi forward pass */
653    for (i=1;i<len;i++)
654    {
655       int curr0, curr1;
656       int from0, from1;
657
658       from0 = cost0;
659       from1 = cost1 + lambda;
660       if (from0 < from1)
661       {
662          curr0 = from0;
663          path0[i]= 0;
664       } else {
665          curr0 = from1;
666          path0[i]= 1;
667       }
668
669       from0 = cost0 + lambda;
670       from1 = cost1;
671       if (from0 < from1)
672       {
673          curr1 = from0;
674          path1[i]= 0;
675       } else {
676          curr1 = from1;
677          path1[i]= 1;
678       }
679       cost0 = curr0 + abs(metric[i]-2*tf_select_table[LM][4*isTransient+2*tf_select+0]);
680       cost1 = curr1 + abs(metric[i]-2*tf_select_table[LM][4*isTransient+2*tf_select+1]);
681    }
682    tf_res[len-1] = cost0 < cost1 ? 0 : 1;
683    /* Viterbi backward pass to check the decisions */
684    for (i=len-2;i>=0;i--)
685    {
686       if (tf_res[i+1] == 1)
687          tf_res[i] = path1[i+1];
688       else
689          tf_res[i] = path0[i+1];
690    }
691    /*printf("%d %f\n", *tf_sum, tf_estimate);*/
692    RESTORE_STACK;
693 #ifdef FUZZING
694    tf_select = rand()&0x1;
695    tf_res[0] = rand()&0x1;
696    for (i=1;i<len;i++)
697       tf_res[i] = tf_res[i-1] ^ ((rand()&0xF) == 0);
698 #endif
699    return tf_select;
700 }
701
702 static void tf_encode(int start, int end, int isTransient, int *tf_res, int LM, int tf_select, ec_enc *enc)
703 {
704    int curr, i;
705    int tf_select_rsv;
706    int tf_changed;
707    int logp;
708    opus_uint32 budget;
709    opus_uint32 tell;
710    budget = enc->storage*8;
711    tell = ec_tell(enc);
712    logp = isTransient ? 2 : 4;
713    /* Reserve space to code the tf_select decision. */
714    tf_select_rsv = LM>0 && tell+logp+1 <= budget;
715    budget -= tf_select_rsv;
716    curr = tf_changed = 0;
717    for (i=start;i<end;i++)
718    {
719       if (tell+logp<=budget)
720       {
721          ec_enc_bit_logp(enc, tf_res[i] ^ curr, logp);
722          tell = ec_tell(enc);
723          curr = tf_res[i];
724          tf_changed |= curr;
725       }
726       else
727          tf_res[i] = curr;
728       logp = isTransient ? 4 : 5;
729    }
730    /* Only code tf_select if it would actually make a difference. */
731    if (tf_select_rsv &&
732          tf_select_table[LM][4*isTransient+0+tf_changed]!=
733          tf_select_table[LM][4*isTransient+2+tf_changed])
734       ec_enc_bit_logp(enc, tf_select, 1);
735    else
736       tf_select = 0;
737    for (i=start;i<end;i++)
738       tf_res[i] = tf_select_table[LM][4*isTransient+2*tf_select+tf_res[i]];
739    /*for(i=0;i<end;i++)printf("%d ", isTransient ? tf_res[i] : LM+tf_res[i]);printf("\n");*/
740 }
741
742
743 static int alloc_trim_analysis(const CELTMode *m, const celt_norm *X,
744       const opus_val16 *bandLogE, int end, int LM, int C, int N0,
745       AnalysisInfo *analysis, opus_val16 *stereo_saving, opus_val16 tf_estimate,
746       int intensity)
747 {
748    int i;
749    opus_val32 diff=0;
750    int c;
751    int trim_index = 5;
752    opus_val16 trim = QCONST16(5.f, 8);
753    opus_val16 logXC, logXC2;
754    if (C==2)
755    {
756       opus_val16 sum = 0; /* Q10 */
757       opus_val16 minXC; /* Q10 */
758       /* Compute inter-channel correlation for low frequencies */
759       for (i=0;i<8;i++)
760       {
761          int j;
762          opus_val32 partial = 0;
763          for (j=m->eBands[i]<<LM;j<m->eBands[i+1]<<LM;j++)
764             partial = MAC16_16(partial, X[j], X[N0+j]);
765          sum = ADD16(sum, EXTRACT16(SHR32(partial, 18)));
766       }
767       sum = MULT16_16_Q15(QCONST16(1.f/8, 15), sum);
768       sum = MIN16(QCONST16(1.f, 10), ABS16(sum));
769       minXC = sum;
770       for (i=8;i<intensity;i++)
771       {
772          int j;
773          opus_val32 partial = 0;
774          for (j=m->eBands[i]<<LM;j<m->eBands[i+1]<<LM;j++)
775             partial = MAC16_16(partial, X[j], X[N0+j]);
776          minXC = MIN16(minXC, ABS16(EXTRACT16(SHR32(partial, 18))));
777       }
778       minXC = MIN16(QCONST16(1.f, 10), ABS16(minXC));
779       /*printf ("%f\n", sum);*/
780       if (sum > QCONST16(.995f,10))
781          trim_index-=4;
782       else if (sum > QCONST16(.92f,10))
783          trim_index-=3;
784       else if (sum > QCONST16(.85f,10))
785          trim_index-=2;
786       else if (sum > QCONST16(.8f,10))
787          trim_index-=1;
788       /* mid-side savings estimations based on the LF average*/
789       logXC = celt_log2(QCONST32(1.001f, 20)-MULT16_16(sum, sum));
790       /* mid-side savings estimations based on min correlation */
791       logXC2 = MAX16(HALF16(logXC), celt_log2(QCONST32(1.001f, 20)-MULT16_16(minXC, minXC)));
792 #ifdef FIXED_POINT
793       /* Compensate for Q20 vs Q14 input and convert output to Q8 */
794       logXC = PSHR32(logXC-QCONST16(6.f, DB_SHIFT),DB_SHIFT-8);
795       logXC2 = PSHR32(logXC2-QCONST16(6.f, DB_SHIFT),DB_SHIFT-8);
796 #endif
797
798       trim += MAX16(-QCONST16(4.f, 8), MULT16_16_Q15(QCONST16(.75f,15),logXC));
799       *stereo_saving = MIN16(*stereo_saving + QCONST16(0.25f, 8), -HALF16(logXC2));
800    }
801
802    /* Estimate spectral tilt */
803    c=0; do {
804       for (i=0;i<end-1;i++)
805       {
806          diff += bandLogE[i+c*m->nbEBands]*(opus_int32)(2+2*i-end);
807       }
808    } while (++c<C);
809    diff /= C*(end-1);
810    /*printf("%f\n", diff);*/
811    if (diff > QCONST16(2.f, DB_SHIFT))
812       trim_index--;
813    if (diff > QCONST16(8.f, DB_SHIFT))
814       trim_index--;
815    if (diff < -QCONST16(4.f, DB_SHIFT))
816       trim_index++;
817    if (diff < -QCONST16(10.f, DB_SHIFT))
818       trim_index++;
819    trim -= MAX16(-QCONST16(2.f, 8), MIN16(QCONST16(2.f, 8), SHR16(diff+QCONST16(1.f, DB_SHIFT),DB_SHIFT-8)/6 ));
820    trim -= 2*SHR16(tf_estimate, 14-8);
821 #ifndef FIXED_POINT
822    if (analysis->valid)
823    {
824       trim -= MAX16(-QCONST16(2.f, 8), MIN16(QCONST16(2.f, 8), 2*(analysis->tonality_slope+.05f)));
825    }
826 #endif
827
828 #ifdef FIXED_POINT
829    trim_index = PSHR32(trim, 8);
830 #else
831    trim_index = (int)floor(.5f+trim);
832 #endif
833    if (trim_index<0)
834       trim_index = 0;
835    if (trim_index>10)
836       trim_index = 10;
837    /*printf("%d\n", trim_index);*/
838 #ifdef FUZZING
839    trim_index = rand()%11;
840 #endif
841    return trim_index;
842 }
843
844 static int stereo_analysis(const CELTMode *m, const celt_norm *X,
845       int LM, int N0)
846 {
847    int i;
848    int thetas;
849    opus_val32 sumLR = EPSILON, sumMS = EPSILON;
850
851    /* Use the L1 norm to model the entropy of the L/R signal vs the M/S signal */
852    for (i=0;i<13;i++)
853    {
854       int j;
855       for (j=m->eBands[i]<<LM;j<m->eBands[i+1]<<LM;j++)
856       {
857          opus_val32 L, R, M, S;
858          /* We cast to 32-bit first because of the -32768 case */
859          L = EXTEND32(X[j]);
860          R = EXTEND32(X[N0+j]);
861          M = ADD32(L, R);
862          S = SUB32(L, R);
863          sumLR = ADD32(sumLR, ADD32(ABS32(L), ABS32(R)));
864          sumMS = ADD32(sumMS, ADD32(ABS32(M), ABS32(S)));
865       }
866    }
867    sumMS = MULT16_32_Q15(QCONST16(0.707107f, 15), sumMS);
868    thetas = 13;
869    /* We don't need thetas for lower bands with LM<=1 */
870    if (LM<=1)
871       thetas -= 8;
872    return MULT16_32_Q15((m->eBands[13]<<(LM+1))+thetas, sumMS)
873          > MULT16_32_Q15(m->eBands[13]<<(LM+1), sumLR);
874 }
875
876 static opus_val16 dynalloc_analysis(const opus_val16 *bandLogE, const opus_val16 *bandLogE2,
877       int nbEBands, int start, int end, int C, int *offsets, int lsb_depth, const opus_int16 *logN,
878       int isTransient, int vbr, int constrained_vbr, const opus_int16 *eBands, int LM,
879       int effectiveBytes, opus_int32 *tot_boost_, int lfe)
880 {
881    int i, c;
882    opus_int32 tot_boost=0;
883    opus_val16 maxDepth;
884    VARDECL(opus_val16, follower);
885    VARDECL(opus_val16, noise_floor);
886    SAVE_STACK;
887    ALLOC(follower, C*nbEBands, opus_val16);
888    ALLOC(noise_floor, C*nbEBands, opus_val16);
889    for (i=0;i<nbEBands;i++)
890       offsets[i] = 0;
891    /* Dynamic allocation code */
892    maxDepth=-QCONST16(32.f, DB_SHIFT);
893    for (i=0;i<end;i++)
894    {
895       /* Noise floor must take into account eMeans, the depth, the width of the bands
896          and the preemphasis filter (approx. square of bark band ID) */
897       noise_floor[i] = MULT16_16(QCONST16(0.0625f, DB_SHIFT),logN[i])
898             +QCONST16(.5f,DB_SHIFT)+SHL16(9-lsb_depth,DB_SHIFT)-SHL16(eMeans[i],6)
899             +MULT16_16(QCONST16(.0062,DB_SHIFT),(i+5)*(i+5));
900    }
901    c=0;do
902    {
903       for (i=0;i<end;i++)
904          maxDepth = MAX16(maxDepth, bandLogE[c*nbEBands+i]-noise_floor[i]);
905    } while (++c<C);
906    /* Make sure that dynamic allocation can't make us bust the budget */
907    if (effectiveBytes > 50 && LM>=1 && !lfe)
908    {
909       int last=0;
910       c=0;do
911       {
912          follower[c*nbEBands] = bandLogE2[c*nbEBands];
913          for (i=1;i<end;i++)
914          {
915             /* The last band to be at least 3 dB higher than the previous one
916                is the last we'll consider. Otherwise, we run into problems on
917                bandlimited signals. */
918             if (bandLogE2[c*nbEBands+i] > bandLogE2[c*nbEBands+i-1]+QCONST16(.5f,DB_SHIFT))
919                last=i;
920             follower[c*nbEBands+i] = MIN16(follower[c*nbEBands+i-1]+QCONST16(1.5f,DB_SHIFT), bandLogE2[c*nbEBands+i]);
921          }
922          for (i=last-1;i>=0;i--)
923             follower[c*nbEBands+i] = MIN16(follower[c*nbEBands+i], MIN16(follower[c*nbEBands+i+1]+QCONST16(2.f,DB_SHIFT), bandLogE2[c*nbEBands+i]));
924          for (i=0;i<end;i++)
925             follower[c*nbEBands+i] = MAX16(follower[c*nbEBands+i], noise_floor[i]);
926       } while (++c<C);
927       if (C==2)
928       {
929          for (i=start;i<end;i++)
930          {
931             /* Consider 24 dB "cross-talk" */
932             follower[nbEBands+i] = MAX16(follower[nbEBands+i], follower[         i]-QCONST16(4.f,DB_SHIFT));
933             follower[         i] = MAX16(follower[         i], follower[nbEBands+i]-QCONST16(4.f,DB_SHIFT));
934             follower[i] = HALF16(MAX16(0, bandLogE[i]-follower[i]) + MAX16(0, bandLogE[nbEBands+i]-follower[nbEBands+i]));
935          }
936       } else {
937          for (i=start;i<end;i++)
938          {
939             follower[i] = MAX16(0, bandLogE[i]-follower[i]);
940          }
941       }
942       /* For non-transient CBR/CVBR frames, halve the dynalloc contribution */
943       if ((!vbr || constrained_vbr)&&!isTransient)
944       {
945          for (i=start;i<end;i++)
946             follower[i] = HALF16(follower[i]);
947       }
948       for (i=start;i<end;i++)
949       {
950          int width;
951          int boost;
952          int boost_bits;
953
954          if (i<8)
955             follower[i] *= 2;
956          if (i>=12)
957             follower[i] = HALF16(follower[i]);
958          follower[i] = MIN16(follower[i], QCONST16(4, DB_SHIFT));
959
960          width = C*(eBands[i+1]-eBands[i])<<LM;
961          if (width<6)
962          {
963             boost = (int)SHR32(EXTEND32(follower[i]),DB_SHIFT);
964             boost_bits = boost*width<<BITRES;
965          } else if (width > 48) {
966             boost = (int)SHR32(EXTEND32(follower[i])*8,DB_SHIFT);
967             boost_bits = (boost*width<<BITRES)/8;
968          } else {
969             boost = (int)SHR32(EXTEND32(follower[i])*width/6,DB_SHIFT);
970             boost_bits = boost*6<<BITRES;
971          }
972          /* For CBR and non-transient CVBR frames, limit dynalloc to 1/4 of the bits */
973          if ((!vbr || (constrained_vbr&&!isTransient))
974                && (tot_boost+boost_bits)>>BITRES>>3 > effectiveBytes/4)
975          {
976             offsets[i] = 0;
977             break;
978          } else {
979             offsets[i] = boost;
980             tot_boost += boost_bits;
981          }
982       }
983    }
984    *tot_boost_ = tot_boost;
985    RESTORE_STACK;
986    return maxDepth;
987 }
988
989
990 static int run_prefilter(CELTEncoder *st, celt_sig *in, celt_sig *prefilter_mem, int CC, int N,
991       int prefilter_tapset, int *pitch, opus_val16 *gain, int *qgain, int enabled, int nbAvailableBytes)
992 {
993    int c;
994    VARDECL(celt_sig, _pre);
995    celt_sig *pre[2];
996    const CELTMode *mode;
997    int pitch_index;
998    opus_val16 gain1;
999    opus_val16 pf_threshold;
1000    int pf_on;
1001    int qg;
1002    SAVE_STACK;
1003
1004    mode = st->mode;
1005    ALLOC(_pre, CC*(N+COMBFILTER_MAXPERIOD), celt_sig);
1006
1007    pre[0] = _pre;
1008    pre[1] = _pre + (N+COMBFILTER_MAXPERIOD);
1009
1010
1011    c=0; do {
1012       OPUS_COPY(pre[c], prefilter_mem+c*COMBFILTER_MAXPERIOD, COMBFILTER_MAXPERIOD);
1013       OPUS_COPY(pre[c]+COMBFILTER_MAXPERIOD, in+c*(N+st->overlap)+st->overlap, N);
1014    } while (++c<CC);
1015
1016    if (enabled)
1017    {
1018       VARDECL(opus_val16, pitch_buf);
1019       ALLOC(pitch_buf, (COMBFILTER_MAXPERIOD+N)>>1, opus_val16);
1020
1021       pitch_downsample(pre, pitch_buf, COMBFILTER_MAXPERIOD+N, CC);
1022       /* Don't search for the fir last 1.5 octave of the range because
1023          there's too many false-positives due to short-term correlation */
1024       pitch_search(pitch_buf+(COMBFILTER_MAXPERIOD>>1), pitch_buf, N,
1025             COMBFILTER_MAXPERIOD-3*COMBFILTER_MINPERIOD, &pitch_index);
1026       pitch_index = COMBFILTER_MAXPERIOD-pitch_index;
1027
1028       gain1 = remove_doubling(pitch_buf, COMBFILTER_MAXPERIOD, COMBFILTER_MINPERIOD,
1029             N, &pitch_index, st->prefilter_period, st->prefilter_gain);
1030       if (pitch_index > COMBFILTER_MAXPERIOD-2)
1031          pitch_index = COMBFILTER_MAXPERIOD-2;
1032       gain1 = MULT16_16_Q15(QCONST16(.7f,15),gain1);
1033       /*printf("%d %d %f %f\n", pitch_change, pitch_index, gain1, st->analysis.tonality);*/
1034       if (st->loss_rate>2)
1035          gain1 = HALF32(gain1);
1036       if (st->loss_rate>4)
1037          gain1 = HALF32(gain1);
1038       if (st->loss_rate>8)
1039          gain1 = 0;
1040    } else {
1041       gain1 = 0;
1042       pitch_index = COMBFILTER_MINPERIOD;
1043    }
1044
1045    /* Gain threshold for enabling the prefilter/postfilter */
1046    pf_threshold = QCONST16(.2f,15);
1047
1048    /* Adjusting the threshold based on rate and continuity */
1049    if (abs(pitch_index-st->prefilter_period)*10>pitch_index)
1050       pf_threshold += QCONST16(.2f,15);
1051    if (nbAvailableBytes<25)
1052       pf_threshold += QCONST16(.1f,15);
1053    if (nbAvailableBytes<35)
1054       pf_threshold += QCONST16(.1f,15);
1055    if (st->prefilter_gain > QCONST16(.4f,15))
1056       pf_threshold -= QCONST16(.1f,15);
1057    if (st->prefilter_gain > QCONST16(.55f,15))
1058       pf_threshold -= QCONST16(.1f,15);
1059
1060    /* Hard threshold at 0.2 */
1061    pf_threshold = MAX16(pf_threshold, QCONST16(.2f,15));
1062    if (gain1<pf_threshold)
1063    {
1064       gain1 = 0;
1065       pf_on = 0;
1066       qg = 0;
1067    } else {
1068       /*This block is not gated by a total bits check only because
1069         of the nbAvailableBytes check above.*/
1070       if (ABS16(gain1-st->prefilter_gain)<QCONST16(.1f,15))
1071          gain1=st->prefilter_gain;
1072
1073 #ifdef FIXED_POINT
1074       qg = ((gain1+1536)>>10)/3-1;
1075 #else
1076       qg = (int)floor(.5f+gain1*32/3)-1;
1077 #endif
1078       qg = IMAX(0, IMIN(7, qg));
1079       gain1 = QCONST16(0.09375f,15)*(qg+1);
1080       pf_on = 1;
1081    }
1082    /*printf("%d %f\n", pitch_index, gain1);*/
1083
1084    c=0; do {
1085       int offset = mode->shortMdctSize-st->overlap;
1086       st->prefilter_period=IMAX(st->prefilter_period, COMBFILTER_MINPERIOD);
1087       OPUS_COPY(in+c*(N+st->overlap), st->in_mem+c*(st->overlap), st->overlap);
1088       if (offset)
1089          comb_filter(in+c*(N+st->overlap)+st->overlap, pre[c]+COMBFILTER_MAXPERIOD,
1090                st->prefilter_period, st->prefilter_period, offset, -st->prefilter_gain, -st->prefilter_gain,
1091                st->prefilter_tapset, st->prefilter_tapset, NULL, 0);
1092
1093       comb_filter(in+c*(N+st->overlap)+st->overlap+offset, pre[c]+COMBFILTER_MAXPERIOD+offset,
1094             st->prefilter_period, pitch_index, N-offset, -st->prefilter_gain, -gain1,
1095             st->prefilter_tapset, prefilter_tapset, mode->window, st->overlap);
1096       OPUS_COPY(st->in_mem+c*(st->overlap), in+c*(N+st->overlap)+N, st->overlap);
1097
1098       if (N>COMBFILTER_MAXPERIOD)
1099       {
1100          OPUS_MOVE(prefilter_mem+c*COMBFILTER_MAXPERIOD, pre[c]+N, COMBFILTER_MAXPERIOD);
1101       } else {
1102          OPUS_MOVE(prefilter_mem+c*COMBFILTER_MAXPERIOD, prefilter_mem+c*COMBFILTER_MAXPERIOD+N, COMBFILTER_MAXPERIOD-N);
1103          OPUS_MOVE(prefilter_mem+c*COMBFILTER_MAXPERIOD+COMBFILTER_MAXPERIOD-N, pre[c]+COMBFILTER_MAXPERIOD, N);
1104       }
1105    } while (++c<CC);
1106
1107    RESTORE_STACK;
1108    *gain = gain1;
1109    *pitch = pitch_index;
1110    *qgain = qg;
1111    return pf_on;
1112 }
1113
1114 static int compute_vbr(const CELTMode *mode, AnalysisInfo *analysis, opus_int32 base_target,
1115       int LM, opus_int32 bitrate, int lastCodedBands, int C, int intensity,
1116       int constrained_vbr, opus_val16 stereo_saving, int tot_boost,
1117       opus_val16 tf_estimate, int pitch_change, opus_val16 maxDepth,
1118       int variable_duration, int lfe, int has_surround_mask, opus_val16 surround_masking)
1119 {
1120    /* The target rate in 8th bits per frame */
1121    opus_int32 target;
1122    int coded_bins;
1123    int coded_bands;
1124    int tf_calibration;
1125    int nbEBands;
1126    const opus_int16 *eBands;
1127
1128    nbEBands = mode->nbEBands;
1129    eBands = mode->eBands;
1130
1131    coded_bands = lastCodedBands ? lastCodedBands : nbEBands;
1132    coded_bins = eBands[coded_bands]<<LM;
1133    if (C==2)
1134       coded_bins += eBands[IMIN(intensity, coded_bands)]<<LM;
1135
1136    target = base_target;
1137
1138    /*printf("%f %f %f %f %d %d ", st->analysis.activity, st->analysis.tonality, tf_estimate, st->stereo_saving, tot_boost, coded_bands);*/
1139 #ifndef FIXED_POINT
1140    if (analysis->valid && analysis->activity<.4)
1141       target -= (opus_int32)((coded_bins<<BITRES)*(.4f-analysis->activity));
1142 #endif
1143    /* Stereo savings */
1144    if (C==2)
1145    {
1146       int coded_stereo_bands;
1147       int coded_stereo_dof;
1148       opus_val16 max_frac;
1149       coded_stereo_bands = IMIN(intensity, coded_bands);
1150       coded_stereo_dof = (eBands[coded_stereo_bands]<<LM)-coded_stereo_bands;
1151       /* Maximum fraction of the bits we can save if the signal is mono. */
1152       max_frac = DIV32_16(MULT16_16(QCONST16(0.8f, 15), coded_stereo_dof), coded_bins);
1153       /*printf("%d %d %d ", coded_stereo_dof, coded_bins, tot_boost);*/
1154       target -= (opus_int32)MIN32(MULT16_32_Q15(max_frac,target),
1155                       SHR32(MULT16_16(stereo_saving-QCONST16(0.1f,8),(coded_stereo_dof<<BITRES)),8));
1156    }
1157    /* Boost the rate according to dynalloc (minus the dynalloc average for calibration). */
1158    target += tot_boost-(16<<LM);
1159    /* Apply transient boost, compensating for average boost. */
1160    tf_calibration = variable_duration ? QCONST16(0.02f,14) : QCONST16(0.04f,14);
1161    target += (opus_int32)SHL32(MULT16_32_Q15(tf_estimate-tf_calibration, target),1);
1162
1163 #ifndef FIXED_POINT
1164    /* Apply tonality boost */
1165    if (analysis->valid) {
1166       opus_int32 tonal_target;
1167       float tonal;
1168
1169       /* Tonality boost (compensating for the average). */
1170       tonal = MAX16(0.f,analysis->tonality-.15f)-0.09f;
1171       tonal_target = target + (opus_int32)((coded_bins<<BITRES)*1.2f*tonal);
1172       if (pitch_change)
1173          tonal_target +=  (opus_int32)((coded_bins<<BITRES)*.8f);
1174       /*printf("%f %f ", st->analysis.tonality, tonal);*/
1175       target = tonal_target;
1176    }
1177 #endif
1178
1179    if (has_surround_mask&&!lfe)
1180    {
1181       opus_int32 surround_target = target + SHR32(MULT16_16(surround_masking,coded_bins<<BITRES), DB_SHIFT);
1182       /*printf("%f %d %d %d %d %d %d ", surround_masking, coded_bins, st->end, st->intensity, surround_target, target, st->bitrate);*/
1183       target = IMAX(target/4, surround_target);
1184    }
1185
1186    {
1187       opus_int32 floor_depth;
1188       int bins;
1189       bins = eBands[nbEBands-2]<<LM;
1190       /*floor_depth = SHR32(MULT16_16((C*bins<<BITRES),celt_log2(SHL32(MAX16(1,sample_max),13))), DB_SHIFT);*/
1191       floor_depth = (opus_int32)SHR32(MULT16_16((C*bins<<BITRES),maxDepth), DB_SHIFT);
1192       floor_depth = IMAX(floor_depth, target>>2);
1193       target = IMIN(target, floor_depth);
1194       /*printf("%f %d\n", maxDepth, floor_depth);*/
1195    }
1196
1197    if ((!has_surround_mask||lfe) && (constrained_vbr || bitrate<64000))
1198    {
1199       opus_val16 rate_factor;
1200 #ifdef FIXED_POINT
1201       rate_factor = MAX16(0,(bitrate-32000));
1202 #else
1203       rate_factor = MAX16(0,(1.f/32768)*(bitrate-32000));
1204 #endif
1205       if (constrained_vbr)
1206          rate_factor = MIN16(rate_factor, QCONST16(0.67f, 15));
1207       target = base_target + (opus_int32)MULT16_32_Q15(rate_factor, target-base_target);
1208
1209    }
1210    /* Don't allow more than doubling the rate */
1211    target = IMIN(2*base_target, target);
1212
1213    return target;
1214 }
1215
1216 int celt_encode_with_ec(CELTEncoder * OPUS_RESTRICT st, const opus_val16 * pcm, int frame_size, unsigned char *compressed, int nbCompressedBytes, ec_enc *enc)
1217 {
1218    int i, c, N;
1219    opus_int32 bits;
1220    ec_enc _enc;
1221    VARDECL(celt_sig, in);
1222    VARDECL(celt_sig, freq);
1223    VARDECL(celt_norm, X);
1224    VARDECL(celt_ener, bandE);
1225    VARDECL(opus_val16, bandLogE);
1226    VARDECL(opus_val16, bandLogE2);
1227    VARDECL(int, fine_quant);
1228    VARDECL(opus_val16, error);
1229    VARDECL(int, pulses);
1230    VARDECL(int, cap);
1231    VARDECL(int, offsets);
1232    VARDECL(int, fine_priority);
1233    VARDECL(int, tf_res);
1234    VARDECL(unsigned char, collapse_masks);
1235    celt_sig *prefilter_mem;
1236    opus_val16 *oldBandE, *oldLogE, *oldLogE2;
1237    int shortBlocks=0;
1238    int isTransient=0;
1239    const int CC = st->channels;
1240    const int C = st->stream_channels;
1241    int LM, M;
1242    int tf_select;
1243    int nbFilledBytes, nbAvailableBytes;
1244    int effEnd;
1245    int codedBands;
1246    int tf_sum;
1247    int alloc_trim;
1248    int pitch_index=COMBFILTER_MINPERIOD;
1249    opus_val16 gain1 = 0;
1250    int dual_stereo=0;
1251    int effectiveBytes;
1252    int dynalloc_logp;
1253    opus_int32 vbr_rate;
1254    opus_int32 total_bits;
1255    opus_int32 total_boost;
1256    opus_int32 balance;
1257    opus_int32 tell;
1258    int prefilter_tapset=0;
1259    int pf_on;
1260    int anti_collapse_rsv;
1261    int anti_collapse_on=0;
1262    int silence=0;
1263    int tf_chan = 0;
1264    opus_val16 tf_estimate;
1265    int pitch_change=0;
1266    opus_int32 tot_boost;
1267    opus_val32 sample_max;
1268    opus_val16 maxDepth;
1269    const OpusCustomMode *mode;
1270    int nbEBands;
1271    int overlap;
1272    const opus_int16 *eBands;
1273    int secondMdct;
1274    int signalBandwidth;
1275    int transient_got_disabled=0;
1276    opus_val16 surround_masking=0;
1277    ALLOC_STACK;
1278
1279    mode = st->mode;
1280    nbEBands = mode->nbEBands;
1281    overlap = mode->overlap;
1282    eBands = mode->eBands;
1283    tf_estimate = 0;
1284    if (nbCompressedBytes<2 || pcm==NULL)
1285      return OPUS_BAD_ARG;
1286
1287    frame_size *= st->upsample;
1288    for (LM=0;LM<=mode->maxLM;LM++)
1289       if (mode->shortMdctSize<<LM==frame_size)
1290          break;
1291    if (LM>mode->maxLM)
1292       return OPUS_BAD_ARG;
1293    M=1<<LM;
1294    N = M*mode->shortMdctSize;
1295
1296    prefilter_mem = st->in_mem+CC*(st->overlap);
1297    oldBandE = (opus_val16*)(st->in_mem+CC*(st->overlap+COMBFILTER_MAXPERIOD));
1298    oldLogE = oldBandE + CC*nbEBands;
1299    oldLogE2 = oldLogE + CC*nbEBands;
1300
1301    if (enc==NULL)
1302    {
1303       tell=1;
1304       nbFilledBytes=0;
1305    } else {
1306       tell=ec_tell(enc);
1307       nbFilledBytes=(tell+4)>>3;
1308    }
1309
1310 #ifdef CUSTOM_MODES
1311    if (st->signalling && enc==NULL)
1312    {
1313       int tmp = (mode->effEBands-st->end)>>1;
1314       st->end = IMAX(1, mode->effEBands-tmp);
1315       compressed[0] = tmp<<5;
1316       compressed[0] |= LM<<3;
1317       compressed[0] |= (C==2)<<2;
1318       /* Convert "standard mode" to Opus header */
1319       if (mode->Fs==48000 && mode->shortMdctSize==120)
1320       {
1321          int c0 = toOpus(compressed[0]);
1322          if (c0<0)
1323             return OPUS_BAD_ARG;
1324          compressed[0] = c0;
1325       }
1326       compressed++;
1327       nbCompressedBytes--;
1328    }
1329 #else
1330    celt_assert(st->signalling==0);
1331 #endif
1332
1333    /* Can't produce more than 1275 output bytes */
1334    nbCompressedBytes = IMIN(nbCompressedBytes,1275);
1335    nbAvailableBytes = nbCompressedBytes - nbFilledBytes;
1336
1337    if (st->vbr && st->bitrate!=OPUS_BITRATE_MAX)
1338    {
1339       opus_int32 den=mode->Fs>>BITRES;
1340       vbr_rate=(st->bitrate*frame_size+(den>>1))/den;
1341 #ifdef CUSTOM_MODES
1342       if (st->signalling)
1343          vbr_rate -= 8<<BITRES;
1344 #endif
1345       effectiveBytes = vbr_rate>>(3+BITRES);
1346    } else {
1347       opus_int32 tmp;
1348       vbr_rate = 0;
1349       tmp = st->bitrate*frame_size;
1350       if (tell>1)
1351          tmp += tell;
1352       if (st->bitrate!=OPUS_BITRATE_MAX)
1353          nbCompressedBytes = IMAX(2, IMIN(nbCompressedBytes,
1354                (tmp+4*mode->Fs)/(8*mode->Fs)-!!st->signalling));
1355       effectiveBytes = nbCompressedBytes;
1356    }
1357
1358    if (enc==NULL)
1359    {
1360       ec_enc_init(&_enc, compressed, nbCompressedBytes);
1361       enc = &_enc;
1362    }
1363
1364    if (vbr_rate>0)
1365    {
1366       /* Computes the max bit-rate allowed in VBR mode to avoid violating the
1367           target rate and buffering.
1368          We must do this up front so that bust-prevention logic triggers
1369           correctly if we don't have enough bits. */
1370       if (st->constrained_vbr)
1371       {
1372          opus_int32 vbr_bound;
1373          opus_int32 max_allowed;
1374          /* We could use any multiple of vbr_rate as bound (depending on the
1375              delay).
1376             This is clamped to ensure we use at least two bytes if the encoder
1377              was entirely empty, but to allow 0 in hybrid mode. */
1378          vbr_bound = vbr_rate;
1379          max_allowed = IMIN(IMAX(tell==1?2:0,
1380                (vbr_rate+vbr_bound-st->vbr_reservoir)>>(BITRES+3)),
1381                nbAvailableBytes);
1382          if(max_allowed < nbAvailableBytes)
1383          {
1384             nbCompressedBytes = nbFilledBytes+max_allowed;
1385             nbAvailableBytes = max_allowed;
1386             ec_enc_shrink(enc, nbCompressedBytes);
1387          }
1388       }
1389    }
1390    total_bits = nbCompressedBytes*8;
1391
1392    effEnd = st->end;
1393    if (effEnd > mode->effEBands)
1394       effEnd = mode->effEBands;
1395
1396    ALLOC(in, CC*(N+st->overlap), celt_sig);
1397
1398    sample_max=MAX32(st->overlap_max, celt_maxabs16(pcm, C*(N-overlap)/st->upsample));
1399    st->overlap_max=celt_maxabs16(pcm+C*(N-overlap)/st->upsample, C*overlap/st->upsample);
1400    sample_max=MAX32(sample_max, st->overlap_max);
1401 #ifdef FIXED_POINT
1402    silence = (sample_max==0);
1403 #else
1404    silence = (sample_max <= (opus_val16)1/(1<<st->lsb_depth));
1405 #endif
1406 #ifdef FUZZING
1407    if ((rand()&0x3F)==0)
1408       silence = 1;
1409 #endif
1410    if (tell==1)
1411       ec_enc_bit_logp(enc, silence, 15);
1412    else
1413       silence=0;
1414    if (silence)
1415    {
1416       /*In VBR mode there is no need to send more than the minimum. */
1417       if (vbr_rate>0)
1418       {
1419          effectiveBytes=nbCompressedBytes=IMIN(nbCompressedBytes, nbFilledBytes+2);
1420          total_bits=nbCompressedBytes*8;
1421          nbAvailableBytes=2;
1422          ec_enc_shrink(enc, nbCompressedBytes);
1423       }
1424       /* Pretend we've filled all the remaining bits with zeros
1425             (that's what the initialiser did anyway) */
1426       tell = nbCompressedBytes*8;
1427       enc->nbits_total+=tell-ec_tell(enc);
1428    }
1429    c=0; do {
1430       preemphasis(pcm+c, in+c*(N+st->overlap)+st->overlap, N, CC, st->upsample,
1431                   mode->preemph, st->preemph_memE+c, st->clip);
1432    } while (++c<CC);
1433
1434
1435
1436    /* Find pitch period and gain */
1437    {
1438       int enabled;
1439       int qg;
1440       enabled = (st->lfe || nbAvailableBytes>12*C) && st->start==0 && !silence && !st->disable_pf
1441             && st->complexity >= 5 && !(st->consec_transient && LM!=3 && st->variable_duration);
1442
1443       prefilter_tapset = st->tapset_decision;
1444       pf_on = run_prefilter(st, in, prefilter_mem, CC, N, prefilter_tapset, &pitch_index, &gain1, &qg, enabled, nbAvailableBytes);
1445       if ((gain1 > QCONST16(.4f,15) || st->prefilter_gain > QCONST16(.4f,15)) && (!st->analysis.valid || st->analysis.tonality > .3)
1446             && (pitch_index > 1.26*st->prefilter_period || pitch_index < .79*st->prefilter_period))
1447          pitch_change = 1;
1448       if (pf_on==0)
1449       {
1450          if(st->start==0 && tell+16<=total_bits)
1451             ec_enc_bit_logp(enc, 0, 1);
1452       } else {
1453          /*This block is not gated by a total bits check only because
1454            of the nbAvailableBytes check above.*/
1455          int octave;
1456          ec_enc_bit_logp(enc, 1, 1);
1457          pitch_index += 1;
1458          octave = EC_ILOG(pitch_index)-5;
1459          ec_enc_uint(enc, octave, 6);
1460          ec_enc_bits(enc, pitch_index-(16<<octave), 4+octave);
1461          pitch_index -= 1;
1462          ec_enc_bits(enc, qg, 3);
1463          ec_enc_icdf(enc, prefilter_tapset, tapset_icdf, 2);
1464       }
1465    }
1466
1467    isTransient = 0;
1468    shortBlocks = 0;
1469    if (st->complexity >= 1 && !st->lfe)
1470    {
1471       isTransient = transient_analysis(in, N+st->overlap, CC,
1472             &tf_estimate, &tf_chan);
1473    }
1474    if (LM>0 && ec_tell(enc)+3<=total_bits)
1475    {
1476       if (isTransient)
1477          shortBlocks = M;
1478    } else {
1479       isTransient = 0;
1480       transient_got_disabled=1;
1481    }
1482
1483    ALLOC(freq, CC*N, celt_sig); /**< Interleaved signal MDCTs */
1484    ALLOC(bandE,nbEBands*CC, celt_ener);
1485    ALLOC(bandLogE,nbEBands*CC, opus_val16);
1486
1487    secondMdct = shortBlocks && st->complexity>=8;
1488    ALLOC(bandLogE2, C*nbEBands, opus_val16);
1489    if (secondMdct)
1490    {
1491       compute_mdcts(mode, 0, in, freq, C, CC, LM, st->upsample);
1492       compute_band_energies(mode, freq, bandE, effEnd, C, M);
1493       amp2Log2(mode, effEnd, st->end, bandE, bandLogE2, C);
1494       for (i=0;i<C*nbEBands;i++)
1495          bandLogE2[i] += HALF16(SHL16(LM, DB_SHIFT));
1496    }
1497
1498    compute_mdcts(mode, shortBlocks, in, freq, C, CC, LM, st->upsample);
1499    if (CC==2&&C==1)
1500       tf_chan = 0;
1501    compute_band_energies(mode, freq, bandE, effEnd, C, M);
1502
1503    if (st->lfe)
1504    {
1505       for (i=2;i<st->end;i++)
1506          bandE[i] = IMIN(bandE[i], MULT16_32_Q15(QCONST16(1e-4f,15),bandE[0]));
1507    }
1508    amp2Log2(mode, effEnd, st->end, bandE, bandLogE, C);
1509    if (st->energy_save)
1510    {
1511       opus_val16 offset = shortBlocks?HALF16(SHL16(LM, DB_SHIFT)):0;
1512 #ifdef FIXED_POINT
1513       /* Compensate for the 1/8 gain we apply in the fixed-point downshift to avoid overflows. */
1514       offset -= QCONST16(3.0f, DB_SHIFT);
1515 #endif
1516       for(i=0;i<C*nbEBands;i++)
1517          st->energy_save[i]=bandLogE[i]-offset;
1518       st->energy_save=NULL;
1519    }
1520    /* This computes how much masking takes place between surround channels */
1521    if (st->energy_mask&&!st->lfe)
1522    {
1523       opus_val32 mask_avg=0;
1524       opus_val16 offset = shortBlocks?HALF16(SHL16(LM, DB_SHIFT)):0;
1525       for (c=0;c<C;c++)
1526       {
1527          opus_val16 followE, followMask;
1528          followE = followMask = -QCONST16(14.f, DB_SHIFT);
1529          for(i=0;i<st->end;i++)
1530          {
1531             /* We use a simple follower to approximate the masking spreading function. */
1532             followE = MAX16(followE-QCONST16(1.f, DB_SHIFT), bandLogE[nbEBands*c+i]-offset);
1533             followMask = MAX16(followMask-QCONST16(1.f, DB_SHIFT), st->energy_mask[nbEBands*c+i]);
1534             mask_avg += followE-followMask;
1535          }
1536       }
1537       surround_masking = DIV32_16(mask_avg,C*st->end) + QCONST16(.0f, DB_SHIFT);
1538       surround_masking = MIN16(MAX16(surround_masking,-QCONST16(1.5f, DB_SHIFT)), 0);
1539    }
1540    /*for (i=0;i<21;i++)
1541       printf("%f ", bandLogE[i]);
1542    printf("\n");*/
1543
1544    if (!secondMdct)
1545    {
1546       for (i=0;i<C*nbEBands;i++)
1547          bandLogE2[i] = bandLogE[i];
1548    }
1549
1550    /* Last chance to catch any transient we might have missed in the
1551       time-domain analysis */
1552    if (LM>0 && ec_tell(enc)+3<=total_bits && !isTransient && st->complexity>=5 && !st->lfe)
1553    {
1554       if (patch_transient_decision(bandLogE, oldBandE, nbEBands, st->end, C))
1555       {
1556          isTransient = 1;
1557          shortBlocks = M;
1558          compute_mdcts(mode, shortBlocks, in, freq, C, CC, LM, st->upsample);
1559          compute_band_energies(mode, freq, bandE, effEnd, C, M);
1560          amp2Log2(mode, effEnd, st->end, bandE, bandLogE, C);
1561          /* Compensate for the scaling of short vs long mdcts */
1562          for (i=0;i<C*nbEBands;i++)
1563             bandLogE2[i] += HALF16(SHL16(LM, DB_SHIFT));
1564          tf_estimate = QCONST16(.2,14);
1565       }
1566    }
1567
1568    if (LM>0 && ec_tell(enc)+3<=total_bits)
1569       ec_enc_bit_logp(enc, isTransient, 3);
1570
1571    ALLOC(X, C*N, celt_norm);         /**< Interleaved normalised MDCTs */
1572
1573    /* Band normalisation */
1574    normalise_bands(mode, freq, X, bandE, effEnd, C, M);
1575
1576    ALLOC(tf_res, nbEBands, int);
1577    /* Disable variable tf resolution for hybrid and at very low bitrate */
1578    if (effectiveBytes>=15*C && st->start==0 && st->complexity>=2 && !st->lfe)
1579    {
1580       int lambda;
1581       if (effectiveBytes<40)
1582          lambda = 12;
1583       else if (effectiveBytes<60)
1584          lambda = 6;
1585       else if (effectiveBytes<100)
1586          lambda = 4;
1587       else
1588          lambda = 3;
1589       lambda*=2;
1590       tf_select = tf_analysis(mode, effEnd, isTransient, tf_res, lambda, X, N, LM, &tf_sum, tf_estimate, tf_chan);
1591       for (i=effEnd;i<st->end;i++)
1592          tf_res[i] = tf_res[effEnd-1];
1593    } else {
1594       tf_sum = 0;
1595       for (i=0;i<st->end;i++)
1596          tf_res[i] = isTransient;
1597       tf_select=0;
1598    }
1599
1600    ALLOC(error, C*nbEBands, opus_val16);
1601    quant_coarse_energy(mode, st->start, st->end, effEnd, bandLogE,
1602          oldBandE, total_bits, error, enc,
1603          C, LM, nbAvailableBytes, st->force_intra,
1604          &st->delayedIntra, st->complexity >= 4, st->loss_rate, st->lfe);
1605
1606    tf_encode(st->start, st->end, isTransient, tf_res, LM, tf_select, enc);
1607
1608    if (ec_tell(enc)+4<=total_bits)
1609    {
1610       if (st->lfe)
1611       {
1612          st->tapset_decision = 0;
1613          st->spread_decision = SPREAD_NORMAL;
1614       } else if (shortBlocks || st->complexity < 3 || nbAvailableBytes < 10*C || st->start != 0)
1615       {
1616          if (st->complexity == 0)
1617             st->spread_decision = SPREAD_NONE;
1618          else
1619             st->spread_decision = SPREAD_NORMAL;
1620       } else {
1621          /* Disable new spreading+tapset estimator until we can show it works
1622             better than the old one. So far it seems like spreading_decision()
1623             works best. */
1624          if (0&&st->analysis.valid)
1625          {
1626             static const opus_val16 spread_thresholds[3] = {-QCONST16(.6f, 15), -QCONST16(.2f, 15), -QCONST16(.07f, 15)};
1627             static const opus_val16 spread_histeresis[3] = {QCONST16(.15f, 15), QCONST16(.07f, 15), QCONST16(.02f, 15)};
1628             static const opus_val16 tapset_thresholds[2] = {QCONST16(.0f, 15), QCONST16(.15f, 15)};
1629             static const opus_val16 tapset_histeresis[2] = {QCONST16(.1f, 15), QCONST16(.05f, 15)};
1630             st->spread_decision = hysteresis_decision(-st->analysis.tonality, spread_thresholds, spread_histeresis, 3, st->spread_decision);
1631             st->tapset_decision = hysteresis_decision(st->analysis.tonality_slope, tapset_thresholds, tapset_histeresis, 2, st->tapset_decision);
1632          } else {
1633             st->spread_decision = spreading_decision(mode, X,
1634                   &st->tonal_average, st->spread_decision, &st->hf_average,
1635                   &st->tapset_decision, pf_on&&!shortBlocks, effEnd, C, M);
1636          }
1637          /*printf("%d %d\n", st->tapset_decision, st->spread_decision);*/
1638          /*printf("%f %d %f %d\n\n", st->analysis.tonality, st->spread_decision, st->analysis.tonality_slope, st->tapset_decision);*/
1639       }
1640       ec_enc_icdf(enc, st->spread_decision, spread_icdf, 5);
1641    }
1642
1643    ALLOC(offsets, nbEBands, int);
1644
1645    maxDepth = dynalloc_analysis(bandLogE, bandLogE2, nbEBands, st->start, st->end, C, offsets,
1646          st->lsb_depth, mode->logN, isTransient, st->vbr, st->constrained_vbr,
1647          eBands, LM, effectiveBytes, &tot_boost, st->lfe);
1648    /* For LFE, everything interesting is in the first band */
1649    if (st->lfe)
1650       offsets[0] = IMIN(8, effectiveBytes/3);
1651    ALLOC(cap, nbEBands, int);
1652    init_caps(mode,cap,LM,C);
1653
1654    dynalloc_logp = 6;
1655    total_bits<<=BITRES;
1656    total_boost = 0;
1657    tell = ec_tell_frac(enc);
1658    for (i=st->start;i<st->end;i++)
1659    {
1660       int width, quanta;
1661       int dynalloc_loop_logp;
1662       int boost;
1663       int j;
1664       width = C*(eBands[i+1]-eBands[i])<<LM;
1665       /* quanta is 6 bits, but no more than 1 bit/sample
1666          and no less than 1/8 bit/sample */
1667       quanta = IMIN(width<<BITRES, IMAX(6<<BITRES, width));
1668       dynalloc_loop_logp = dynalloc_logp;
1669       boost = 0;
1670       for (j = 0; tell+(dynalloc_loop_logp<<BITRES) < total_bits-total_boost
1671             && boost < cap[i]; j++)
1672       {
1673          int flag;
1674          flag = j<offsets[i];
1675          ec_enc_bit_logp(enc, flag, dynalloc_loop_logp);
1676          tell = ec_tell_frac(enc);
1677          if (!flag)
1678             break;
1679          boost += quanta;
1680          total_boost += quanta;
1681          dynalloc_loop_logp = 1;
1682       }
1683       /* Making dynalloc more likely */
1684       if (j)
1685          dynalloc_logp = IMAX(2, dynalloc_logp-1);
1686       offsets[i] = boost;
1687    }
1688
1689    if (C==2)
1690    {
1691       int effectiveRate;
1692
1693       static const opus_val16 intensity_thresholds[21]=
1694       /* 0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19  20  off*/
1695         { 16,21,23,25,27,29,31,33,35,38,42,46,50,54,58,63,68,75,84,102,130};
1696       static const opus_val16 intensity_histeresis[21]=
1697         {  2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 4, 5, 6,  8, 12};
1698
1699       /* Always use MS for 2.5 ms frames until we can do a better analysis */
1700       if (LM!=0)
1701          dual_stereo = stereo_analysis(mode, X, LM, N);
1702
1703       /* Account for coarse energy */
1704       effectiveRate = (8*effectiveBytes - 80)>>LM;
1705
1706       /* effectiveRate in kb/s */
1707       effectiveRate = 2*effectiveRate/5;
1708
1709       st->intensity = hysteresis_decision((opus_val16)effectiveRate, intensity_thresholds, intensity_histeresis, 21, st->intensity);
1710       st->intensity = IMIN(st->end,IMAX(st->start, st->intensity));
1711    }
1712
1713    alloc_trim = 5;
1714    if (tell+(6<<BITRES) <= total_bits - total_boost)
1715    {
1716       if (st->lfe)
1717          alloc_trim = 5;
1718       else
1719          alloc_trim = alloc_trim_analysis(mode, X, bandLogE,
1720             st->end, LM, C, N, &st->analysis, &st->stereo_saving, tf_estimate, st->intensity);
1721       ec_enc_icdf(enc, alloc_trim, trim_icdf, 7);
1722       tell = ec_tell_frac(enc);
1723    }
1724
1725    /* Variable bitrate */
1726    if (vbr_rate>0)
1727    {
1728      opus_val16 alpha;
1729      opus_int32 delta;
1730      /* The target rate in 8th bits per frame */
1731      opus_int32 target, base_target;
1732      opus_int32 min_allowed;
1733      int lm_diff = mode->maxLM - LM;
1734
1735      /* Don't attempt to use more than 510 kb/s, even for frames smaller than 20 ms.
1736         The CELT allocator will just not be able to use more than that anyway. */
1737      nbCompressedBytes = IMIN(nbCompressedBytes,1275>>(3-LM));
1738      base_target = vbr_rate - ((40*C+20)<<BITRES);
1739
1740      if (st->constrained_vbr)
1741         base_target += (st->vbr_offset>>lm_diff);
1742
1743      target = compute_vbr(mode, &st->analysis, base_target, LM, st->bitrate,
1744            st->lastCodedBands, C, st->intensity, st->constrained_vbr,
1745            st->stereo_saving, tot_boost, tf_estimate, pitch_change, maxDepth,
1746            st->variable_duration, st->lfe, st->energy_mask!=NULL, surround_masking);
1747
1748      /* The current offset is removed from the target and the space used
1749         so far is added*/
1750      target=target+tell;
1751      /* In VBR mode the frame size must not be reduced so much that it would
1752          result in the encoder running out of bits.
1753         The margin of 2 bytes ensures that none of the bust-prevention logic
1754          in the decoder will have triggered so far. */
1755      min_allowed = ((tell+total_boost+(1<<(BITRES+3))-1)>>(BITRES+3)) + 2 - nbFilledBytes;
1756
1757      nbAvailableBytes = (target+(1<<(BITRES+2)))>>(BITRES+3);
1758      nbAvailableBytes = IMAX(min_allowed,nbAvailableBytes);
1759      nbAvailableBytes = IMIN(nbCompressedBytes,nbAvailableBytes+nbFilledBytes) - nbFilledBytes;
1760
1761      /* By how much did we "miss" the target on that frame */
1762      delta = target - vbr_rate;
1763
1764      target=nbAvailableBytes<<(BITRES+3);
1765
1766      /*If the frame is silent we don't adjust our drift, otherwise
1767        the encoder will shoot to very high rates after hitting a
1768        span of silence, but we do allow the bitres to refill.
1769        This means that we'll undershoot our target in CVBR/VBR modes
1770        on files with lots of silence. */
1771      if(silence)
1772      {
1773        nbAvailableBytes = 2;
1774        target = 2*8<<BITRES;
1775        delta = 0;
1776      }
1777
1778      if (st->vbr_count < 970)
1779      {
1780         st->vbr_count++;
1781         alpha = celt_rcp(SHL32(EXTEND32(st->vbr_count+20),16));
1782      } else
1783         alpha = QCONST16(.001f,15);
1784      /* How many bits have we used in excess of what we're allowed */
1785      if (st->constrained_vbr)
1786         st->vbr_reservoir += target - vbr_rate;
1787      /*printf ("%d\n", st->vbr_reservoir);*/
1788
1789      /* Compute the offset we need to apply in order to reach the target */
1790      if (st->constrained_vbr)
1791      {
1792         st->vbr_drift += (opus_int32)MULT16_32_Q15(alpha,(delta*(1<<lm_diff))-st->vbr_offset-st->vbr_drift);
1793         st->vbr_offset = -st->vbr_drift;
1794      }
1795      /*printf ("%d\n", st->vbr_drift);*/
1796
1797      if (st->constrained_vbr && st->vbr_reservoir < 0)
1798      {
1799         /* We're under the min value -- increase rate */
1800         int adjust = (-st->vbr_reservoir)/(8<<BITRES);
1801         /* Unless we're just coding silence */
1802         nbAvailableBytes += silence?0:adjust;
1803         st->vbr_reservoir = 0;
1804         /*printf ("+%d\n", adjust);*/
1805      }
1806      nbCompressedBytes = IMIN(nbCompressedBytes,nbAvailableBytes+nbFilledBytes);
1807      /*printf("%d\n", nbCompressedBytes*50*8);*/
1808      /* This moves the raw bits to take into account the new compressed size */
1809      ec_enc_shrink(enc, nbCompressedBytes);
1810    }
1811
1812    /* Bit allocation */
1813    ALLOC(fine_quant, nbEBands, int);
1814    ALLOC(pulses, nbEBands, int);
1815    ALLOC(fine_priority, nbEBands, int);
1816
1817    /* bits =           packet size                    - where we are - safety*/
1818    bits = (((opus_int32)nbCompressedBytes*8)<<BITRES) - ec_tell_frac(enc) - 1;
1819    anti_collapse_rsv = isTransient&&LM>=2&&bits>=((LM+2)<<BITRES) ? (1<<BITRES) : 0;
1820    bits -= anti_collapse_rsv;
1821    signalBandwidth = st->end-1;
1822 #ifndef FIXED_POINT
1823    if (st->analysis.valid)
1824    {
1825       int min_bandwidth;
1826       if (st->bitrate < (opus_int32)32000*C)
1827          min_bandwidth = 13;
1828       else if (st->bitrate < (opus_int32)48000*C)
1829          min_bandwidth = 16;
1830       else if (st->bitrate < (opus_int32)60000*C)
1831          min_bandwidth = 18;
1832       else  if (st->bitrate < (opus_int32)80000*C)
1833          min_bandwidth = 19;
1834       else
1835          min_bandwidth = 20;
1836       signalBandwidth = IMAX(st->analysis.bandwidth, min_bandwidth);
1837    }
1838 #endif
1839    if (st->lfe)
1840       signalBandwidth = 1;
1841    codedBands = compute_allocation(mode, st->start, st->end, offsets, cap,
1842          alloc_trim, &st->intensity, &dual_stereo, bits, &balance, pulses,
1843          fine_quant, fine_priority, C, LM, enc, 1, st->lastCodedBands, signalBandwidth);
1844    if (st->lastCodedBands)
1845       st->lastCodedBands = IMIN(st->lastCodedBands+1,IMAX(st->lastCodedBands-1,codedBands));
1846    else
1847       st->lastCodedBands = codedBands;
1848
1849    quant_fine_energy(mode, st->start, st->end, oldBandE, error, fine_quant, enc, C);
1850
1851 #ifdef MEASURE_NORM_MSE
1852    float X0[3000];
1853    float bandE0[60];
1854    c=0; do
1855       for (i=0;i<N;i++)
1856          X0[i+c*N] = X[i+c*N];
1857    while (++c<C);
1858    for (i=0;i<C*nbEBands;i++)
1859       bandE0[i] = bandE[i];
1860 #endif
1861
1862    /* Residual quantisation */
1863    ALLOC(collapse_masks, C*nbEBands, unsigned char);
1864    quant_all_bands(1, mode, st->start, st->end, X, C==2 ? X+N : NULL, collapse_masks,
1865          bandE, pulses, shortBlocks, st->spread_decision, dual_stereo, st->intensity, tf_res,
1866          nbCompressedBytes*(8<<BITRES)-anti_collapse_rsv, balance, enc, LM, codedBands, &st->rng);
1867
1868    if (anti_collapse_rsv > 0)
1869    {
1870       anti_collapse_on = st->consec_transient<2;
1871 #ifdef FUZZING
1872       anti_collapse_on = rand()&0x1;
1873 #endif
1874       ec_enc_bits(enc, anti_collapse_on, 1);
1875    }
1876    quant_energy_finalise(mode, st->start, st->end, oldBandE, error, fine_quant, fine_priority, nbCompressedBytes*8-ec_tell(enc), enc, C);
1877
1878    if (silence)
1879    {
1880       for (i=0;i<C*nbEBands;i++)
1881          oldBandE[i] = -QCONST16(28.f,DB_SHIFT);
1882    }
1883
1884 #ifdef RESYNTH
1885    /* Re-synthesis of the coded audio if required */
1886    {
1887       celt_sig *out_mem[2];
1888
1889       log2Amp(mode, st->start, st->end, bandE, oldBandE, C);
1890       if (silence)
1891       {
1892          for (i=0;i<C*nbEBands;i++)
1893             bandE[i] = 0;
1894       }
1895
1896 #ifdef MEASURE_NORM_MSE
1897       measure_norm_mse(mode, X, X0, bandE, bandE0, M, N, C);
1898 #endif
1899       if (anti_collapse_on)
1900       {
1901          anti_collapse(mode, X, collapse_masks, LM, C, N,
1902                st->start, st->end, oldBandE, oldLogE, oldLogE2, pulses, st->rng);
1903       }
1904
1905       /* Synthesis */
1906       denormalise_bands(mode, X, freq, bandE, st->start, effEnd, C, M);
1907
1908       c=0; do {
1909          OPUS_MOVE(st->syn_mem[c], st->syn_mem[c]+N, 2*MAX_PERIOD-N+overlap/2);
1910       } while (++c<CC);
1911
1912       if (CC==2&&C==1)
1913       {
1914          for (i=0;i<N;i++)
1915             freq[N+i] = freq[i];
1916       }
1917
1918       c=0; do {
1919          out_mem[c] = st->syn_mem[c]+2*MAX_PERIOD-N;
1920       } while (++c<CC);
1921
1922       compute_inv_mdcts(mode, shortBlocks, freq, out_mem, CC, LM);
1923
1924       c=0; do {
1925          st->prefilter_period=IMAX(st->prefilter_period, COMBFILTER_MINPERIOD);
1926          st->prefilter_period_old=IMAX(st->prefilter_period_old, COMBFILTER_MINPERIOD);
1927          comb_filter(out_mem[c], out_mem[c], st->prefilter_period_old, st->prefilter_period, mode->shortMdctSize,
1928                st->prefilter_gain_old, st->prefilter_gain, st->prefilter_tapset_old, st->prefilter_tapset,
1929                mode->window, st->overlap);
1930          if (LM!=0)
1931             comb_filter(out_mem[c]+mode->shortMdctSize, out_mem[c]+mode->shortMdctSize, st->prefilter_period, pitch_index, N-mode->shortMdctSize,
1932                   st->prefilter_gain, gain1, st->prefilter_tapset, prefilter_tapset,
1933                   mode->window, overlap);
1934       } while (++c<CC);
1935
1936       /* We reuse freq[] as scratch space for the de-emphasis */
1937       deemphasis(out_mem, (opus_val16*)pcm, N, CC, st->upsample, mode->preemph, st->preemph_memD, freq);
1938       st->prefilter_period_old = st->prefilter_period;
1939       st->prefilter_gain_old = st->prefilter_gain;
1940       st->prefilter_tapset_old = st->prefilter_tapset;
1941    }
1942 #endif
1943
1944    st->prefilter_period = pitch_index;
1945    st->prefilter_gain = gain1;
1946    st->prefilter_tapset = prefilter_tapset;
1947 #ifdef RESYNTH
1948    if (LM!=0)
1949    {
1950       st->prefilter_period_old = st->prefilter_period;
1951       st->prefilter_gain_old = st->prefilter_gain;
1952       st->prefilter_tapset_old = st->prefilter_tapset;
1953    }
1954 #endif
1955
1956    if (CC==2&&C==1) {
1957       for (i=0;i<nbEBands;i++)
1958          oldBandE[nbEBands+i]=oldBandE[i];
1959    }
1960
1961    if (!isTransient)
1962    {
1963       for (i=0;i<CC*nbEBands;i++)
1964          oldLogE2[i] = oldLogE[i];
1965       for (i=0;i<CC*nbEBands;i++)
1966          oldLogE[i] = oldBandE[i];
1967    } else {
1968       for (i=0;i<CC*nbEBands;i++)
1969          oldLogE[i] = MIN16(oldLogE[i], oldBandE[i]);
1970    }
1971    /* In case start or end were to change */
1972    c=0; do
1973    {
1974       for (i=0;i<st->start;i++)
1975       {
1976          oldBandE[c*nbEBands+i]=0;
1977          oldLogE[c*nbEBands+i]=oldLogE2[c*nbEBands+i]=-QCONST16(28.f,DB_SHIFT);
1978       }
1979       for (i=st->end;i<nbEBands;i++)
1980       {
1981          oldBandE[c*nbEBands+i]=0;
1982          oldLogE[c*nbEBands+i]=oldLogE2[c*nbEBands+i]=-QCONST16(28.f,DB_SHIFT);
1983       }
1984    } while (++c<CC);
1985
1986    if (isTransient || transient_got_disabled)
1987       st->consec_transient++;
1988    else
1989       st->consec_transient=0;
1990    st->rng = enc->rng;
1991
1992    /* If there's any room left (can only happen for very high rates),
1993       it's already filled with zeros */
1994    ec_enc_done(enc);
1995
1996 #ifdef CUSTOM_MODES
1997    if (st->signalling)
1998       nbCompressedBytes++;
1999 #endif
2000
2001    RESTORE_STACK;
2002    if (ec_get_error(enc))
2003       return OPUS_INTERNAL_ERROR;
2004    else
2005       return nbCompressedBytes;
2006 }
2007
2008
2009 #ifdef CUSTOM_MODES
2010
2011 #ifdef FIXED_POINT
2012 int opus_custom_encode(CELTEncoder * OPUS_RESTRICT st, const opus_int16 * pcm, int frame_size, unsigned char *compressed, int nbCompressedBytes)
2013 {
2014    return celt_encode_with_ec(st, pcm, frame_size, compressed, nbCompressedBytes, NULL);
2015 }
2016
2017 #ifndef DISABLE_FLOAT_API
2018 int opus_custom_encode_float(CELTEncoder * OPUS_RESTRICT st, const float * pcm, int frame_size, unsigned char *compressed, int nbCompressedBytes)
2019 {
2020    int j, ret, C, N;
2021    VARDECL(opus_int16, in);
2022    ALLOC_STACK;
2023
2024    if (pcm==NULL)
2025       return OPUS_BAD_ARG;
2026
2027    C = st->channels;
2028    N = frame_size;
2029    ALLOC(in, C*N, opus_int16);
2030
2031    for (j=0;j<C*N;j++)
2032      in[j] = FLOAT2INT16(pcm[j]);
2033
2034    ret=celt_encode_with_ec(st,in,frame_size,compressed,nbCompressedBytes, NULL);
2035 #ifdef RESYNTH
2036    for (j=0;j<C*N;j++)
2037       ((float*)pcm)[j]=in[j]*(1.f/32768.f);
2038 #endif
2039    RESTORE_STACK;
2040    return ret;
2041 }
2042 #endif /* DISABLE_FLOAT_API */
2043 #else
2044
2045 int opus_custom_encode(CELTEncoder * OPUS_RESTRICT st, const opus_int16 * pcm, int frame_size, unsigned char *compressed, int nbCompressedBytes)
2046 {
2047    int j, ret, C, N;
2048    VARDECL(celt_sig, in);
2049    ALLOC_STACK;
2050
2051    if (pcm==NULL)
2052       return OPUS_BAD_ARG;
2053
2054    C=st->channels;
2055    N=frame_size;
2056    ALLOC(in, C*N, celt_sig);
2057    for (j=0;j<C*N;j++) {
2058      in[j] = SCALEOUT(pcm[j]);
2059    }
2060
2061    ret = celt_encode_with_ec(st,in,frame_size,compressed,nbCompressedBytes, NULL);
2062 #ifdef RESYNTH
2063    for (j=0;j<C*N;j++)
2064       ((opus_int16*)pcm)[j] = FLOAT2INT16(in[j]);
2065 #endif
2066    RESTORE_STACK;
2067    return ret;
2068 }
2069
2070 int opus_custom_encode_float(CELTEncoder * OPUS_RESTRICT st, const float * pcm, int frame_size, unsigned char *compressed, int nbCompressedBytes)
2071 {
2072    return celt_encode_with_ec(st, pcm, frame_size, compressed, nbCompressedBytes, NULL);
2073 }
2074
2075 #endif
2076
2077 #endif /* CUSTOM_MODES */
2078
2079 int opus_custom_encoder_ctl(CELTEncoder * OPUS_RESTRICT st, int request, ...)
2080 {
2081    va_list ap;
2082
2083    va_start(ap, request);
2084    switch (request)
2085    {
2086       case OPUS_SET_COMPLEXITY_REQUEST:
2087       {
2088          int value = va_arg(ap, opus_int32);
2089          if (value<0 || value>10)
2090             goto bad_arg;
2091          st->complexity = value;
2092       }
2093       break;
2094       case CELT_SET_START_BAND_REQUEST:
2095       {
2096          opus_int32 value = va_arg(ap, opus_int32);
2097          if (value<0 || value>=st->mode->nbEBands)
2098             goto bad_arg;
2099          st->start = value;
2100       }
2101       break;
2102       case CELT_SET_END_BAND_REQUEST:
2103       {
2104          opus_int32 value = va_arg(ap, opus_int32);
2105          if (value<1 || value>st->mode->nbEBands)
2106             goto bad_arg;
2107          st->end = value;
2108       }
2109       break;
2110       case CELT_SET_PREDICTION_REQUEST:
2111       {
2112          int value = va_arg(ap, opus_int32);
2113          if (value<0 || value>2)
2114             goto bad_arg;
2115          st->disable_pf = value<=1;
2116          st->force_intra = value==0;
2117       }
2118       break;
2119       case OPUS_SET_PACKET_LOSS_PERC_REQUEST:
2120       {
2121          int value = va_arg(ap, opus_int32);
2122          if (value<0 || value>100)
2123             goto bad_arg;
2124          st->loss_rate = value;
2125       }
2126       break;
2127       case OPUS_SET_VBR_CONSTRAINT_REQUEST:
2128       {
2129          opus_int32 value = va_arg(ap, opus_int32);
2130          st->constrained_vbr = value;
2131       }
2132       break;
2133       case OPUS_SET_VBR_REQUEST:
2134       {
2135          opus_int32 value = va_arg(ap, opus_int32);
2136          st->vbr = value;
2137       }
2138       break;
2139       case OPUS_SET_BITRATE_REQUEST:
2140       {
2141          opus_int32 value = va_arg(ap, opus_int32);
2142          if (value<=500 && value!=OPUS_BITRATE_MAX)
2143             goto bad_arg;
2144          value = IMIN(value, 260000*st->channels);
2145          st->bitrate = value;
2146       }
2147       break;
2148       case CELT_SET_CHANNELS_REQUEST:
2149       {
2150          opus_int32 value = va_arg(ap, opus_int32);
2151          if (value<1 || value>2)
2152             goto bad_arg;
2153          st->stream_channels = value;
2154       }
2155       break;
2156       case OPUS_SET_LSB_DEPTH_REQUEST:
2157       {
2158           opus_int32 value = va_arg(ap, opus_int32);
2159           if (value<8 || value>24)
2160              goto bad_arg;
2161           st->lsb_depth=value;
2162       }
2163       break;
2164       case OPUS_GET_LSB_DEPTH_REQUEST:
2165       {
2166           opus_int32 *value = va_arg(ap, opus_int32*);
2167           *value=st->lsb_depth;
2168       }
2169       break;
2170       case OPUS_SET_EXPERT_FRAME_DURATION_REQUEST:
2171       {
2172           opus_int32 value = va_arg(ap, opus_int32);
2173           st->variable_duration = value;
2174       }
2175       break;
2176       case OPUS_RESET_STATE:
2177       {
2178          int i;
2179          opus_val16 *oldBandE, *oldLogE, *oldLogE2;
2180          oldBandE = (opus_val16*)(st->in_mem+st->channels*(st->overlap+COMBFILTER_MAXPERIOD));
2181          oldLogE = oldBandE + st->channels*st->mode->nbEBands;
2182          oldLogE2 = oldLogE + st->channels*st->mode->nbEBands;
2183          OPUS_CLEAR((char*)&st->ENCODER_RESET_START,
2184                opus_custom_encoder_get_size(st->mode, st->channels)-
2185                ((char*)&st->ENCODER_RESET_START - (char*)st));
2186          for (i=0;i<st->channels*st->mode->nbEBands;i++)
2187             oldLogE[i]=oldLogE2[i]=-QCONST16(28.f,DB_SHIFT);
2188          st->vbr_offset = 0;
2189          st->delayedIntra = 1;
2190          st->spread_decision = SPREAD_NORMAL;
2191          st->tonal_average = 256;
2192          st->hf_average = 0;
2193          st->tapset_decision = 0;
2194       }
2195       break;
2196 #ifdef CUSTOM_MODES
2197       case CELT_SET_INPUT_CLIPPING_REQUEST:
2198       {
2199          opus_int32 value = va_arg(ap, opus_int32);
2200          st->clip = value;
2201       }
2202       break;
2203 #endif
2204       case CELT_SET_SIGNALLING_REQUEST:
2205       {
2206          opus_int32 value = va_arg(ap, opus_int32);
2207          st->signalling = value;
2208       }
2209       break;
2210       case CELT_SET_ANALYSIS_REQUEST:
2211       {
2212          AnalysisInfo *info = va_arg(ap, AnalysisInfo *);
2213          if (info)
2214             OPUS_COPY(&st->analysis, info, 1);
2215       }
2216       break;
2217       case CELT_GET_MODE_REQUEST:
2218       {
2219          const CELTMode ** value = va_arg(ap, const CELTMode**);
2220          if (value==0)
2221             goto bad_arg;
2222          *value=st->mode;
2223       }
2224       break;
2225       case OPUS_GET_FINAL_RANGE_REQUEST:
2226       {
2227          opus_uint32 * value = va_arg(ap, opus_uint32 *);
2228          if (value==0)
2229             goto bad_arg;
2230          *value=st->rng;
2231       }
2232       break;
2233       case OPUS_SET_LFE_REQUEST:
2234       {
2235           opus_int32 value = va_arg(ap, opus_int32);
2236           st->lfe = value;
2237       }
2238       break;
2239       case OPUS_SET_ENERGY_SAVE_REQUEST:
2240       {
2241           opus_val16 *value = va_arg(ap, opus_val16*);
2242           st->energy_save=value;
2243       }
2244       break;
2245       case OPUS_SET_ENERGY_MASK_REQUEST:
2246       {
2247           opus_val16 *value = va_arg(ap, opus_val16*);
2248           st->energy_mask = value;
2249       }
2250       break;
2251       default:
2252          goto bad_request;
2253    }
2254    va_end(ap);
2255    return OPUS_OK;
2256 bad_arg:
2257    va_end(ap);
2258    return OPUS_BAD_ARG;
2259 bad_request:
2260    va_end(ap);
2261    return OPUS_UNIMPLEMENTED;
2262 }