/* Copyright (c) 2007-2008 CSIRO

   Copyright (c) 2007-2010 Xiph.Org Foundation

   Copyright (c) 2008 Gregory Maxwell

   Written by Jean-Marc Valin and Gregory Maxwell */

/*

   Redistribution and use in source and binary forms, with or without

   modification, are permitted provided that the following conditions

   are met:

   - Redistributions of source code must retain the above copyright

   notice, this list of conditions and the following disclaimer.

   - Redistributions in binary form must reproduce the above copyright

   notice, this list of conditions and the following disclaimer in the

   documentation and/or other materials provided with the distribution.

   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS

   ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT

   LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR

   A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER

   OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,

   EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,

   PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR

   PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF

   LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING

   NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS

   SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

*/

#ifdef HAVE_CONFIG_H

#include "config.h"

#endif

#define CELT_ENCODER_C

#include "cpu_support.h"

#include "os_support.h"

#include "mdct.h"

#include <math.h>

#include "celt.h"

#include "pitch.h"

#include "bands.h"

#include "modes.h"

#include "entcode.h"

#include "quant_bands.h"

#include "rate.h"

#include "stack_alloc.h"

#include "mathops.h"

#include "float_cast.h"

#include <stdarg.h>

#include "celt_lpc.h"

#include "vq.h"

#ifndef M_PI

#define M_PI 3.141592653

#endif

/** Encoder state

 @brief Encoder state

 */

struct OpusCustomEncoder {

   const OpusCustomMode *mode;     /**< Mode used by the encoder */

   int channels;

   int stream_channels;

   int force_intra;

   int clip;

   int disable_pf;

   int complexity;

   int upsample;

   int start, end;

   opus_int32 bitrate;

   int vbr;

   int signalling;

   int constrained_vbr;      /* If zero, VBR can do whatever it likes with the rate */

   int loss_rate;

   int lsb_depth;

   int lfe;

   int disable_inv;

   int arch;

#ifdef ENABLE_QEXT

   int enable_qext;

   int qext_scale;

#endif

   /* Everything beyond this point gets cleared on a reset */

#define ENCODER_RESET_START rng

   opus_uint32 rng;

   int spread_decision;

   opus_val32 delayedIntra;

   int tonal_average;

   int lastCodedBands;

   int hf_average;

   int tapset_decision;

   int prefilter_period;

   opus_val16 prefilter_gain;

   int prefilter_tapset;

#ifdef RESYNTH

   int prefilter_period_old;

   opus_val16 prefilter_gain_old;

   int prefilter_tapset_old;

#endif

   int consec_transient;

   AnalysisInfo analysis;

   SILKInfo silk_info;

   opus_val32 preemph_memE[2];

   opus_val32 preemph_memD[2];

   /* VBR-related parameters */

   opus_int32 vbr_reservoir;

   opus_int32 vbr_drift;

   opus_int32 vbr_offset;

   opus_int32 vbr_count;

   opus_val32 overlap_max;

   opus_val16 stereo_saving;

   int intensity;

   celt_glog *energy_mask;

   celt_glog spec_avg;

#ifdef RESYNTH

#ifdef ENABLE_QEXT

   /* +MAX_PERIOD/2 to make space for overlap */

   celt_sig syn_mem[2][2*DEC_PITCH_BUF_SIZE+MAX_PERIOD];

#else

   /* +MAX_PERIOD/2 to make space for overlap */

   celt_sig syn_mem[2][DEC_PITCH_BUF_SIZE+MAX_PERIOD/2];

#endif

#endif

   celt_sig in_mem[1]; /* Size = channels*mode->overlap */

   /* celt_sig prefilter_mem[],  Size = channels*COMBFILTER_MAXPERIOD */

   /* celt_glog oldBandE[],     Size = channels*mode->nbEBands */

   /* celt_glog oldLogE[],      Size = channels*mode->nbEBands */

   /* celt_glog oldLogE2[],     Size = channels*mode->nbEBands */

   /* celt_glog energyError[],  Size = channels*mode->nbEBands */

};

int celt_encoder_get_size(int channels)

{

#ifdef ENABLE_QEXT

   CELTMode *mode = opus_custom_mode_create(96000, 1920, NULL);

#else

   CELTMode *mode = opus_custom_mode_create(48000, 960, NULL);

#endif

   return opus_custom_encoder_get_size(mode, channels);

}

OPUS_CUSTOM_NOSTATIC int opus_custom_encoder_get_size(const CELTMode *mode, int channels)

