Surround masking rewrite
[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_mask;
115    opus_val16 spec_avg;
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 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(31.9f, 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             opus_int32 cap = ((effectiveBytes/4)<<BITRES<<3);
977             offsets[i] = cap-tot_boost;
978             tot_boost = cap;
979             break;
980          } else {
981             offsets[i] = boost;
982             tot_boost += boost_bits;
983          }
984       }
985    }
986    *tot_boost_ = tot_boost;
987    RESTORE_STACK;
988    return maxDepth;
989 }
990
991
992 static int run_prefilter(CELTEncoder *st, celt_sig *in, celt_sig *prefilter_mem, int CC, int N,
993       int prefilter_tapset, int *pitch, opus_val16 *gain, int *qgain, int enabled, int nbAvailableBytes)
994 {
995    int c;
996    VARDECL(celt_sig, _pre);
997    celt_sig *pre[2];
998    const CELTMode *mode;
999    int pitch_index;
1000    opus_val16 gain1;
1001    opus_val16 pf_threshold;
1002    int pf_on;
1003    int qg;
1004    SAVE_STACK;
1005
1006    mode = st->mode;
1007    ALLOC(_pre, CC*(N+COMBFILTER_MAXPERIOD), celt_sig);
1008
1009    pre[0] = _pre;
1010    pre[1] = _pre + (N+COMBFILTER_MAXPERIOD);
1011
1012
1013    c=0; do {
1014       OPUS_COPY(pre[c], prefilter_mem+c*COMBFILTER_MAXPERIOD, COMBFILTER_MAXPERIOD);
1015       OPUS_COPY(pre[c]+COMBFILTER_MAXPERIOD, in+c*(N+st->overlap)+st->overlap, N);
1016    } while (++c<CC);
1017
1018    if (enabled)
1019    {
1020       VARDECL(opus_val16, pitch_buf);
1021       ALLOC(pitch_buf, (COMBFILTER_MAXPERIOD+N)>>1, opus_val16);
1022
1023       pitch_downsample(pre, pitch_buf, COMBFILTER_MAXPERIOD+N, CC);
1024       /* Don't search for the fir last 1.5 octave of the range because
1025          there's too many false-positives due to short-term correlation */
1026       pitch_search(pitch_buf+(COMBFILTER_MAXPERIOD>>1), pitch_buf, N,
1027             COMBFILTER_MAXPERIOD-3*COMBFILTER_MINPERIOD, &pitch_index);
1028       pitch_index = COMBFILTER_MAXPERIOD-pitch_index;
1029
1030       gain1 = remove_doubling(pitch_buf, COMBFILTER_MAXPERIOD, COMBFILTER_MINPERIOD,
1031             N, &pitch_index, st->prefilter_period, st->prefilter_gain);
1032       if (pitch_index > COMBFILTER_MAXPERIOD-2)
1033          pitch_index = COMBFILTER_MAXPERIOD-2;
1034       gain1 = MULT16_16_Q15(QCONST16(.7f,15),gain1);
1035       /*printf("%d %d %f %f\n", pitch_change, pitch_index, gain1, st->analysis.tonality);*/
1036       if (st->loss_rate>2)
1037          gain1 = HALF32(gain1);
1038       if (st->loss_rate>4)
1039          gain1 = HALF32(gain1);
1040       if (st->loss_rate>8)
1041          gain1 = 0;
1042    } else {
1043       gain1 = 0;
1044       pitch_index = COMBFILTER_MINPERIOD;
1045    }
1046
1047    /* Gain threshold for enabling the prefilter/postfilter */
1048    pf_threshold = QCONST16(.2f,15);
1049
1050    /* Adjusting the threshold based on rate and continuity */
1051    if (abs(pitch_index-st->prefilter_period)*10>pitch_index)
1052       pf_threshold += QCONST16(.2f,15);
1053    if (nbAvailableBytes<25)
1054       pf_threshold += QCONST16(.1f,15);
1055    if (nbAvailableBytes<35)
1056       pf_threshold += QCONST16(.1f,15);
1057    if (st->prefilter_gain > QCONST16(.4f,15))
1058       pf_threshold -= QCONST16(.1f,15);
1059    if (st->prefilter_gain > QCONST16(.55f,15))
1060       pf_threshold -= QCONST16(.1f,15);
1061
1062    /* Hard threshold at 0.2 */
1063    pf_threshold = MAX16(pf_threshold, QCONST16(.2f,15));
1064    if (gain1<pf_threshold)
1065    {
1066       gain1 = 0;
1067       pf_on = 0;
1068       qg = 0;
1069    } else {
1070       /*This block is not gated by a total bits check only because
1071         of the nbAvailableBytes check above.*/
1072       if (ABS16(gain1-st->prefilter_gain)<QCONST16(.1f,15))
1073          gain1=st->prefilter_gain;
1074
1075 #ifdef FIXED_POINT
1076       qg = ((gain1+1536)>>10)/3-1;
1077 #else
1078       qg = (int)floor(.5f+gain1*32/3)-1;
1079 #endif
1080       qg = IMAX(0, IMIN(7, qg));
1081       gain1 = QCONST16(0.09375f,15)*(qg+1);
1082       pf_on = 1;
1083    }
1084    /*printf("%d %f\n", pitch_index, gain1);*/
1085
1086    c=0; do {
1087       int offset = mode->shortMdctSize-st->overlap;
1088       st->prefilter_period=IMAX(st->prefilter_period, COMBFILTER_MINPERIOD);
1089       OPUS_COPY(in+c*(N+st->overlap), st->in_mem+c*(st->overlap), st->overlap);
1090       if (offset)
1091          comb_filter(in+c*(N+st->overlap)+st->overlap, pre[c]+COMBFILTER_MAXPERIOD,
1092                st->prefilter_period, st->prefilter_period, offset, -st->prefilter_gain, -st->prefilter_gain,
1093                st->prefilter_tapset, st->prefilter_tapset, NULL, 0);
1094
1095       comb_filter(in+c*(N+st->overlap)+st->overlap+offset, pre[c]+COMBFILTER_MAXPERIOD+offset,
1096             st->prefilter_period, pitch_index, N-offset, -st->prefilter_gain, -gain1,
1097             st->prefilter_tapset, prefilter_tapset, mode->window, st->overlap);
1098       OPUS_COPY(st->in_mem+c*(st->overlap), in+c*(N+st->overlap)+N, st->overlap);
1099
1100       if (N>COMBFILTER_MAXPERIOD)
1101       {
1102          OPUS_MOVE(prefilter_mem+c*COMBFILTER_MAXPERIOD, pre[c]+N, COMBFILTER_MAXPERIOD);
1103       } else {
1104          OPUS_MOVE(prefilter_mem+c*COMBFILTER_MAXPERIOD, prefilter_mem+c*COMBFILTER_MAXPERIOD+N, COMBFILTER_MAXPERIOD-N);
1105          OPUS_MOVE(prefilter_mem+c*COMBFILTER_MAXPERIOD+COMBFILTER_MAXPERIOD-N, pre[c]+COMBFILTER_MAXPERIOD, N);
1106       }
1107    } while (++c<CC);
1108
1109    RESTORE_STACK;
1110    *gain = gain1;
1111    *pitch = pitch_index;
1112    *qgain = qg;
1113    return pf_on;
1114 }
1115
1116 static int compute_vbr(const CELTMode *mode, AnalysisInfo *analysis, opus_int32 base_target,
1117       int LM, opus_int32 bitrate, int lastCodedBands, int C, int intensity,
1118       int constrained_vbr, opus_val16 stereo_saving, int tot_boost,
1119       opus_val16 tf_estimate, int pitch_change, opus_val16 maxDepth,
1120       int variable_duration, int lfe, int has_surround_mask, opus_val16 surround_masking,
1121       opus_val16 temporal_vbr)
1122 {
1123    /* The target rate in 8th bits per frame */
1124    opus_int32 target;
1125    int coded_bins;
1126    int coded_bands;
1127    opus_val16 tf_calibration;
1128    int nbEBands;
1129    const opus_int16 *eBands;
1130
1131    nbEBands = mode->nbEBands;
1132    eBands = mode->eBands;
1133
1134    coded_bands = lastCodedBands ? lastCodedBands : nbEBands;
1135    coded_bins = eBands[coded_bands]<<LM;
1136    if (C==2)
1137       coded_bins += eBands[IMIN(intensity, coded_bands)]<<LM;
1138
1139    target = base_target;
1140
1141    /*printf("%f %f %f %f %d %d ", st->analysis.activity, st->analysis.tonality, tf_estimate, st->stereo_saving, tot_boost, coded_bands);*/
1142 #ifndef FIXED_POINT
1143    if (analysis->valid && analysis->activity<.4)
1144       target -= (opus_int32)((coded_bins<<BITRES)*(.4f-analysis->activity));
1145 #endif
1146    /* Stereo savings */
1147    if (C==2)
1148    {
1149       int coded_stereo_bands;
1150       int coded_stereo_dof;
1151       opus_val16 max_frac;
1152       coded_stereo_bands = IMIN(intensity, coded_bands);
1153       coded_stereo_dof = (eBands[coded_stereo_bands]<<LM)-coded_stereo_bands;
1154       /* Maximum fraction of the bits we can save if the signal is mono. */
1155       max_frac = DIV32_16(MULT16_16(QCONST16(0.8f, 15), coded_stereo_dof), coded_bins);
1156       /*printf("%d %d %d ", coded_stereo_dof, coded_bins, tot_boost);*/
1157       target -= (opus_int32)MIN32(MULT16_32_Q15(max_frac,target),
1158                       SHR32(MULT16_16(stereo_saving-QCONST16(0.1f,8),(coded_stereo_dof<<BITRES)),8));
1159    }
1160    /* Boost the rate according to dynalloc (minus the dynalloc average for calibration). */
1161    target += tot_boost-(16<<LM);
1162    /* Apply transient boost, compensating for average boost. */
1163    tf_calibration = variable_duration==OPUS_FRAMESIZE_VARIABLE ?
1164                     QCONST16(0.02f,14) : QCONST16(0.04f,14);
1165    target += (opus_int32)SHL32(MULT16_32_Q15(tf_estimate-tf_calibration, target),1);
1166
1167 #ifndef FIXED_POINT
1168    /* Apply tonality boost */
1169    if (analysis->valid && !lfe)
1170    {
1171       opus_int32 tonal_target;
1172       float tonal;
1173
1174       /* Tonality boost (compensating for the average). */
1175       tonal = MAX16(0.f,analysis->tonality-.15f)-0.09f;
1176       tonal_target = target + (opus_int32)((coded_bins<<BITRES)*1.2f*tonal);
1177       if (pitch_change)
1178          tonal_target +=  (opus_int32)((coded_bins<<BITRES)*.8f);
1179       /*printf("%f %f ", st->analysis.tonality, tonal);*/
1180       target = tonal_target;
1181    }
1182 #endif
1183
1184    if (has_surround_mask&&!lfe)
1185    {
1186       opus_int32 surround_target = target + (opus_int32)SHR32(MULT16_16(surround_masking,coded_bins<<BITRES), DB_SHIFT);
1187       /*printf("%f %d %d %d %d %d %d ", surround_masking, coded_bins, st->end, st->intensity, surround_target, target, st->bitrate);*/
1188       target = IMAX(target/4, surround_target);
1189    }
1190
1191    {
1192       opus_int32 floor_depth;
1193       int bins;
1194       bins = eBands[nbEBands-2]<<LM;
1195       /*floor_depth = SHR32(MULT16_16((C*bins<<BITRES),celt_log2(SHL32(MAX16(1,sample_max),13))), DB_SHIFT);*/
1196       floor_depth = (opus_int32)SHR32(MULT16_16((C*bins<<BITRES),maxDepth), DB_SHIFT);
1197       floor_depth = IMAX(floor_depth, target>>2);
1198       target = IMIN(target, floor_depth);
1199       /*printf("%f %d\n", maxDepth, floor_depth);*/
1200    }
1201
1202    if ((!has_surround_mask||lfe) && (constrained_vbr || bitrate<64000))
1203    {
1204       opus_val16 rate_factor;
1205 #ifdef FIXED_POINT
1206       rate_factor = MAX16(0,(bitrate-32000));
1207 #else
1208       rate_factor = MAX16(0,(1.f/32768)*(bitrate-32000));
1209 #endif
1210       if (constrained_vbr)
1211          rate_factor = MIN16(rate_factor, QCONST16(0.67f, 15));
1212       target = base_target + (opus_int32)MULT16_32_Q15(rate_factor, target-base_target);
1213
1214    }
1215
1216    if (!has_surround_mask && tf_estimate < QCONST16(.2f, 14))
1217    {
1218       opus_val16 amount;
1219       opus_val16 tvbr_factor;
1220       amount = MULT16_16_Q15(QCONST16(.0000031f, 30), IMAX(0, IMIN(32000, 96000-bitrate)));
1221       tvbr_factor = SHR32(MULT16_16(temporal_vbr, amount), DB_SHIFT);
1222       target += (opus_int32)MULT16_32_Q15(tvbr_factor, target);
1223    }
1224
1225    /* Don't allow more than doubling the rate */
1226    target = IMIN(2*base_target, target);
1227
1228    return target;
1229 }
1230
1231 int celt_encode_with_ec(CELTEncoder * OPUS_RESTRICT st, const opus_val16 * pcm, int frame_size, unsigned char *compressed, int nbCompressedBytes, ec_enc *enc)
1232 {
1233    int i, c, N;
1234    opus_int32 bits;
1235    ec_enc _enc;
1236    VARDECL(celt_sig, in);
1237    VARDECL(celt_sig, freq);
1238    VARDECL(celt_norm, X);
1239    VARDECL(celt_ener, bandE);
1240    VARDECL(opus_val16, bandLogE);
1241    VARDECL(opus_val16, bandLogE2);
1242    VARDECL(int, fine_quant);
1243    VARDECL(opus_val16, error);
1244    VARDECL(int, pulses);
1245    VARDECL(int, cap);
1246    VARDECL(int, offsets);
1247    VARDECL(int, fine_priority);
1248    VARDECL(int, tf_res);
1249    VARDECL(unsigned char, collapse_masks);
1250    celt_sig *prefilter_mem;
1251    opus_val16 *oldBandE, *oldLogE, *oldLogE2;
1252    int shortBlocks=0;
1253    int isTransient=0;
1254    const int CC = st->channels;
1255    const int C = st->stream_channels;
1256    int LM, M;
1257    int tf_select;
1258    int nbFilledBytes, nbAvailableBytes;
1259    int effEnd;
1260    int codedBands;
1261    int tf_sum;
1262    int alloc_trim;
1263    int pitch_index=COMBFILTER_MINPERIOD;
1264    opus_val16 gain1 = 0;
1265    int dual_stereo=0;
1266    int effectiveBytes;
1267    int dynalloc_logp;
1268    opus_int32 vbr_rate;
1269    opus_int32 total_bits;
1270    opus_int32 total_boost;
1271    opus_int32 balance;
1272    opus_int32 tell;
1273    int prefilter_tapset=0;
1274    int pf_on;
1275    int anti_collapse_rsv;
1276    int anti_collapse_on=0;
1277    int silence=0;
1278    int tf_chan = 0;
1279    opus_val16 tf_estimate;
1280    int pitch_change=0;
1281    opus_int32 tot_boost;
1282    opus_val32 sample_max;
1283    opus_val16 maxDepth;
1284    const OpusCustomMode *mode;
1285    int nbEBands;
1286    int overlap;
1287    const opus_int16 *eBands;
1288    int secondMdct;
1289    int signalBandwidth;
1290    int transient_got_disabled=0;
1291    opus_val16 surround_masking=0;
1292    opus_val16 temporal_vbr=0;
1293    ALLOC_STACK;
1294
1295    mode = st->mode;
1296    nbEBands = mode->nbEBands;
1297    overlap = mode->overlap;
1298    eBands = mode->eBands;
1299    tf_estimate = 0;
1300    if (nbCompressedBytes<2 || pcm==NULL)
1301      return OPUS_BAD_ARG;
1302
1303    frame_size *= st->upsample;
1304    for (LM=0;LM<=mode->maxLM;LM++)
1305       if (mode->shortMdctSize<<LM==frame_size)
1306          break;
1307    if (LM>mode->maxLM)
1308       return OPUS_BAD_ARG;
1309    M=1<<LM;
1310    N = M*mode->shortMdctSize;
1311
1312    prefilter_mem = st->in_mem+CC*(st->overlap);
1313    oldBandE = (opus_val16*)(st->in_mem+CC*(st->overlap+COMBFILTER_MAXPERIOD));
1314    oldLogE = oldBandE + CC*nbEBands;
1315    oldLogE2 = oldLogE + CC*nbEBands;
1316
1317    if (enc==NULL)
1318    {
1319       tell=1;
1320       nbFilledBytes=0;
1321    } else {
1322       tell=ec_tell(enc);
1323       nbFilledBytes=(tell+4)>>3;
1324    }
1325
1326 #ifdef CUSTOM_MODES
1327    if (st->signalling && enc==NULL)
1328    {
1329       int tmp = (mode->effEBands-st->end)>>1;
1330       st->end = IMAX(1, mode->effEBands-tmp);
1331       compressed[0] = tmp<<5;
1332       compressed[0] |= LM<<3;
1333       compressed[0] |= (C==2)<<2;
1334       /* Convert "standard mode" to Opus header */
1335       if (mode->Fs==48000 && mode->shortMdctSize==120)
1336       {
1337          int c0 = toOpus(compressed[0]);
1338          if (c0<0)
1339             return OPUS_BAD_ARG;
1340          compressed[0] = c0;
1341       }
1342       compressed++;
1343       nbCompressedBytes--;
1344    }
1345 #else
1346    celt_assert(st->signalling==0);
1347 #endif
1348
1349    /* Can't produce more than 1275 output bytes */
1350    nbCompressedBytes = IMIN(nbCompressedBytes,1275);
1351    nbAvailableBytes = nbCompressedBytes - nbFilledBytes;
1352
1353    if (st->vbr && st->bitrate!=OPUS_BITRATE_MAX)
1354    {
1355       opus_int32 den=mode->Fs>>BITRES;
1356       vbr_rate=(st->bitrate*frame_size+(den>>1))/den;
1357 #ifdef CUSTOM_MODES
1358       if (st->signalling)
1359          vbr_rate -= 8<<BITRES;
1360 #endif
1361       effectiveBytes = vbr_rate>>(3+BITRES);
1362    } else {
1363       opus_int32 tmp;
1364       vbr_rate = 0;
1365       tmp = st->bitrate*frame_size;
1366       if (tell>1)
1367          tmp += tell;
1368       if (st->bitrate!=OPUS_BITRATE_MAX)
1369          nbCompressedBytes = IMAX(2, IMIN(nbCompressedBytes,
1370                (tmp+4*mode->Fs)/(8*mode->Fs)-!!st->signalling));
1371       effectiveBytes = nbCompressedBytes;
1372    }
1373
1374    if (enc==NULL)
1375    {
1376       ec_enc_init(&_enc, compressed, nbCompressedBytes);
1377       enc = &_enc;
1378    }
1379
1380    if (vbr_rate>0)
1381    {
1382       /* Computes the max bit-rate allowed in VBR mode to avoid violating the
1383           target rate and buffering.
1384          We must do this up front so that bust-prevention logic triggers
1385           correctly if we don't have enough bits. */
1386       if (st->constrained_vbr)
1387       {
1388          opus_int32 vbr_bound;
1389          opus_int32 max_allowed;
1390          /* We could use any multiple of vbr_rate as bound (depending on the
1391              delay).
1392             This is clamped to ensure we use at least two bytes if the encoder
1393              was entirely empty, but to allow 0 in hybrid mode. */
1394          vbr_bound = vbr_rate;
1395          max_allowed = IMIN(IMAX(tell==1?2:0,
1396                (vbr_rate+vbr_bound-st->vbr_reservoir)>>(BITRES+3)),
1397                nbAvailableBytes);
1398          if(max_allowed < nbAvailableBytes)
1399          {
1400             nbCompressedBytes = nbFilledBytes+max_allowed;
1401             nbAvailableBytes = max_allowed;
1402             ec_enc_shrink(enc, nbCompressedBytes);
1403          }
1404       }
1405    }
1406    total_bits = nbCompressedBytes*8;
1407
1408    effEnd = st->end;
1409    if (effEnd > mode->effEBands)
1410       effEnd = mode->effEBands;
1411
1412    ALLOC(in, CC*(N+st->overlap), celt_sig);
1413
1414    sample_max=MAX32(st->overlap_max, celt_maxabs16(pcm, C*(N-overlap)/st->upsample));
1415    st->overlap_max=celt_maxabs16(pcm+C*(N-overlap)/st->upsample, C*overlap/st->upsample);
1416    sample_max=MAX32(sample_max, st->overlap_max);
1417 #ifdef FIXED_POINT
1418    silence = (sample_max==0);
1419 #else
1420    silence = (sample_max <= (opus_val16)1/(1<<st->lsb_depth));
1421 #endif
1422 #ifdef FUZZING
1423    if ((rand()&0x3F)==0)
1424       silence = 1;
1425 #endif
1426    if (tell==1)
1427       ec_enc_bit_logp(enc, silence, 15);
1428    else
1429       silence=0;
1430    if (silence)
1431    {
1432       /*In VBR mode there is no need to send more than the minimum. */
1433       if (vbr_rate>0)
1434       {
1435          effectiveBytes=nbCompressedBytes=IMIN(nbCompressedBytes, nbFilledBytes+2);
1436          total_bits=nbCompressedBytes*8;
1437          nbAvailableBytes=2;
1438          ec_enc_shrink(enc, nbCompressedBytes);
1439       }
1440       /* Pretend we've filled all the remaining bits with zeros
1441             (that's what the initialiser did anyway) */
1442       tell = nbCompressedBytes*8;
1443       enc->nbits_total+=tell-ec_tell(enc);
1444    }
1445    c=0; do {
1446       preemphasis(pcm+c, in+c*(N+st->overlap)+st->overlap, N, CC, st->upsample,
1447                   mode->preemph, st->preemph_memE+c, st->clip);
1448    } while (++c<CC);
1449
1450
1451
1452    /* Find pitch period and gain */
1453    {
1454       int enabled;
1455       int qg;
1456       enabled = (st->lfe || nbAvailableBytes>12*C) && st->start==0 && !silence && !st->disable_pf
1457             && st->complexity >= 5 && !(st->consec_transient && LM!=3 && st->variable_duration==OPUS_FRAMESIZE_VARIABLE);
1458
1459       prefilter_tapset = st->tapset_decision;
1460       pf_on = run_prefilter(st, in, prefilter_mem, CC, N, prefilter_tapset, &pitch_index, &gain1, &qg, enabled, nbAvailableBytes);
1461       if ((gain1 > QCONST16(.4f,15) || st->prefilter_gain > QCONST16(.4f,15)) && (!st->analysis.valid || st->analysis.tonality > .3)
1462             && (pitch_index > 1.26*st->prefilter_period || pitch_index < .79*st->prefilter_period))
1463          pitch_change = 1;
1464       if (pf_on==0)
1465       {
1466          if(st->start==0 && tell+16<=total_bits)
1467             ec_enc_bit_logp(enc, 0, 1);
1468       } else {
1469          /*This block is not gated by a total bits check only because
1470            of the nbAvailableBytes check above.*/
1471          int octave;
1472          ec_enc_bit_logp(enc, 1, 1);
1473          pitch_index += 1;
1474          octave = EC_ILOG(pitch_index)-5;
1475          ec_enc_uint(enc, octave, 6);
1476          ec_enc_bits(enc, pitch_index-(16<<octave), 4+octave);
1477          pitch_index -= 1;
1478          ec_enc_bits(enc, qg, 3);
1479          ec_enc_icdf(enc, prefilter_tapset, tapset_icdf, 2);
1480       }
1481    }
1482
1483    isTransient = 0;
1484    shortBlocks = 0;
1485    if (st->complexity >= 1 && !st->lfe)
1486    {
1487       isTransient = transient_analysis(in, N+st->overlap, CC,
1488             &tf_estimate, &tf_chan);
1489    }
1490    if (LM>0 && ec_tell(enc)+3<=total_bits)
1491    {
1492       if (isTransient)
1493          shortBlocks = M;
1494    } else {
1495       isTransient = 0;
1496       transient_got_disabled=1;
1497    }
1498
1499    ALLOC(freq, CC*N, celt_sig); /**< Interleaved signal MDCTs */
1500    ALLOC(bandE,nbEBands*CC, celt_ener);
1501    ALLOC(bandLogE,nbEBands*CC, opus_val16);
1502
1503    secondMdct = shortBlocks && st->complexity>=8;
1504    ALLOC(bandLogE2, C*nbEBands, opus_val16);
1505    if (secondMdct)
1506    {
1507       compute_mdcts(mode, 0, in, freq, C, CC, LM, st->upsample);
1508       compute_band_energies(mode, freq, bandE, effEnd, C, M);
1509       amp2Log2(mode, effEnd, st->end, bandE, bandLogE2, C);
1510       for (i=0;i<C*nbEBands;i++)
1511          bandLogE2[i] += HALF16(SHL16(LM, DB_SHIFT));
1512    }
1513
1514    compute_mdcts(mode, shortBlocks, in, freq, C, CC, LM, st->upsample);
1515    if (CC==2&&C==1)
1516       tf_chan = 0;
1517    compute_band_energies(mode, freq, bandE, effEnd, C, M);
1518
1519    if (st->lfe)
1520    {
1521       for (i=2;i<st->end;i++)
1522       {
1523          bandE[i] = IMIN(bandE[i], MULT16_32_Q15(QCONST16(1e-4f,15),bandE[0]));
1524          bandE[i] = MAX32(bandE[i], EPSILON);
1525       }
1526    }
1527    amp2Log2(mode, effEnd, st->end, bandE, bandLogE, C);
1528    /* This computes how much masking takes place between surround channels */
1529    if (st->energy_mask&&!st->lfe)
1530    {
1531       opus_val32 mask_avg=0;
1532       for (c=0;c<C;c++)
1533       {
1534          for(i=0;i<st->end;i++)
1535          {
1536             mask_avg += st->energy_mask[nbEBands*c+i];
1537          }
1538       }
1539       surround_masking = DIV32_16(mask_avg,C*st->end);
1540       surround_masking = MIN16(MAX16(surround_masking, -QCONST16(2.f, DB_SHIFT)), QCONST16(.2f, DB_SHIFT));
1541       surround_masking -= HALF16(HALF16(surround_masking));
1542    }
1543    /* Temporal VBR (but not for LFE) */
1544    if (!st->lfe)
1545    {
1546       opus_val16 follow=-QCONST16(10.0f,DB_SHIFT);
1547       float frame_avg=0;
1548       opus_val16 offset = shortBlocks?HALF16(SHL16(LM, DB_SHIFT)):0;
1549       for(i=st->start;i<st->end;i++)
1550       {
1551          follow = MAX16(follow-QCONST16(1.f, DB_SHIFT), bandLogE[i]-offset);
1552          if (C==2)
1553             follow = MAX16(follow, bandLogE[i+nbEBands]-offset);
1554          frame_avg += follow;
1555       }
1556       frame_avg /= (st->end-st->start);
1557       temporal_vbr = SUB16(frame_avg,st->spec_avg);
1558       temporal_vbr = MIN16(QCONST16(3.f, DB_SHIFT), MAX16(-QCONST16(1.5f, DB_SHIFT), temporal_vbr));
1559       st->spec_avg += MULT16_16_Q15(QCONST16(.02f, 15), temporal_vbr);
1560    }
1561    /*for (i=0;i<21;i++)
1562       printf("%f ", bandLogE[i]);
1563    printf("\n");*/
1564
1565    if (!secondMdct)
1566    {
1567       for (i=0;i<C*nbEBands;i++)
1568          bandLogE2[i] = bandLogE[i];
1569    }
1570
1571    /* Last chance to catch any transient we might have missed in the
1572       time-domain analysis */
1573    if (LM>0 && ec_tell(enc)+3<=total_bits && !isTransient && st->complexity>=5 && !st->lfe)
1574    {
1575       if (patch_transient_decision(bandLogE, oldBandE, nbEBands, st->end, C))
1576       {
1577          isTransient = 1;
1578          shortBlocks = M;
1579          compute_mdcts(mode, shortBlocks, in, freq, C, CC, LM, st->upsample);
1580          compute_band_energies(mode, freq, bandE, effEnd, C, M);
1581          amp2Log2(mode, effEnd, st->end, bandE, bandLogE, C);
1582          /* Compensate for the scaling of short vs long mdcts */
1583          for (i=0;i<C*nbEBands;i++)
1584             bandLogE2[i] += HALF16(SHL16(LM, DB_SHIFT));
1585          tf_estimate = QCONST16(.2f,14);
1586       }
1587    }
1588
1589    if (LM>0 && ec_tell(enc)+3<=total_bits)
1590       ec_enc_bit_logp(enc, isTransient, 3);
1591
1592    ALLOC(X, C*N, celt_norm);         /**< Interleaved normalised MDCTs */
1593
1594    /* Band normalisation */
1595    normalise_bands(mode, freq, X, bandE, effEnd, C, M);
1596
1597    ALLOC(tf_res, nbEBands, int);
1598    /* Disable variable tf resolution for hybrid and at very low bitrate */
1599    if (effectiveBytes>=15*C && st->start==0 && st->complexity>=2 && !st->lfe)
1600    {
1601       int lambda;
1602       if (effectiveBytes<40)
1603          lambda = 12;
1604       else if (effectiveBytes<60)
1605          lambda = 6;
1606       else if (effectiveBytes<100)
1607          lambda = 4;
1608       else
1609          lambda = 3;
1610       lambda*=2;
1611       tf_select = tf_analysis(mode, effEnd, isTransient, tf_res, lambda, X, N, LM, &tf_sum, tf_estimate, tf_chan);
1612       for (i=effEnd;i<st->end;i++)
1613          tf_res[i] = tf_res[effEnd-1];
1614    } else {
1615       tf_sum = 0;
1616       for (i=0;i<st->end;i++)
1617          tf_res[i] = isTransient;
1618       tf_select=0;
1619    }
1620
1621    ALLOC(error, C*nbEBands, opus_val16);
1622    quant_coarse_energy(mode, st->start, st->end, effEnd, bandLogE,
1623          oldBandE, total_bits, error, enc,
1624          C, LM, nbAvailableBytes, st->force_intra,
1625          &st->delayedIntra, st->complexity >= 4, st->loss_rate, st->lfe);
1626
1627    tf_encode(st->start, st->end, isTransient, tf_res, LM, tf_select, enc);
1628
1629    if (ec_tell(enc)+4<=total_bits)
1630    {
1631       if (st->lfe)
1632       {
1633          st->tapset_decision = 0;
1634          st->spread_decision = SPREAD_NORMAL;
1635       } else if (shortBlocks || st->complexity < 3 || nbAvailableBytes < 10*C || st->start != 0)
1636       {
1637          if (st->complexity == 0)
1638             st->spread_decision = SPREAD_NONE;
1639          else
1640             st->spread_decision = SPREAD_NORMAL;
1641       } else {
1642          /* Disable new spreading+tapset estimator until we can show it works
1643             better than the old one. So far it seems like spreading_decision()
1644             works best. */
1645          if (0&&st->analysis.valid)
1646          {
1647             static const opus_val16 spread_thresholds[3] = {-QCONST16(.6f, 15), -QCONST16(.2f, 15), -QCONST16(.07f, 15)};
1648             static const opus_val16 spread_histeresis[3] = {QCONST16(.15f, 15), QCONST16(.07f, 15), QCONST16(.02f, 15)};
1649             static const opus_val16 tapset_thresholds[2] = {QCONST16(.0f, 15), QCONST16(.15f, 15)};
1650             static const opus_val16 tapset_histeresis[2] = {QCONST16(.1f, 15), QCONST16(.05f, 15)};
1651             st->spread_decision = hysteresis_decision(-st->analysis.tonality, spread_thresholds, spread_histeresis, 3, st->spread_decision);
1652             st->tapset_decision = hysteresis_decision(st->analysis.tonality_slope, tapset_thresholds, tapset_histeresis, 2, st->tapset_decision);
1653          } else {
1654             st->spread_decision = spreading_decision(mode, X,
1655                   &st->tonal_average, st->spread_decision, &st->hf_average,
1656                   &st->tapset_decision, pf_on&&!shortBlocks, effEnd, C, M);
1657          }
1658          /*printf("%d %d\n", st->tapset_decision, st->spread_decision);*/
1659          /*printf("%f %d %f %d\n\n", st->analysis.tonality, st->spread_decision, st->analysis.tonality_slope, st->tapset_decision);*/
1660       }
1661       ec_enc_icdf(enc, st->spread_decision, spread_icdf, 5);
1662    }
1663
1664    ALLOC(offsets, nbEBands, int);
1665
1666    maxDepth = dynalloc_analysis(bandLogE, bandLogE2, nbEBands, st->start, st->end, C, offsets,
1667          st->lsb_depth, mode->logN, isTransient, st->vbr, st->constrained_vbr,
1668          eBands, LM, effectiveBytes, &tot_boost, st->lfe);
1669    /* For LFE, everything interesting is in the first band */
1670    if (st->lfe)
1671       offsets[0] = IMIN(8, effectiveBytes/3);
1672    ALLOC(cap, nbEBands, int);
1673    init_caps(mode,cap,LM,C);
1674
1675    dynalloc_logp = 6;
1676    total_bits<<=BITRES;
1677    total_boost = 0;
1678    tell = ec_tell_frac(enc);
1679    for (i=st->start;i<st->end;i++)
1680    {
1681       int width, quanta;
1682       int dynalloc_loop_logp;
1683       int boost;
1684       int j;
1685       width = C*(eBands[i+1]-eBands[i])<<LM;
1686       /* quanta is 6 bits, but no more than 1 bit/sample
1687          and no less than 1/8 bit/sample */
1688       quanta = IMIN(width<<BITRES, IMAX(6<<BITRES, width));
1689       dynalloc_loop_logp = dynalloc_logp;
1690       boost = 0;
1691       for (j = 0; tell+(dynalloc_loop_logp<<BITRES) < total_bits-total_boost
1692             && boost < cap[i]; j++)
1693       {
1694          int flag;
1695          flag = j<offsets[i];
1696          ec_enc_bit_logp(enc, flag, dynalloc_loop_logp);
1697          tell = ec_tell_frac(enc);
1698          if (!flag)
1699             break;
1700          boost += quanta;
1701          total_boost += quanta;
1702          dynalloc_loop_logp = 1;
1703       }
1704       /* Making dynalloc more likely */
1705       if (j)
1706          dynalloc_logp = IMAX(2, dynalloc_logp-1);
1707       offsets[i] = boost;
1708    }
1709
1710    if (C==2)
1711    {
1712       int effectiveRate;
1713
1714       static const opus_val16 intensity_thresholds[21]=
1715       /* 0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19  20  off*/
1716         { 16,21,23,25,27,29,31,33,35,38,42,46,50,54,58,63,68,75,84,102,130};
1717       static const opus_val16 intensity_histeresis[21]=
1718         {  2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 4, 5, 6,  8, 12};
1719
1720       /* Always use MS for 2.5 ms frames until we can do a better analysis */
1721       if (LM!=0)
1722          dual_stereo = stereo_analysis(mode, X, LM, N);
1723
1724       /* Account for coarse energy */
1725       effectiveRate = (8*effectiveBytes - 80)>>LM;
1726
1727       /* effectiveRate in kb/s */
1728       effectiveRate = 2*effectiveRate/5;
1729
1730       st->intensity = hysteresis_decision((opus_val16)effectiveRate, intensity_thresholds, intensity_histeresis, 21, st->intensity);
1731       st->intensity = IMIN(st->end,IMAX(st->start, st->intensity));
1732    }
1733
1734    alloc_trim = 5;
1735    if (tell+(6<<BITRES) <= total_bits - total_boost)
1736    {
1737       if (st->lfe)
1738          alloc_trim = 5;
1739       else
1740          alloc_trim = alloc_trim_analysis(mode, X, bandLogE,
1741             st->end, LM, C, N, &st->analysis, &st->stereo_saving, tf_estimate, st->intensity);
1742       ec_enc_icdf(enc, alloc_trim, trim_icdf, 7);
1743       tell = ec_tell_frac(enc);
1744    }
1745
1746    /* Variable bitrate */
1747    if (vbr_rate>0)
1748    {
1749      opus_val16 alpha;
1750      opus_int32 delta;
1751      /* The target rate in 8th bits per frame */
1752      opus_int32 target, base_target;
1753      opus_int32 min_allowed;
1754      int lm_diff = mode->maxLM - LM;
1755
1756      /* Don't attempt to use more than 510 kb/s, even for frames smaller than 20 ms.
1757         The CELT allocator will just not be able to use more than that anyway. */
1758      nbCompressedBytes = IMIN(nbCompressedBytes,1275>>(3-LM));
1759      base_target = vbr_rate - ((40*C+20)<<BITRES);
1760
1761      if (st->constrained_vbr)
1762         base_target += (st->vbr_offset>>lm_diff);
1763
1764      target = compute_vbr(mode, &st->analysis, base_target, LM, st->bitrate,
1765            st->lastCodedBands, C, st->intensity, st->constrained_vbr,
1766            st->stereo_saving, tot_boost, tf_estimate, pitch_change, maxDepth,
1767            st->variable_duration, st->lfe, st->energy_mask!=NULL, surround_masking,
1768            temporal_vbr);
1769
1770      /* The current offset is removed from the target and the space used
1771         so far is added*/
1772      target=target+tell;
1773      /* In VBR mode the frame size must not be reduced so much that it would
1774          result in the encoder running out of bits.
1775         The margin of 2 bytes ensures that none of the bust-prevention logic
1776          in the decoder will have triggered so far. */
1777      min_allowed = ((tell+total_boost+(1<<(BITRES+3))-1)>>(BITRES+3)) + 2 - nbFilledBytes;
1778
1779      nbAvailableBytes = (target+(1<<(BITRES+2)))>>(BITRES+3);
1780      nbAvailableBytes = IMAX(min_allowed,nbAvailableBytes);
1781      nbAvailableBytes = IMIN(nbCompressedBytes,nbAvailableBytes+nbFilledBytes) - nbFilledBytes;
1782
1783      /* By how much did we "miss" the target on that frame */
1784      delta = target - vbr_rate;
1785
1786      target=nbAvailableBytes<<(BITRES+3);
1787
1788      /*If the frame is silent we don't adjust our drift, otherwise
1789        the encoder will shoot to very high rates after hitting a
1790        span of silence, but we do allow the bitres to refill.
1791        This means that we'll undershoot our target in CVBR/VBR modes
1792        on files with lots of silence. */
1793      if(silence)
1794      {
1795        nbAvailableBytes = 2;
1796        target = 2*8<<BITRES;
1797        delta = 0;
1798      }
1799
1800      if (st->vbr_count < 970)
1801      {
1802         st->vbr_count++;
1803         alpha = celt_rcp(SHL32(EXTEND32(st->vbr_count+20),16));
1804      } else
1805         alpha = QCONST16(.001f,15);
1806      /* How many bits have we used in excess of what we're allowed */
1807      if (st->constrained_vbr)
1808         st->vbr_reservoir += target - vbr_rate;
1809      /*printf ("%d\n", st->vbr_reservoir);*/
1810
1811      /* Compute the offset we need to apply in order to reach the target */
1812      if (st->constrained_vbr)
1813      {
1814         st->vbr_drift += (opus_int32)MULT16_32_Q15(alpha,(delta*(1<<lm_diff))-st->vbr_offset-st->vbr_drift);
1815         st->vbr_offset = -st->vbr_drift;
1816      }
1817      /*printf ("%d\n", st->vbr_drift);*/
1818
1819      if (st->constrained_vbr && st->vbr_reservoir < 0)
1820      {
1821         /* We're under the min value -- increase rate */
1822         int adjust = (-st->vbr_reservoir)/(8<<BITRES);
1823         /* Unless we're just coding silence */
1824         nbAvailableBytes += silence?0:adjust;
1825         st->vbr_reservoir = 0;
1826         /*printf ("+%d\n", adjust);*/
1827      }
1828      nbCompressedBytes = IMIN(nbCompressedBytes,nbAvailableBytes+nbFilledBytes);
1829      /*printf("%d\n", nbCompressedBytes*50*8);*/
1830      /* This moves the raw bits to take into account the new compressed size */
1831      ec_enc_shrink(enc, nbCompressedBytes);
1832    }
1833
1834    /* Bit allocation */
1835    ALLOC(fine_quant, nbEBands, int);
1836    ALLOC(pulses, nbEBands, int);
1837    ALLOC(fine_priority, nbEBands, int);
1838
1839    /* bits =           packet size                    - where we are - safety*/
1840    bits = (((opus_int32)nbCompressedBytes*8)<<BITRES) - ec_tell_frac(enc) - 1;
1841    anti_collapse_rsv = isTransient&&LM>=2&&bits>=((LM+2)<<BITRES) ? (1<<BITRES) : 0;
1842    bits -= anti_collapse_rsv;
1843    signalBandwidth = st->end-1;
1844 #ifndef FIXED_POINT
1845    if (st->analysis.valid)
1846    {
1847       int min_bandwidth;
1848       if (st->bitrate < (opus_int32)32000*C)
1849          min_bandwidth = 13;
1850       else if (st->bitrate < (opus_int32)48000*C)
1851          min_bandwidth = 16;
1852       else if (st->bitrate < (opus_int32)60000*C)
1853          min_bandwidth = 18;
1854       else  if (st->bitrate < (opus_int32)80000*C)
1855          min_bandwidth = 19;
1856       else
1857          min_bandwidth = 20;
1858       signalBandwidth = IMAX(st->analysis.bandwidth, min_bandwidth);
1859    }
1860 #endif
1861    if (st->lfe)
1862       signalBandwidth = 1;
1863    codedBands = compute_allocation(mode, st->start, st->end, offsets, cap,
1864          alloc_trim, &st->intensity, &dual_stereo, bits, &balance, pulses,
1865          fine_quant, fine_priority, C, LM, enc, 1, st->lastCodedBands, signalBandwidth);
1866    if (st->lastCodedBands)
1867       st->lastCodedBands = IMIN(st->lastCodedBands+1,IMAX(st->lastCodedBands-1,codedBands));
1868    else
1869       st->lastCodedBands = codedBands;
1870
1871    quant_fine_energy(mode, st->start, st->end, oldBandE, error, fine_quant, enc, C);
1872
1873    /* Residual quantisation */
1874    ALLOC(collapse_masks, C*nbEBands, unsigned char);
1875    quant_all_bands(1, mode, st->start, st->end, X, C==2 ? X+N : NULL, collapse_masks,
1876          bandE, pulses, shortBlocks, st->spread_decision, dual_stereo, st->intensity, tf_res,
1877          nbCompressedBytes*(8<<BITRES)-anti_collapse_rsv, balance, enc, LM, codedBands, &st->rng);
1878
1879    if (anti_collapse_rsv > 0)
1880    {
1881       anti_collapse_on = st->consec_transient<2;
1882 #ifdef FUZZING
1883       anti_collapse_on = rand()&0x1;
1884 #endif
1885       ec_enc_bits(enc, anti_collapse_on, 1);
1886    }
1887    quant_energy_finalise(mode, st->start, st->end, oldBandE, error, fine_quant, fine_priority, nbCompressedBytes*8-ec_tell(enc), enc, C);
1888
1889    if (silence)
1890    {
1891       for (i=0;i<C*nbEBands;i++)
1892          oldBandE[i] = -QCONST16(28.f,DB_SHIFT);
1893    }
1894
1895 #ifdef RESYNTH
1896    /* Re-synthesis of the coded audio if required */
1897    {
1898       celt_sig *out_mem[2];
1899
1900       if (anti_collapse_on)
1901       {
1902          anti_collapse(mode, X, collapse_masks, LM, C, N,
1903                st->start, st->end, oldBandE, oldLogE, oldLogE2, pulses, st->rng);
1904       }
1905
1906       if (silence)
1907       {
1908          for (i=0;i<C*N;i++)
1909             freq[i] = 0;
1910       } else {
1911          /* Synthesis */
1912          denormalise_bands(mode, X, freq, oldBandE, st->start, effEnd, C, M);
1913       }
1914
1915       c=0; do {
1916          OPUS_MOVE(st->syn_mem[c], st->syn_mem[c]+N, 2*MAX_PERIOD-N+overlap/2);
1917       } while (++c<CC);
1918
1919       if (CC==2&&C==1)
1920       {
1921          for (i=0;i<N;i++)
1922             freq[N+i] = freq[i];
1923       }
1924
1925       c=0; do {
1926          out_mem[c] = st->syn_mem[c]+2*MAX_PERIOD-N;
1927       } while (++c<CC);
1928
1929       compute_inv_mdcts(mode, shortBlocks, freq, out_mem, CC, LM);
1930
1931       c=0; do {
1932          st->prefilter_period=IMAX(st->prefilter_period, COMBFILTER_MINPERIOD);
1933          st->prefilter_period_old=IMAX(st->prefilter_period_old, COMBFILTER_MINPERIOD);
1934          comb_filter(out_mem[c], out_mem[c], st->prefilter_period_old, st->prefilter_period, mode->shortMdctSize,
1935                st->prefilter_gain_old, st->prefilter_gain, st->prefilter_tapset_old, st->prefilter_tapset,
1936                mode->window, st->overlap);
1937          if (LM!=0)
1938             comb_filter(out_mem[c]+mode->shortMdctSize, out_mem[c]+mode->shortMdctSize, st->prefilter_period, pitch_index, N-mode->shortMdctSize,
1939                   st->prefilter_gain, gain1, st->prefilter_tapset, prefilter_tapset,
1940                   mode->window, overlap);
1941       } while (++c<CC);
1942
1943       /* We reuse freq[] as scratch space for the de-emphasis */
1944       deemphasis(out_mem, (opus_val16*)pcm, N, CC, st->upsample, mode->preemph, st->preemph_memD, freq);
1945       st->prefilter_period_old = st->prefilter_period;
1946       st->prefilter_gain_old = st->prefilter_gain;
1947       st->prefilter_tapset_old = st->prefilter_tapset;
1948    }
1949 #endif
1950
1951    st->prefilter_period = pitch_index;
1952    st->prefilter_gain = gain1;
1953    st->prefilter_tapset = prefilter_tapset;
1954 #ifdef RESYNTH
1955    if (LM!=0)
1956    {
1957       st->prefilter_period_old = st->prefilter_period;
1958       st->prefilter_gain_old = st->prefilter_gain;
1959       st->prefilter_tapset_old = st->prefilter_tapset;
1960    }
1961 #endif
1962
1963    if (CC==2&&C==1) {
1964       for (i=0;i<nbEBands;i++)
1965          oldBandE[nbEBands+i]=oldBandE[i];
1966    }
1967
1968    if (!isTransient)
1969    {
1970       for (i=0;i<CC*nbEBands;i++)
1971          oldLogE2[i] = oldLogE[i];
1972       for (i=0;i<CC*nbEBands;i++)
1973          oldLogE[i] = oldBandE[i];
1974    } else {
1975       for (i=0;i<CC*nbEBands;i++)
1976          oldLogE[i] = MIN16(oldLogE[i], oldBandE[i]);
1977    }
1978    /* In case start or end were to change */
1979    c=0; do
1980    {
1981       for (i=0;i<st->start;i++)
1982       {
1983          oldBandE[c*nbEBands+i]=0;
1984          oldLogE[c*nbEBands+i]=oldLogE2[c*nbEBands+i]=-QCONST16(28.f,DB_SHIFT);
1985       }
1986       for (i=st->end;i<nbEBands;i++)
1987       {
1988          oldBandE[c*nbEBands+i]=0;
1989          oldLogE[c*nbEBands+i]=oldLogE2[c*nbEBands+i]=-QCONST16(28.f,DB_SHIFT);
1990       }
1991    } while (++c<CC);
1992
1993    if (isTransient || transient_got_disabled)
1994       st->consec_transient++;
1995    else
1996       st->consec_transient=0;
1997    st->rng = enc->rng;
1998
1999    /* If there's any room left (can only happen for very high rates),
2000       it's already filled with zeros */
2001    ec_enc_done(enc);
2002
2003 #ifdef CUSTOM_MODES
2004    if (st->signalling)
2005       nbCompressedBytes++;
2006 #endif
2007
2008    RESTORE_STACK;
2009    if (ec_get_error(enc))
2010       return OPUS_INTERNAL_ERROR;
2011    else
2012       return nbCompressedBytes;
2013 }
2014
2015
2016 #ifdef CUSTOM_MODES
2017
2018 #ifdef FIXED_POINT
2019 int opus_custom_encode(CELTEncoder * OPUS_RESTRICT st, const opus_int16 * pcm, int frame_size, unsigned char *compressed, int nbCompressedBytes)
2020 {
2021    return celt_encode_with_ec(st, pcm, frame_size, compressed, nbCompressedBytes, NULL);
2022 }
2023
2024 #ifndef DISABLE_FLOAT_API
2025 int opus_custom_encode_float(CELTEncoder * OPUS_RESTRICT st, const float * pcm, int frame_size, unsigned char *compressed, int nbCompressedBytes)
2026 {
2027    int j, ret, C, N;
2028    VARDECL(opus_int16, in);
2029    ALLOC_STACK;
2030
2031    if (pcm==NULL)
2032       return OPUS_BAD_ARG;
2033
2034    C = st->channels;
2035    N = frame_size;
2036    ALLOC(in, C*N, opus_int16);
2037
2038    for (j=0;j<C*N;j++)
2039      in[j] = FLOAT2INT16(pcm[j]);
2040
2041    ret=celt_encode_with_ec(st,in,frame_size,compressed,nbCompressedBytes, NULL);
2042 #ifdef RESYNTH
2043    for (j=0;j<C*N;j++)
2044       ((float*)pcm)[j]=in[j]*(1.f/32768.f);
2045 #endif
2046    RESTORE_STACK;
2047    return ret;
2048 }
2049 #endif /* DISABLE_FLOAT_API */
2050 #else
2051
2052 int opus_custom_encode(CELTEncoder * OPUS_RESTRICT st, const opus_int16 * pcm, int frame_size, unsigned char *compressed, int nbCompressedBytes)
2053 {
2054    int j, ret, C, N;
2055    VARDECL(celt_sig, in);
2056    ALLOC_STACK;
2057
2058    if (pcm==NULL)
2059       return OPUS_BAD_ARG;
2060
2061    C=st->channels;
2062    N=frame_size;
2063    ALLOC(in, C*N, celt_sig);
2064    for (j=0;j<C*N;j++) {
2065      in[j] = SCALEOUT(pcm[j]);
2066    }
2067
2068    ret = celt_encode_with_ec(st,in,frame_size,compressed,nbCompressedBytes, NULL);
2069 #ifdef RESYNTH
2070    for (j=0;j<C*N;j++)
2071       ((opus_int16*)pcm)[j] = FLOAT2INT16(in[j]);
2072 #endif
2073    RESTORE_STACK;
2074    return ret;
2075 }
2076
2077 int opus_custom_encode_float(CELTEncoder * OPUS_RESTRICT st, const float * pcm, int frame_size, unsigned char *compressed, int nbCompressedBytes)
2078 {
2079    return celt_encode_with_ec(st, pcm, frame_size, compressed, nbCompressedBytes, NULL);
2080 }
2081
2082 #endif
2083
2084 #endif /* CUSTOM_MODES */
2085
2086 int opus_custom_encoder_ctl(CELTEncoder * OPUS_RESTRICT st, int request, ...)
2087 {
2088    va_list ap;
2089
2090    va_start(ap, request);
2091    switch (request)
2092    {
2093       case OPUS_SET_COMPLEXITY_REQUEST:
2094       {
2095          int value = va_arg(ap, opus_int32);
2096          if (value<0 || value>10)
2097             goto bad_arg;
2098          st->complexity = value;
2099       }
2100       break;
2101       case CELT_SET_START_BAND_REQUEST:
2102       {
2103          opus_int32 value = va_arg(ap, opus_int32);
2104          if (value<0 || value>=st->mode->nbEBands)
2105             goto bad_arg;
2106          st->start = value;
2107       }
2108       break;
2109       case CELT_SET_END_BAND_REQUEST:
2110       {
2111          opus_int32 value = va_arg(ap, opus_int32);
2112          if (value<1 || value>st->mode->nbEBands)
2113             goto bad_arg;
2114          st->end = value;
2115       }
2116       break;
2117       case CELT_SET_PREDICTION_REQUEST:
2118       {
2119          int value = va_arg(ap, opus_int32);
2120          if (value<0 || value>2)
2121             goto bad_arg;
2122          st->disable_pf = value<=1;
2123          st->force_intra = value==0;
2124       }
2125       break;
2126       case OPUS_SET_PACKET_LOSS_PERC_REQUEST:
2127       {
2128          int value = va_arg(ap, opus_int32);
2129          if (value<0 || value>100)
2130             goto bad_arg;
2131          st->loss_rate = value;
2132       }
2133       break;
2134       case OPUS_SET_VBR_CONSTRAINT_REQUEST:
2135       {
2136          opus_int32 value = va_arg(ap, opus_int32);
2137          st->constrained_vbr = value;
2138       }
2139       break;
2140       case OPUS_SET_VBR_REQUEST:
2141       {
2142          opus_int32 value = va_arg(ap, opus_int32);
2143          st->vbr = value;
2144       }
2145       break;
2146       case OPUS_SET_BITRATE_REQUEST:
2147       {
2148          opus_int32 value = va_arg(ap, opus_int32);
2149          if (value<=500 && value!=OPUS_BITRATE_MAX)
2150             goto bad_arg;
2151          value = IMIN(value, 260000*st->channels);
2152          st->bitrate = value;
2153       }
2154       break;
2155       case CELT_SET_CHANNELS_REQUEST:
2156       {
2157          opus_int32 value = va_arg(ap, opus_int32);
2158          if (value<1 || value>2)
2159             goto bad_arg;
2160          st->stream_channels = value;
2161       }
2162       break;
2163       case OPUS_SET_LSB_DEPTH_REQUEST:
2164       {
2165           opus_int32 value = va_arg(ap, opus_int32);
2166           if (value<8 || value>24)
2167              goto bad_arg;
2168           st->lsb_depth=value;
2169       }
2170       break;
2171       case OPUS_GET_LSB_DEPTH_REQUEST:
2172       {
2173           opus_int32 *value = va_arg(ap, opus_int32*);
2174           *value=st->lsb_depth;
2175       }
2176       break;
2177       case OPUS_SET_EXPERT_FRAME_DURATION_REQUEST:
2178       {
2179           opus_int32 value = va_arg(ap, opus_int32);
2180           st->variable_duration = value;
2181       }
2182       break;
2183       case OPUS_RESET_STATE:
2184       {
2185          int i;
2186          opus_val16 *oldBandE, *oldLogE, *oldLogE2;
2187          oldBandE = (opus_val16*)(st->in_mem+st->channels*(st->overlap+COMBFILTER_MAXPERIOD));
2188          oldLogE = oldBandE + st->channels*st->mode->nbEBands;
2189          oldLogE2 = oldLogE + st->channels*st->mode->nbEBands;
2190          OPUS_CLEAR((char*)&st->ENCODER_RESET_START,
2191                opus_custom_encoder_get_size(st->mode, st->channels)-
2192                ((char*)&st->ENCODER_RESET_START - (char*)st));
2193          for (i=0;i<st->channels*st->mode->nbEBands;i++)
2194             oldLogE[i]=oldLogE2[i]=-QCONST16(28.f,DB_SHIFT);
2195          st->vbr_offset = 0;
2196          st->delayedIntra = 1;
2197          st->spread_decision = SPREAD_NORMAL;
2198          st->tonal_average = 256;
2199          st->hf_average = 0;
2200          st->tapset_decision = 0;
2201       }
2202       break;
2203 #ifdef CUSTOM_MODES
2204       case CELT_SET_INPUT_CLIPPING_REQUEST:
2205       {
2206          opus_int32 value = va_arg(ap, opus_int32);
2207          st->clip = value;
2208       }
2209       break;
2210 #endif
2211       case CELT_SET_SIGNALLING_REQUEST:
2212       {
2213          opus_int32 value = va_arg(ap, opus_int32);
2214          st->signalling = value;
2215       }
2216       break;
2217       case CELT_SET_ANALYSIS_REQUEST:
2218       {
2219          AnalysisInfo *info = va_arg(ap, AnalysisInfo *);
2220          if (info)
2221             OPUS_COPY(&st->analysis, info, 1);
2222       }
2223       break;
2224       case CELT_GET_MODE_REQUEST:
2225       {
2226          const CELTMode ** value = va_arg(ap, const CELTMode**);
2227          if (value==0)
2228             goto bad_arg;
2229          *value=st->mode;
2230       }
2231       break;
2232       case OPUS_GET_FINAL_RANGE_REQUEST:
2233       {
2234          opus_uint32 * value = va_arg(ap, opus_uint32 *);
2235          if (value==0)
2236             goto bad_arg;
2237          *value=st->rng;
2238       }
2239       break;
2240       case OPUS_SET_LFE_REQUEST:
2241       {
2242           opus_int32 value = va_arg(ap, opus_int32);
2243           st->lfe = value;
2244       }
2245       break;
2246       case OPUS_SET_ENERGY_MASK_REQUEST:
2247       {
2248           opus_val16 *value = va_arg(ap, opus_val16*);
2249           st->energy_mask = value;
2250       }
2251       break;
2252       default:
2253          goto bad_request;
2254    }
2255    va_end(ap);
2256    return OPUS_OK;
2257 bad_arg:
2258    va_end(ap);
2259    return OPUS_BAD_ARG;
2260 bad_request:
2261    va_end(ap);
2262    return OPUS_UNIMPLEMENTED;
2263 }