{

   int extra=0;

   int size;

#ifdef ENABLE_QEXT

   int qext_scale;

   extra = channels*NB_QEXT_BANDS*sizeof(celt_glog);

   if (mode->Fs == 96000 && (mode->shortMdctSize==240 || mode->shortMdctSize==180)) {

      qext_scale = 2;

   } else qext_scale = 1;

#endif

   size = sizeof(struct CELTEncoder)

         + (channels*mode->overlap-1)*sizeof(celt_sig)    /* celt_sig in_mem[channels*mode->overlap]; */

         + channels*QEXT_SCALE(COMBFILTER_MAXPERIOD)*sizeof(celt_sig) /* celt_sig prefilter_mem[channels*COMBFILTER_MAXPERIOD]; */

         + 4*channels*mode->nbEBands*sizeof(celt_glog)    /* celt_glog oldBandE[channels*mode->nbEBands]; */

                                                          /* celt_glog oldLogE[channels*mode->nbEBands]; */

                                                          /* celt_glog oldLogE2[channels*mode->nbEBands]; */

                                                          /* celt_glog energyError[channels*mode->nbEBands]; */

         + extra;

   return size;

}

#if defined(CUSTOM_MODES) || defined(ENABLE_OPUS_CUSTOM_API)

CELTEncoder *opus_custom_encoder_create(const CELTMode *mode, int channels, int *error)

{

   int ret;

   CELTEncoder *st = (CELTEncoder *)opus_alloc(opus_custom_encoder_get_size(mode, channels));

   /* init will handle the NULL case */

   ret = opus_custom_encoder_init(st, mode, channels);

   if (ret != OPUS_OK)

   {

      opus_custom_encoder_destroy(st);

      st = NULL;

   }

   if (error)

      *error = ret;

   return st;

}

#endif /* CUSTOM_MODES */

static int opus_custom_encoder_init_arch(CELTEncoder *st, const CELTMode *mode,

                                         int channels, int arch)

{

   if (channels < 0 || channels > 2)

      return OPUS_BAD_ARG;

   if (st==NULL || mode==NULL)

      return OPUS_ALLOC_FAIL;

   OPUS_CLEAR((char*)st, opus_custom_encoder_get_size(mode, channels));

   st->mode = mode;

   st->stream_channels = st->channels = channels;

   st->upsample = 1;

   st->start = 0;

   st->end = st->mode->effEBands;

   st->signalling = 1;

   st->arch = arch;

   st->constrained_vbr = 1;

   st->clip = 1;

   st->bitrate = OPUS_BITRATE_MAX;

   st->vbr = 0;

   st->force_intra  = 0;

   st->complexity = 5;

   st->lsb_depth=24;

#ifdef ENABLE_QEXT

   if (st->mode->Fs == 96000 && (mode->shortMdctSize==240 || mode->shortMdctSize==180)) st->qext_scale = 2;

   else st->qext_scale = 1;

#endif

   opus_custom_encoder_ctl(st, OPUS_RESET_STATE);

   return OPUS_OK;

}

#if defined(CUSTOM_MODES) || defined(ENABLE_OPUS_CUSTOM_API)

int opus_custom_encoder_init(CELTEncoder *st, const CELTMode *mode, int channels)

{

   return opus_custom_encoder_init_arch(st, mode, channels, opus_select_arch());

}

#endif

int celt_encoder_init(CELTEncoder *st, opus_int32 sampling_rate, int channels,

                      int arch)

{

   int ret;

#ifdef ENABLE_QEXT

   if (sampling_rate==96000) {

      st->upsample = 1;

      return opus_custom_encoder_init_arch(st,

              opus_custom_mode_create(96000, 1920, NULL), channels, arch);

   }

#endif

   ret = opus_custom_encoder_init_arch(st,

           opus_custom_mode_create(48000, 960, NULL), channels, arch);

   if (ret != OPUS_OK)

      return ret;

   st->upsample = resampling_factor(sampling_rate);

   return OPUS_OK;

}

#if defined(CUSTOM_MODES) || defined(ENABLE_OPUS_CUSTOM_API)

void opus_custom_encoder_destroy(CELTEncoder *st)

{

   opus_free(st);

}

#endif /* CUSTOM_MODES */

static int transient_analysis(const opus_val32 * OPUS_RESTRICT in, int len, int C,

                              opus_val16 *tf_estimate, int *tf_chan, int allow_weak_transients,

                              int *weak_transient, opus_val16 tone_freq, opus_val32 toneishness)

{

   int i;

   VARDECL(opus_val16, tmp);

   opus_val32 mem0,mem1;

   int is_transient = 0;

   opus_int32 mask_metric = 0;

   int c;

   opus_val16 tf_max;

   int len2;

   /* Forward masking: 6.7 dB/ms. */

#ifdef FIXED_POINT

   int forward_shift = 4;

#else

   opus_val16 forward_decay = QCONST16(.0625f,15);

#endif

   /* Table of 6*64/x, trained on real data to minimize the average error */

   static const unsigned char inv_table[128] = {

         255,255,156,110, 86, 70, 59, 51, 45, 40, 37, 33, 31, 28, 26, 25,

          23, 22, 21, 20, 19, 18, 17, 16, 16, 15, 15, 14, 13, 13, 12, 12,

          12, 12, 11, 11, 11, 10, 10, 10,  9,  9,  9,  9,  9,  9,  8,  8,

           8,  8,  8,  7,  7,  7,  7,  7,  7,  6,  6,  6,  6,  6,  6,  6,

           6,  6,  6,  6,  6,  6,  6,  6,  6,  5,  5,  5,  5,  5,  5,  5,

           5,  5,  5,  5,  5,  4,  4,  4,  4,  4,  4,  4,  4,  4,  4,  4,

           4,  4,  4,  4,  4,  4,  4,  4,  4,  4,  4,  4,  4,  4,  3,  3,

           3,  3,  3,  3,  3,  3,  3,  3,  3,  3,  3,  3,  3,  3,  3,  2,

   };

   SAVE_STACK;

#ifdef FIXED_POINT

   int in_shift = IMAX(0, celt_ilog2(1+celt_maxabs32(in, C*len))-14);

#endif

   ALLOC(tmp, len, opus_val16);

   *weak_transient = 0;

   /* For lower bitrates, let's be more conservative and have a forward masking

      decay of 3.3 dB/ms. This avoids having to code transients at very low

      bitrate (mostly for hybrid), which can result in unstable energy and/or

      partial collapse. */

   if (allow_weak_transients)

   {

#ifdef FIXED_POINT

      forward_shift = 5;

#else

      forward_decay = QCONST16(.03125f,15);

#endif

   }

   len2=len/2;

   for (c=0;c<C;c++)

   {

      opus_val32 mean;

      opus_int32 unmask=0;

      opus_val32 norm;

      opus_val16 maxE;

      mem0=0;

      mem1=0;

      /* High-pass filter: (1 - 2*z^-1 + z^-2) / (1 - z^-1 + .5*z^-2) */

      for (i=0;i<len;i++)

      {

#ifndef FIXED_POINT

         float mem00;

#endif

         opus_val32 x,y;

         x = SHR32(in[i+c*len],in_shift);

         y = ADD32(mem0, x);

#ifdef FIXED_POINT

         mem0 = mem1 + y - SHL32(x,1);

         mem1 = x - SHR32(y,1);

#else

         /* Original code:

         mem0 = mem1 + y - 2*x;

         mem1 = x - .5f*y;

         Modified code to shorten dependency chains: */

         mem00=mem0;

         mem0 = mem0 - x + .5f*mem1;

         mem1 =  x - mem00;

#endif

         tmp[i] = SROUND16(y, 2);

         /*printf("%f ", tmp[i]);*/

      }

      /*printf("\n");*/

      /* First few samples are bad because we don't propagate the memory */

      OPUS_CLEAR(tmp, 12);

#ifdef FIXED_POINT

      /* Normalize tmp to max range */

      {

         int shift=0;

         shift = 14-celt_ilog2(MAX16(1, celt_maxabs16(tmp, len)));

         if (shift!=0)

         {

            for (i=0;i<len;i++)

               tmp[i] = SHL16(tmp[i], shift);

         }

      }

#endif

      mean=0;

      mem0=0;

      /* Grouping by two to reduce complexity */

      /* Forward pass to compute the post-echo threshold*/

      for (i=0;i<len2;i++)

      {

         opus_val32 x2 = PSHR32(MULT16_16(tmp[2*i],tmp[2*i]) + MULT16_16(tmp[2*i+1],tmp[2*i+1]),4);

         mean += PSHR32(x2, 12);

#ifdef FIXED_POINT

         /* FIXME: Use PSHR16() instead */

         mem0 = mem0 + PSHR32(x2-mem0,forward_shift);

         tmp[i] = PSHR32(mem0, 12);

#else

         mem0 = x2 + (1.f-forward_decay)*mem0;

         tmp[i] = forward_decay*mem0;

#endif

      }

      mem0=0;

      maxE=0;

      /* Backward pass to compute the pre-echo threshold */

      for (i=len2-1;i>=0;i--)

      {

         /* Backward masking: 13.9 dB/ms. */

#ifdef FIXED_POINT

         /* FIXME: Use PSHR16() instead */

         mem0 = mem0 + PSHR32(SHL32(tmp[i],4)-mem0,3);

         tmp[i] = PSHR32(mem0, 4);

         maxE = MAX16(maxE, tmp[i]);

#else

         mem0 = tmp[i] + 0.875f*mem0;

         tmp[i] = 0.125f*mem0;

         maxE = MAX16(maxE, 0.125f*mem0);

#endif

      }

      /*for (i=0;i<len2;i++)printf("%f ", tmp[i]/mean);printf("\n");*/

      /* Compute the ratio of the "frame energy" over the harmonic mean of the energy.

         This essentially corresponds to a bitrate-normalized temporal noise-to-mask

         ratio */

      /* As a compromise with the old transient detector, frame energy is the

         geometric mean of the energy and half the max */

#ifdef FIXED_POINT

      /* Costs two sqrt() to avoid overflows */

      mean = MULT16_16(celt_sqrt(mean), celt_sqrt(MULT16_16(maxE,len2>>1)));

#else

      mean = celt_sqrt(mean * maxE*.5*len2);

#endif

      /* Inverse of the mean energy in Q15+6 */

      norm = SHL32(EXTEND32(len2),6+14)/ADD32(EPSILON,SHR32(mean,1));

      /* Compute harmonic mean discarding the unreliable boundaries

         The data is smooth, so we only take 1/4th of the samples */

      unmask=0;

      /* We should never see NaNs here. If we find any, then something really bad happened and we better abort

         before it does any damage later on. If these asserts are disabled (no hardening), then the table

         lookup a few lines below (id = ...) is likely to crash dur to an out-of-bounds read. DO NOT FIX

         that crash on NaN since it could result in a worse issue later on. */

      celt_assert(!celt_isnan(tmp[0]));

      celt_assert(!celt_isnan(norm));

      for (i=12;i<len2-5;i+=4)

      {

         int id;

#ifdef FIXED_POINT

         id = MAX32(0,MIN32(127,MULT16_32_Q15(tmp[i]+EPSILON,norm))); /* Do not round to nearest */

#else

         id = (int)MAX32(0,MIN32(127,floor(64*norm*(tmp[i]+EPSILON)))); /* Do not round to nearest */

#endif

         unmask += inv_table[id];

      }

      /*printf("%d\n", unmask);*/

      /* Normalize, compensate for the 1/4th of the sample and the factor of 6 in the inverse table */

      unmask = 64*unmask*4/(6*(len2-17));

      if (unmask>mask_metric)

      {

         *tf_chan = c;

         mask_metric = unmask;

      }

   }

   is_transient = mask_metric>200;

   /* Prevent the transient detector from confusing the partial cycle of a

      very low frequency tone with a transient. */

   if (toneishness > QCONST32(.98f, 29) && tone_freq < QCONST16(0.026f, 13))

   {

      is_transient = 0;

      mask_metric = 0;

   }

   /* For low bitrates, define "weak transients" that need to be

      handled differently to avoid partial collapse. */

   if (allow_weak_transients && is_transient && mask_metric<600) {

      is_transient = 0;

      *weak_transient = 1;

   }

   /* Arbitrary metric for VBR boost */

   tf_max = MAX16(0,celt_sqrt(27*mask_metric)-42);

   /* *tf_estimate = 1 + MIN16(1, sqrt(MAX16(0, tf_max-30))/20); */

   *tf_estimate = celt_sqrt(MAX32(0, SHL32(MULT16_16(QCONST16(0.0069,14),MIN16(163,tf_max)),14)-QCONST32(0.139,28)));

   /*printf("%d %f\n", tf_max, mask_metric);*/

   RESTORE_STACK;

#ifdef FUZZING

   is_transient = rand()&0x1;

#endif

   /*printf("%d %f %d\n", is_transient, (float)*tf_estimate, tf_max);*/

   return is_transient;

}

/* Looks for sudden increases of energy to decide whether we need to patch

   the transient decision */

static int patch_transient_decision(celt_glog *newE, celt_glog *oldE, int nbEBands,

      int start, int end, int C)

{

   int i, c;

   opus_val32 mean_diff=0;

   celt_glog spread_old[26];

   /* Apply an aggressive (-6 dB/Bark) spreading function to the old frame to

      avoid false detection caused by irrelevant bands */

   if (C==1)

   {

      spread_old[start] = oldE[start];

      for (i=start+1;i<end;i++)

         spread_old[i] = MAXG(spread_old[i-1]-GCONST(1.0f), oldE[i]);

   } else {

      spread_old[start] = MAXG(oldE[start],oldE[start+nbEBands]);

      for (i=start+1;i<end;i++)

         spread_old[i] = MAXG(spread_old[i-1]-GCONST(1.0f),

                               MAXG(oldE[i],oldE[i+nbEBands]));

   }

   for (i=end-2;i>=start;i--)

      spread_old[i] = MAXG(spread_old[i], spread_old[i+1]-GCONST(1.0f));

   /* Compute mean increase */

   c=0; do {

      for (i=IMAX(2,start);i<end-1;i++)

      {

         opus_val16 x1, x2;

         x1 = MAXG(0, newE[i + c*nbEBands]);

         x2 = MAXG(0, spread_old[i]);

         mean_diff = ADD32(mean_diff, MAXG(0, SUB32(x1, x2)));

      }

   } while (++c<C);

   mean_diff = DIV32(mean_diff, C*(end-1-IMAX(2,start)));

   /*printf("%f %f %d\n", mean_diff, max_diff, count);*/

   return mean_diff > GCONST(1.f);

}

/** Apply window and compute the MDCT for all sub-frames and

    all channels in a frame */

static void compute_mdcts(const CELTMode *mode, int shortBlocks, celt_sig * OPUS_RESTRICT in,

                          celt_sig * OPUS_RESTRICT out, int C, int CC, int LM, int upsample,

                          int arch)

{

   const int overlap = mode->overlap;

   int N;

   int B;

   int shift;

   int i, b, c;

   if (shortBlocks)

   {

      B = shortBlocks;

      N = mode->shortMdctSize;

      shift = mode->maxLM;

   } else {

      B = 1;

      N = mode->shortMdctSize<<LM;

      shift = mode->maxLM-LM;

   }

   c=0; do {

      for (b=0;b<B;b++)

      {

         /* Interleaving the sub-frames while doing the MDCTs */

         clt_mdct_forward(&mode->mdct, in+c*(B*N+overlap)+b*N,

                          &out[b+c*N*B], mode->window, overlap, shift, B,

                          arch);

      }

   } while (++c<CC);

   if (CC==2&&C==1)

   {

      for (i=0;i<B*N;i++)

         out[i] = ADD32(HALF32(out[i]), HALF32(out[B*N+i]));

   }

   if (upsample != 1)

   {

      c=0; do

      {

         int bound = B*N/upsample;

         for (i=0;i<bound;i++)

            out[c*B*N+i] *= upsample;

         OPUS_CLEAR(&out[c*B*N+bound], B*N-bound);

      } while (++c<C);

   }

}

void celt_preemphasis(const opus_res * OPUS_RESTRICT pcmp, celt_sig * OPUS_RESTRICT inp,

                        int N, int CC, int upsample, const opus_val16 *coef, celt_sig *mem, int clip)

{

   int i;

   opus_val16 coef0;

   celt_sig m;

   int Nu;

   coef0 = coef[0];

   m = *mem;

   /* Fast path for the normal 48kHz case and no clipping */

   if (coef[1] == 0 && upsample == 1 && !clip)

   {

      for (i=0;i<N;i++)

      {

         celt_sig x;

         x = RES2SIG(pcmp[CC*i]);

         /* Apply pre-emphasis */

         inp[i] = x - m;

         m = MULT16_32_Q15(coef0, x);

      }

      *mem = m;

      return;

   }

   Nu = N/upsample;

   if (upsample!=1)

   {

      OPUS_CLEAR(inp, N);

   }

   for (i=0;i<Nu;i++)

      inp[i*upsample] = RES2SIG(pcmp[CC*i]);

#ifndef FIXED_POINT

   if (clip)

   {

      /* Clip input to avoid encoding non-portable files */

      for (i=0;i<Nu;i++)

         inp[i*upsample] = MAX32(-65536.f, MIN32(65536.f,inp[i*upsample]));

   }

#else

   (void)clip; /* Avoids a warning about clip being unused. */

#endif

#if defined(CUSTOM_MODES) || defined(ENABLE_OPUS_CUSTOM_API) || defined(ENABLE_QEXT)

   if (coef[1] != 0)

   {

      opus_val16 coef1 = coef[1];

#if defined(FIXED_POINT) && defined(ENABLE_QEXT)

      /* If we need the extra precision, we use the fact that coef[3] is exact to do a Newton-Raphson

         iteration and get us more precision on coef[2]. */

      opus_val32 coef2_q30 = SHL32(coef[2], 18) + PSHR32(MULT16_16(QCONST32(1.f, 25) - MULT16_16(coef[3], coef[2]), coef[2]), 7);

      celt_assert(SIG_SHIFT == 12);

#else

      opus_val16 coef2 = coef[2];

#endif

      for (i=0;i<N;i++)

      {

         celt_sig x, tmp;

         x = inp[i];

         /* Apply pre-emphasis */

#if defined(FIXED_POINT) && defined(ENABLE_QEXT)

         tmp = SHL32(MULT32_32_Q31(coef2_q30, x), 1);

#else

         tmp = SHL32(MULT16_32_Q15(coef2, x), 15-SIG_SHIFT);

#endif

         inp[i] = tmp + m;

         m = MULT16_32_Q15(coef1, inp[i]) - MULT16_32_Q15(coef0, tmp);

      }

   } else

#endif

   {

      for (i=0;i<N;i++)

      {

         celt_sig x;

         x = inp[i];

         /* Apply pre-emphasis */

         inp[i] = x - m;

         m = MULT16_32_Q15(coef0, x);

      }

   }

   *mem = m;

}

static opus_val32 l1_metric(const celt_norm *tmp, int N, int LM, opus_val16 bias)

{

   int i;

   opus_val32 L1;

   L1 = 0;

   for (i=0;i<N;i++)

      L1 += EXTEND32(ABS16(SHR32(tmp[i], NORM_SHIFT-14)));

   /* When in doubt, prefer good freq resolution */

   L1 = MAC16_32_Q15(L1, LM*bias, L1);

   return L1;

}

static int tf_analysis(const CELTMode *m, int len, int isTransient,

      int *tf_res, int lambda, celt_norm *X, int N0, int LM,

      opus_val16 tf_estimate, int tf_chan, int *importance)

{

   int i;

   VARDECL(int, metric);

   int cost0;

   int cost1;

   VARDECL(int, path0);

   VARDECL(int, path1);

   VARDECL(celt_norm, tmp);

   VARDECL(celt_norm, tmp_1);

   int sel;

   int selcost[2];

   int tf_select=0;

   opus_val16 bias;

   SAVE_STACK;

   bias = MULT16_16_Q14(QCONST16(.04f,15), MAX16(-QCONST16(.25f,14), QCONST16(.5f,14)-tf_estimate));

   /*printf("%f ", bias);*/

   ALLOC(metric, len, int);

   ALLOC(tmp, (m->eBands[len]-m->eBands[len-1])<<LM, celt_norm);

   ALLOC(tmp_1, (m->eBands[len]-m->eBands[len-1])<<LM, celt_norm);

   ALLOC(path0, len, int);

   ALLOC(path1, len, int);

   for (i=0;i<len;i++)

   {

      int k, N;

      int narrow;

      opus_val32 L1, best_L1;

      int best_level=0;

      N = (m->eBands[i+1]-m->eBands[i])<<LM;

      /* band is too narrow to be split down to LM=-1 */

      narrow = (m->eBands[i+1]-m->eBands[i])==1;

      OPUS_COPY(tmp, &X[tf_chan*N0 + (m->eBands[i]<<LM)], N);

      /* Just add the right channel if we're in stereo */

      /*if (C==2)

         for (j=0;j<N;j++)

            tmp[j] = ADD16(SHR16(tmp[j], 1),SHR16(X[N0+j+(m->eBands[i]<<LM)], 1));*/

      L1 = l1_metric(tmp, N, isTransient ? LM : 0, bias);

      best_L1 = L1;

      /* Check the -1 case for transients */

      if (isTransient && !narrow)

      {

         OPUS_COPY(tmp_1, tmp, N);

         haar1(tmp_1, N>>LM, 1<<LM);

         L1 = l1_metric(tmp_1, N, LM+1, bias);

         if (L1<best_L1)

         {

            best_L1 = L1;

            best_level = -1;

         }

      }

      /*printf ("%f ", L1);*/

      for (k=0;k<LM+!(isTransient||narrow);k++)

      {

         int B;

         if (isTransient)

            B = (LM-k-1);

         else

            B = k+1;

         haar1(tmp, N>>k, 1<<k);

         L1 = l1_metric(tmp, N, B, bias);

         if (L1 < best_L1)

         {

            best_L1 = L1;

            best_level = k+1;

         }

      }

      /*printf ("%d ", isTransient ? LM-best_level : best_level);*/

      /* metric is in Q1 to be able to select the mid-point (-0.5) for narrower bands */

      if (isTransient)

         metric[i] = 2*best_level;

      else

         metric[i] = -2*best_level;

      /* For bands that can't be split to -1, set the metric to the half-way point to avoid

         biasing the decision */

      if (narrow && (metric[i]==0 || metric[i]==-2*LM))

         metric[i]-=1;

      /*printf("%d ", metric[i]/2 + (!isTransient)*LM);*/

   }

   /*printf("\n");*/

   /* Search for the optimal tf resolution, including tf_select */

   tf_select = 0;

   for (sel=0;sel<2;sel++)

   {

      cost0 = importance[0]*abs(metric[0]-2*tf_select_table[LM][4*isTransient+2*sel+0]);

      cost1 = importance[0]*abs(metric[0]-2*tf_select_table[LM][4*isTransient+2*sel+1]) + (isTransient ? 0 : lambda);

      for (i=1;i<len;i++)

      {

         int curr0, curr1;

         curr0 = IMIN(cost0, cost1 + lambda);

         curr1 = IMIN(cost0 + lambda, cost1);

         cost0 = curr0 + importance[i]*abs(metric[i]-2*tf_select_table[LM][4*isTransient+2*sel+0]);

         cost1 = curr1 + importance[i]*abs(metric[i]-2*tf_select_table[LM][4*isTransient+2*sel+1]);

      }

      cost0 = IMIN(cost0, cost1);

      selcost[sel]=cost0;

   }

   /* For now, we're conservative and only allow tf_select=1 for transients.

    * If tests confirm it's useful for non-transients, we could allow it. */

   if (selcost[1]<selcost[0] && isTransient)

      tf_select=1;

   cost0 = importance[0]*abs(metric[0]-2*tf_select_table[LM][4*isTransient+2*tf_select+0]);

   cost1 = importance[0]*abs(metric[0]-2*tf_select_table[LM][4*isTransient+2*tf_select+1]) + (isTransient ? 0 : lambda);

   /* Viterbi forward pass */

   for (i=1;i<len;i++)

   {

      int curr0, curr1;

      int from0, from1;

      from0 = cost0;

      from1 = cost1 + lambda;

      if (from0 < from1)

      {

         curr0 = from0;

         path0[i]= 0;

      } else {

         curr0 = from1;

         path0[i]= 1;

      }

      from0 = cost0 + lambda;

      from1 = cost1;

      if (from0 < from1)

      {

         curr1 = from0;

         path1[i]= 0;

      } else {

         curr1 = from1;

         path1[i]= 1;

      }

      cost0 = curr0 + importance[i]*abs(metric[i]-2*tf_select_table[LM][4*isTransient+2*tf_select+0]);

      cost1 = curr1 + importance[i]*abs(metric[i]-2*tf_select_table[LM][4*isTransient+2*tf_select+1]);

   }

   tf_res[len-1] = cost0 < cost1 ? 0 : 1;

   /* Viterbi backward pass to check the decisions */

   for (i=len-2;i>=0;i--)

   {

      if (tf_res[i+1] == 1)

         tf_res[i] = path1[i+1];

      else

         tf_res[i] = path0[i+1];

   }

   /*printf("%d %f\n", *tf_sum, tf_estimate);*/

   RESTORE_STACK;

#ifdef FUZZING

   tf_select = rand()&0x1;

   tf_res[0] = rand()&0x1;

   for (i=1;i<len;i++)

      tf_res[i] = tf_res[i-1] ^ ((rand()&0xF) == 0);

#endif

   return tf_select;

}

static void tf_encode(int start, int end, int isTransient, int *tf_res, int LM, int tf_select, ec_enc *enc)

{

   int curr, i;

   int tf_select_rsv;

   int tf_changed;

   int logp;

   opus_uint32 budget;

   opus_uint32 tell;

   budget = enc->storage*8;

   tell = ec_tell(enc);

   logp = isTransient ? 2 : 4;

   /* Reserve space to code the tf_select decision. */

   tf_select_rsv = LM>0 && tell+logp+1 <= budget;

   budget -= tf_select_rsv;

   curr = tf_changed = 0;

   for (i=start;i<end;i++)

   {

      if (tell+logp<=budget)

      {

         ec_enc_bit_logp(enc, tf_res[i] ^ curr, logp);

         tell = ec_tell(enc);

         curr = tf_res[i];

         tf_changed |= curr;

      }

      else

         tf_res[i] = curr;

      logp = isTransient ? 4 : 5;

   }

   /* Only code tf_select if it would actually make a difference. */

   if (tf_select_rsv &&

         tf_select_table[LM][4*isTransient+0+tf_changed]!=

         tf_select_table[LM][4*isTransient+2+tf_changed])

      ec_enc_bit_logp(enc, tf_select, 1);

   else

      tf_select = 0;

   for (i=start;i<end;i++)

      tf_res[i] = tf_select_table[LM][4*isTransient+2*tf_select+tf_res[i]];

   /*for(i=0;i<end;i++)printf("%d ", isTransient ? tf_res[i] : LM+tf_res[i]);printf("\n");*/

}

static int alloc_trim_analysis(const CELTMode *m, const celt_norm *X,

      const celt_glog *bandLogE, int end, int LM, int C, int N0,

      AnalysisInfo *analysis, opus_val16 *stereo_saving, opus_val16 tf_estimate,

      int intensity, celt_glog surround_trim, opus_int32 equiv_rate, int arch)

{

   int i;

   opus_val32 diff=0;

   int c;

   int trim_index;

   opus_val16 trim = QCONST16(5.f, 8);

   opus_val16 logXC, logXC2;

   /* At low bitrate, reducing the trim seems to help. At higher bitrates, it's less

      clear what's best, so we're keeping it as it was before, at least for now. */

   if (equiv_rate < 64000) {

      trim = QCONST16(4.f, 8);

   } else if (equiv_rate < 80000) {

      opus_int32 frac = (equiv_rate-64000) >> 10;

      trim = QCONST16(4.f, 8) + QCONST16(1.f/16.f, 8)*frac;

   }

   if (C==2)

   {

      opus_val16 sum = 0; /* Q10 */

      opus_val16 minXC; /* Q10 */

      /* Compute inter-channel correlation for low frequencies */

      for (i=0;i<8;i++)

      {

         opus_val32 partial;

         partial = celt_inner_prod_norm_shift(&X[m->eBands[i]<<LM], &X[N0+(m->eBands[i]<<LM)],

               (m->eBands[i+1]-m->eBands[i])<<LM, arch);

         sum = ADD16(sum, EXTRACT16(SHR32(partial, 18)));

      }

      sum = MULT16_16_Q15(QCONST16(1.f/8, 15), sum);

      sum = MIN16(QCONST16(1.f, 10), ABS16(sum));

      minXC = sum;

      for (i=8;i<intensity;i++)

      {

         opus_val32 partial;

         partial = celt_inner_prod_norm_shift(&X[m->eBands[i]<<LM], &X[N0+(m->eBands[i]<<LM)],

               (m->eBands[i+1]-m->eBands[i])<<LM, arch);

         minXC = MIN16(minXC, ABS16(EXTRACT16(SHR32(partial, 18))));

      }

      minXC = MIN16(QCONST16(1.f, 10), ABS16(minXC));

      /*printf ("%f\n", sum);*/

      /* mid-side savings estimations based on the LF average*/

      logXC = celt_log2(QCONST32(1.001f, 20)-MULT16_16(sum, sum));

      /* mid-side savings estimations based on min correlation */

      logXC2 = MAX16(HALF16(logXC), celt_log2(QCONST32(1.001f, 20)-MULT16_16(minXC, minXC)));

#ifdef FIXED_POINT

      /* Compensate for Q20 vs Q14 input and convert output to Q8 */

      logXC = PSHR32(logXC-QCONST16(6.f, 10),10-8);

      logXC2 = PSHR32(logXC2-QCONST16(6.f, 10),10-8);

#endif

      trim += MAX16(-QCONST16(4.f, 8), MULT16_16_Q15(QCONST16(.75f,15),logXC));

      *stereo_saving = MIN16(*stereo_saving + QCONST16(0.25f, 8), -HALF16(logXC2));

   }

   /* Estimate spectral tilt */

   c=0; do {

      for (i=0;i<end-1;i++)

      {

         diff += SHR32(bandLogE[i+c*m->nbEBands], 5)*(opus_int32)(2+2*i-end);