celt.c

/* 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 FOUNDATION 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_C

#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 "plc.h"
#include "vq.h"

static const unsigned char trim_icdf[11] = {126, 124, 119, 109, 87, 41, 19, 9, 4, 2, 0};
/* Probs: NONE: 21.875%, LIGHT: 6.25%, NORMAL: 65.625%, AGGRESSIVE: 6.25% */
static const unsigned char spread_icdf[4] = {25, 23, 2, 0};

static const unsigned char tapset_icdf[3]={2,1,0};

static const unsigned char toOpusTable[20] = {
      0xE0, 0xE8, 0xF0, 0xF8,
      0xC0, 0xC8, 0xD0, 0xD8,
      0xA0, 0xA8, 0xB0, 0xB8,
      0x00, 0x00, 0x00, 0x00,
      0x80, 0x88, 0x90, 0x98,
};

static const unsigned char fromOpusTable[16] = {
      0x80, 0x88, 0x90, 0x98,
      0x40, 0x48, 0x50, 0x58,
      0x20, 0x28, 0x30, 0x38,
      0x00, 0x08, 0x10, 0x18
};

static inline int toOpus(unsigned char c)
{
   int ret=0;
   if (c<0xA0)
      ret = toOpusTable[c>>3];
   if (ret == 0)
      return -1;
   else
      return ret|(c&0x7);
}

static inline int fromOpus(unsigned char c)
{
   if (c<0x80)
      return -1;
   else
      return fromOpusTable[(c>>3)-16] | (c&0x7);
}

#define COMBFILTER_MAXPERIOD 1024
#define COMBFILTER_MINPERIOD 15

static int resampling_factor(opus_int32 rate)
{
   int ret;
   switch (rate)
   {
   case 48000:
      ret = 1;
      break;
   case 24000:
      ret = 2;
      break;
   case 16000:
      ret = 3;
      break;
   case 12000:
      ret = 4;
      break;
   case 8000:
      ret = 6;
      break;
   default:
      ret = 0;
   }
   return ret;
}

/** Encoder state
 @brief Encoder state
 */
struct CELTEncoder {
   const CELTMode *mode;     /**< Mode used by the encoder */
   int overlap;
   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;

   /* 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;

   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;

#ifdef RESYNTH
   celt_sig syn_mem[2][2*MAX_PERIOD];
#endif

   celt_sig in_mem[1]; /* Size = channels*mode->overlap */
   /* celt_sig prefilter_mem[],  Size = channels*COMBFILTER_PERIOD */
   /* celt_sig overlap_mem[],  Size = channels*mode->overlap */
   /* opus_val16 oldEBands[], Size = 2*channels*mode->nbEBands */
};

int celt_encoder_get_size(int channels)
{
   CELTMode *mode = celt_mode_create(48000, 960, NULL);
   return celt_encoder_get_size_custom(mode, channels);
}

int celt_encoder_get_size_custom(const CELTMode *mode, int channels)
{
   int size = sizeof(struct CELTEncoder)
         + (2*channels*mode->overlap-1)*sizeof(celt_sig)
         + channels*COMBFILTER_MAXPERIOD*sizeof(celt_sig)
         + 3*channels*mode->nbEBands*sizeof(opus_val16);
   return size;
}

CELTEncoder *celt_encoder_create(int sampling_rate, int channels, int *error)
{
   CELTEncoder *st;
   st = (CELTEncoder *)celt_alloc(celt_encoder_get_size(channels));
   if (st!=NULL && celt_encoder_init(st, sampling_rate, channels, error)==NULL)
   {
      celt_encoder_destroy(st);
      st = NULL;
   }
   return st;
}

CELTEncoder *celt_encoder_create_custom(const CELTMode *mode, int channels, int *error)
{
   CELTEncoder *st = (CELTEncoder *)celt_alloc(celt_encoder_get_size_custom(mode, channels));
   if (st!=NULL && celt_encoder_init_custom(st, mode, channels, error)==NULL)
   {
      celt_encoder_destroy(st);
      st = NULL;
   }
   return st;
}

CELTEncoder *celt_encoder_init(CELTEncoder *st, int sampling_rate, int channels, int *error)
{
   celt_encoder_init_custom(st, celt_mode_create(48000, 960, NULL), channels, error);
   st->upsample = resampling_factor(sampling_rate);
   if (st->upsample==0)
   {
      if (error)
         *error = CELT_BAD_ARG;
      return NULL;
   }
   return st;
}

CELTEncoder *celt_encoder_init_custom(CELTEncoder *st, const CELTMode *mode, int channels, int *error)
{
   if (channels < 0 || channels > 2)
   {
      if (error)
         *error = CELT_BAD_ARG;
      return NULL;
   }

   if (st==NULL || mode==NULL)
   {
      if (error)
         *error = CELT_ALLOC_FAIL;
      return NULL;
   }

   CELT_MEMSET((char*)st, 0, celt_encoder_get_size_custom(mode, channels));

   st->mode = mode;
   st->overlap = mode->overlap;
   st->stream_channels = st->channels = channels;

   st->upsample = 1;
   st->start = 0;
   st->end = st->mode->effEBands;
   st->signalling = 1;

   st->constrained_vbr = 1;
   st->clip = 1;

   st->bitrate = 255000*channels;
   st->vbr = 0;
   st->vbr_offset = 0;
   st->force_intra  = 0;
   st->delayedIntra = 1;
   st->tonal_average = 256;
   st->spread_decision = SPREAD_NORMAL;
   st->hf_average = 0;
   st->tapset_decision = 0;
   st->complexity = 5;

   if (error)
      *error = CELT_OK;
   return st;
}

void celt_encoder_destroy(CELTEncoder *st)
{
   celt_free(st);
}

static inline opus_int16 FLOAT2INT16(float x)
{
   x = x*CELT_SIG_SCALE;
   x = MAX32(x, -32768);
   x = MIN32(x, 32767);
   return (opus_int16)float2int(x);
}

static inline opus_val16 SIG2WORD16(celt_sig x)
{
#ifdef FIXED_POINT
   x = PSHR32(x, SIG_SHIFT);
   x = MAX32(x, -32768);
   x = MIN32(x, 32767);
   return EXTRACT16(x);
#else
   return (opus_val16)x;
#endif
}

static int transient_analysis(const opus_val32 * restrict in, int len, int C,
                              int overlap)
{
   int i;
   VARDECL(opus_val16, tmp);
   opus_val32 mem0=0,mem1=0;
   int is_transient = 0;
   int block;
   int N;
   VARDECL(opus_val16, bins);
   SAVE_STACK;
   ALLOC(tmp, len, opus_val16);

   block = overlap/2;
   N=len/block;
   ALLOC(bins, N, opus_val16);
   if (C==1)
   {
      for (i=0;i<len;i++)
         tmp[i] = SHR32(in[i],SIG_SHIFT);
   } else {
      for (i=0;i<len;i++)
         tmp[i] = SHR32(ADD32(in[i],in[i+len]), SIG_SHIFT+1);
   }

   /* High-pass filter: (1 - 2*z^-1 + z^-2) / (1 - z^-1 + .5*z^-2) */
   for (i=0;i<len;i++)
   {
      opus_val32 x,y;
      x = tmp[i];
      y = ADD32(mem0, x);
#ifdef FIXED_POINT
      mem0 = mem1 + y - SHL32(x,1);
      mem1 = x - SHR32(y,1);
#else
      mem0 = mem1 + y - 2*x;
      mem1 = x - .5f*y;
#endif
      tmp[i] = EXTRACT16(SHR(y,2));
   }
   /* First few samples are bad because we don't propagate the memory */
   for (i=0;i<12;i++)
      tmp[i] = 0;

   for (i=0;i<N;i++)
   {
      int j;
      opus_val16 max_abs=0;
      for (j=0;j<block;j++)
         max_abs = MAX16(max_abs, ABS16(tmp[i*block+j]));
      bins[i] = max_abs;
   }
   for (i=0;i<N;i++)
   {
      int j;
      int conseq=0;
      opus_val16 t1, t2, t3;

      t1 = MULT16_16_Q15(QCONST16(.15f, 15), bins[i]);
      t2 = MULT16_16_Q15(QCONST16(.4f, 15), bins[i]);
      t3 = MULT16_16_Q15(QCONST16(.15f, 15), bins[i]);
      for (j=0;j<i;j++)
      {
         if (bins[j] < t1)
            conseq++;
         if (bins[j] < t2)
            conseq++;
         else
            conseq = 0;
      }
      if (conseq>=3)
         is_transient=1;
      conseq = 0;
      for (j=i+1;j<N;j++)
      {
         if (bins[j] < t3)
            conseq++;
         else
            conseq = 0;
      }
      if (conseq>=7)
         is_transient=1;
   }
   RESTORE_STACK;
#ifdef FUZZING
   is_transient = rand()&0x1;
#endif
   return is_transient;
}

/** 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 * restrict in, celt_sig * restrict out, int _C, int LM)
{
   const int C = CHANNELS(_C);
   if (C==1 && !shortBlocks)
   {
      const int overlap = OVERLAP(mode);
      clt_mdct_forward(&mode->mdct, in, out, mode->window, overlap, mode->maxLM-LM, 1);
   } else {
      const int overlap = OVERLAP(mode);
      int N = mode->shortMdctSize<<LM;
      int B = 1;
      int b, c;
      if (shortBlocks)
      {
         N = mode->shortMdctSize;
         B = shortBlocks;
      }
      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, shortBlocks ? mode->maxLM : mode->maxLM-LM, B);
         }
      } while (++c<C);
   }
}

/** Compute the IMDCT and apply window for all sub-frames and
    all channels in a frame */
static void compute_inv_mdcts(const CELTMode *mode, int shortBlocks, celt_sig *X,
      celt_sig * restrict out_mem[],
      celt_sig * restrict overlap_mem[], int _C, int LM)
{
   int c;
   const int C = CHANNELS(_C);
   const int N = mode->shortMdctSize<<LM;
   const int overlap = OVERLAP(mode);
   VARDECL(opus_val32, x);
   SAVE_STACK;

   ALLOC(x, N+overlap, opus_val32);
   c=0; do {
      int j;
      int b;
      int N2 = N;
      int B = 1;

      if (shortBlocks)
      {
         N2 = mode->shortMdctSize;
         B = shortBlocks;
      }
      /* Prevents problems from the imdct doing the overlap-add */
      CELT_MEMSET(x, 0, overlap);

      for (b=0;b<B;b++)
      {
         /* IMDCT on the interleaved the sub-frames */
         clt_mdct_backward(&mode->mdct, &X[b+c*N2*B], x+N2*b, mode->window, overlap, shortBlocks ? mode->maxLM : mode->maxLM-LM, B);
      }

      for (j=0;j<overlap;j++)
         out_mem[c][j] = x[j] + overlap_mem[c][j];
      for (;j<N;j++)
         out_mem[c][j] = x[j];
      for (j=0;j<overlap;j++)
         overlap_mem[c][j] = x[N+j];
   } while (++c<C);
   RESTORE_STACK;
}

static void deemphasis(celt_sig *in[], opus_val16 *pcm, int N, int _C, int downsample, const opus_val16 *coef, celt_sig *mem)
{
   const int C = CHANNELS(_C);
   int c;
   int count=0;
   c=0; do {
      int j;
      celt_sig * restrict x;
      opus_val16  * restrict y;
      celt_sig m = mem[c];
      x =in[c];
      y = pcm+c;
      for (j=0;j<N;j++)
      {
         celt_sig tmp = *x + m;
         m = MULT16_32_Q15(coef[0], tmp)
           - MULT16_32_Q15(coef[1], *x);
         tmp = SHL32(MULT16_32_Q15(coef[3], tmp), 2);
         x++;
         /* Technically the store could be moved outside of the if because
            the stores we don't want will just be overwritten */
         if (++count==downsample)
         {
            *y = SCALEOUT(SIG2WORD16(tmp));
            y+=C;
            count=0;
         }
      }
      mem[c] = m;
   } while (++c<C);
}

static void comb_filter(opus_val32 *y, opus_val32 *x, int T0, int T1, int N,
      opus_val16 g0, opus_val16 g1, int tapset0, int tapset1,
      const opus_val16 *window, int overlap)
{
   int i;
   /* printf ("%d %d %f %f\n", T0, T1, g0, g1); */
   opus_val16 g00, g01, g02, g10, g11, g12;
   static const opus_val16 gains[3][3] = {
         {QCONST16(0.3066406250f, 15), QCONST16(0.2170410156f, 15), QCONST16(0.1296386719f, 15)},
         {QCONST16(0.4638671875f, 15), QCONST16(0.2680664062f, 15), QCONST16(0.f, 15)},
         {QCONST16(0.7998046875f, 15), QCONST16(0.1000976562f, 15), QCONST16(0.f, 15)}};
   g00 = MULT16_16_Q15(g0, gains[tapset0][0]);
   g01 = MULT16_16_Q15(g0, gains[tapset0][1]);
   g02 = MULT16_16_Q15(g0, gains[tapset0][2]);
   g10 = MULT16_16_Q15(g1, gains[tapset1][0]);
   g11 = MULT16_16_Q15(g1, gains[tapset1][1]);
   g12 = MULT16_16_Q15(g1, gains[tapset1][2]);
   for (i=0;i<overlap;i++)
   {
      opus_val16 f;
      f = MULT16_16_Q15(window[i],window[i]);
      y[i] = x[i]
               + MULT16_32_Q15(MULT16_16_Q15((Q15ONE-f),g00),x[i-T0])
               + MULT16_32_Q15(MULT16_16_Q15((Q15ONE-f),g01),x[i-T0-1])
               + MULT16_32_Q15(MULT16_16_Q15((Q15ONE-f),g01),x[i-T0+1])
               + MULT16_32_Q15(MULT16_16_Q15((Q15ONE-f),g02),x[i-T0-2])
               + MULT16_32_Q15(MULT16_16_Q15((Q15ONE-f),g02),x[i-T0+2])
               + MULT16_32_Q15(MULT16_16_Q15(f,g10),x[i-T1])
               + MULT16_32_Q15(MULT16_16_Q15(f,g11),x[i-T1-1])
               + MULT16_32_Q15(MULT16_16_Q15(f,g11),x[i-T1+1])
               + MULT16_32_Q15(MULT16_16_Q15(f,g12),x[i-T1-2])
               + MULT16_32_Q15(MULT16_16_Q15(f,g12),x[i-T1+2]);

   }
   for (i=overlap;i<N;i++)
      y[i] = x[i]
               + MULT16_32_Q15(g10,x[i-T1])
               + MULT16_32_Q15(g11,x[i-T1-1])
               + MULT16_32_Q15(g11,x[i-T1+1])
               + MULT16_32_Q15(g12,x[i-T1-2])
               + MULT16_32_Q15(g12,x[i-T1+2]);
}

static const signed char tf_select_table[4][8] = {
      {0, -1, 0, -1,    0,-1, 0,-1},
      {0, -1, 0, -2,    1, 0, 1,-1},
      {0, -2, 0, -3,    2, 0, 1,-1},
      {0, -2, 0, -3,    3, 0, 1,-1},
};

static opus_val32 l1_metric(const celt_norm *tmp, int N, int LM, int width)
{
   int i, j;
   static const opus_val16 sqrtM_1[4] = {Q15ONE, QCONST16(.70710678f,15), QCONST16(0.5f,15), QCONST16(0.35355339f,15)};
   opus_val32 L1;
   opus_val16 bias;
   L1=0;
   for (i=0;i<1<<LM;i++)
   {
      opus_val32 L2 = 0;
      for (j=0;j<N>>LM;j++)
         L2 = MAC16_16(L2, tmp[(j<<LM)+i], tmp[(j<<LM)+i]);
      L1 += celt_sqrt(L2);
   }
   L1 = MULT16_32_Q15(sqrtM_1[LM], L1);
   if (width==1)
      bias = QCONST16(.12f,15)*LM;
   else if (width==2)
      bias = QCONST16(.05f,15)*LM;
   else
      bias = QCONST16(.02f,15)*LM;
   L1 = MAC16_32_Q15(L1, bias, L1);
   return L1;
}

static int tf_analysis(const CELTMode *m, int len, int C, int isTransient,
      int *tf_res, int nbCompressedBytes, celt_norm *X, int N0, int LM,
      int *tf_sum)
{
   int i;
   VARDECL(int, metric);
   int cost0;
   int cost1;
   VARDECL(int, path0);
   VARDECL(int, path1);
   VARDECL(celt_norm, tmp);
   int lambda;
   int tf_select=0;
   SAVE_STACK;

   if (nbCompressedBytes<15*C)
   {
      *tf_sum = 0;
      for (i=0;i<len;i++)
         tf_res[i] = isTransient;
      return 0;
   }
   if (nbCompressedBytes<40)
      lambda = 12;
   else if (nbCompressedBytes<60)
      lambda = 6;
   else if (nbCompressedBytes<100)
      lambda = 4;
   else
      lambda = 3;

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

   *tf_sum = 0;
   for (i=0;i<len;i++)
   {
      int j, k, N;
      opus_val32 L1, best_L1;
      int best_level=0;
      N = (m->eBands[i+1]-m->eBands[i])<<LM;
      for (j=0;j<N;j++)
         tmp[j] = X[j+(m->eBands[i]<<LM)];
      /* Just add the right channel if we're in stereo */
      if (C==2)
         for (j=0;j<N;j++)
            tmp[j] = ADD16(tmp[j],X[N0+j+(m->eBands[i]<<LM)]);
      L1 = l1_metric(tmp, N, isTransient ? LM : 0, N>>LM);
      best_L1 = L1;
      /*printf ("%f ", L1);*/
      for (k=0;k<LM;k++)
      {
         int B;

         if (isTransient)
            B = (LM-k-1);
         else
            B = k+1;

         if (isTransient)
            haar1(tmp, N>>(LM-k), 1<<(LM-k));
         else
            haar1(tmp, N>>k, 1<<k);

         L1 = l1_metric(tmp, N, B, N>>LM);

         if (L1 < best_L1)
         {
            best_L1 = L1;
            best_level = k+1;
         }
      }
      /*printf ("%d ", isTransient ? LM-best_level : best_level);*/
      if (isTransient)
         metric[i] = best_level;
      else
         metric[i] = -best_level;
      *tf_sum += metric[i];
   }
   /*printf("\n");*/
   /* TODO: Detect the extreme transients that require tf_select = 1 */
   tf_select = 0;

   cost0 = 0;
   cost1 = 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 + abs(metric[i]-tf_select_table[LM][4*isTransient+2*tf_select+0]);
      cost1 = curr1 + abs(metric[i]-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];
   }
   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]];
   /*printf("%d %d ", isTransient, tf_select); for(i=0;i<end;i++)printf("%d ", tf_res[i]);printf("\n");*/
}

static void tf_decode(int start, int end, int isTransient, int *tf_res, int LM, ec_dec *dec)
{
   int i, curr, tf_select;
   int tf_select_rsv;
   int tf_changed;
   int logp;
   opus_uint32 budget;
   opus_uint32 tell;

   budget = dec->storage*8;
   tell = ec_tell(dec);
   logp = isTransient ? 2 : 4;
   tf_select_rsv = LM>0 && tell+logp+1<=budget;
   budget -= tf_select_rsv;
   tf_changed = curr = 0;
   for (i=start;i<end;i++)
   {
      if (tell+logp<=budget)
      {
         curr ^= ec_dec_bit_logp(dec, logp);
         tell = ec_tell(dec);
         tf_changed |= curr;
      }
      tf_res[i] = curr;
      logp = isTransient ? 4 : 5;
   }
   tf_select = 0;
   if (tf_select_rsv &&
     tf_select_table[LM][4*isTransient+0+tf_changed] !=
     tf_select_table[LM][4*isTransient+2+tf_changed])
   {
      tf_select = ec_dec_bit_logp(dec, 1);
   }
   for (i=start;i<end;i++)
   {
      tf_res[i] = tf_select_table[LM][4*isTransient+2*tf_select+tf_res[i]];
   }
}

static void init_caps(const CELTMode *m,int *cap,int LM,int C)
{
   int i;
   for (i=0;i<m->nbEBands;i++)
   {
      int N;
      N=(m->eBands[i+1]-m->eBands[i])<<LM;
      cap[i] = (m->cache.caps[m->nbEBands*(2*LM+C-1)+i]+64)*C*N>>2;
   }
}

static int alloc_trim_analysis(const CELTMode *m, const celt_norm *X,
      const opus_val16 *bandLogE, int end, int LM, int C, int N0)
{
   int i;
   opus_val32 diff=0;
   int c;
   int trim_index = 5;
   if (C==2)
   {
      opus_val16 sum = 0; /* Q10 */
      /* Compute inter-channel correlation for low frequencies */
      for (i=0;i<8;i++)
      {
         int j;
         opus_val32 partial = 0;
         for (j=m->eBands[i]<<LM;j<m->eBands[i+1]<<LM;j++)
            partial = MAC16_16(partial, X[j], X[N0+j]);
         sum = ADD16(sum, EXTRACT16(SHR32(partial, 18)));
      }
      sum = MULT16_16_Q15(QCONST16(1.f/8, 15), sum);
      /*printf ("%f\n", sum);*/
      if (sum > QCONST16(.995f,10))
         trim_index-=4;
      else if (sum > QCONST16(.92f,10))
         trim_index-=3;
      else if (sum > QCONST16(.85f,10))
         trim_index-=2;
      else if (sum > QCONST16(.8f,10))
         trim_index-=1;
   }

   /* Estimate spectral tilt */
   c=0; do {
      for (i=0;i<end-1;i++)
      {
         diff += bandLogE[i+c*m->nbEBands]*(opus_int32)(2+2*i-m->nbEBands);
      }
   } while (++c<C);
   /* We divide by two here to avoid making the tilt larger for stereo as a
      result of a bug in the loop above */
   diff /= 2*C*(end-1);
   /*printf("%f\n", diff);*/
   if (diff > QCONST16(2.f, DB_SHIFT))
      trim_index--;
   if (diff > QCONST16(8.f, DB_SHIFT))
      trim_index--;
   if (diff < -QCONST16(4.f, DB_SHIFT))
      trim_index++;
   if (diff < -QCONST16(10.f, DB_SHIFT))
      trim_index++;

   if (trim_index<0)
      trim_index = 0;
   if (trim_index>10)
      trim_index = 10;
#ifdef FUZZING
   trim_index = rand()%11;
#endif
   return trim_index;
}

static int stereo_analysis(const CELTMode *m, const celt_norm *X,
      int LM, int N0)
{
   int i;
   int thetas;
   opus_val32 sumLR = EPSILON, sumMS = EPSILON;

   /* Use the L1 norm to model the entropy of the L/R signal vs the M/S signal */
   for (i=0;i<13;i++)
   {
      int j;
      for (j=m->eBands[i]<<LM;j<m->eBands[i+1]<<LM;j++)
      {
         opus_val16 L, R, M, S;
         L = X[j];
         R = X[N0+j];
         M = L+R;
         S = L-R;
         sumLR += EXTEND32(ABS16(L)) + EXTEND32(ABS16(R));
         sumMS += EXTEND32(ABS16(M)) + EXTEND32(ABS16(S));
      }
   }
   sumMS = MULT16_32_Q15(QCONST16(0.707107f, 15), sumMS);
   thetas = 13;
   /* We don't need thetas for lower bands with LM<=1 */
   if (LM<=1)
      thetas -= 8;
   return MULT16_32_Q15((m->eBands[13]<<(LM+1))+thetas, sumMS)
         > MULT16_32_Q15(m->eBands[13]<<(LM+1), sumLR);
}

#ifdef FIXED_POINT
CELT_STATIC
int celt_encode_with_ec(CELTEncoder * restrict st, const opus_int16 * pcm, int frame_size, unsigned char *compressed, int nbCompressedBytes, ec_enc *enc)
{
#else
CELT_STATIC
int celt_encode_with_ec_float(CELTEncoder * restrict st, const celt_sig * pcm, int frame_size, unsigned char *compressed, int nbCompressedBytes, ec_enc *enc)
{
#endif
   int i, c, N;
   opus_int32 bits;
   ec_enc _enc;
   VARDECL(celt_sig, in);
   VARDECL(celt_sig, freq);
   VARDECL(celt_norm, X);
   VARDECL(celt_ener, bandE);
   VARDECL(opus_val16, bandLogE);
   VARDECL(int, fine_quant);
   VARDECL(opus_val16, error);
   VARDECL(int, pulses);
   VARDECL(int, cap);
   VARDECL(int, offsets);
   VARDECL(int, fine_priority);
   VARDECL(int, tf_res);
   VARDECL(unsigned char, collapse_masks);
   celt_sig *prefilter_mem;
   opus_val16 *oldBandE, *oldLogE, *oldLogE2;
   int shortBlocks=0;
   int isTransient=0;
   int resynth;
   const int CC = CHANNELS(st->channels);
   const int C = CHANNELS(st->stream_channels);
   int LM, M;
   int tf_select;
   int nbFilledBytes, nbAvailableBytes;
   int effEnd;
   int codedBands;
   int tf_sum;
   int alloc_trim;
   int pitch_index=COMBFILTER_MINPERIOD;
   opus_val16 gain1 = 0;
   int intensity=0;
   int dual_stereo=0;
   int effectiveBytes;
   opus_val16 pf_threshold;
   int dynalloc_logp;
   opus_int32 vbr_rate;
   opus_int32 total_bits;
   opus_int32 total_boost;
   opus_int32 balance;
   opus_int32 tell;
   int prefilter_tapset=0;
   int pf_on;
   int anti_collapse_rsv;
   int anti_collapse_on=0;
   int silence=0;
   ALLOC_STACK;

   if (nbCompressedBytes<2 || pcm==NULL)
     return CELT_BAD_ARG;

   frame_size *= st->upsample;
   for (LM=0;LM<=st->mode->maxLM;LM++)
      if (st->mode->shortMdctSize<<LM==frame_size)
         break;
   if (LM>st->mode->maxLM)
      return CELT_BAD_ARG;
   M=1<<LM;
   N = M*st->mode->shortMdctSize;

   prefilter_mem = st->in_mem+CC*(st->overlap);
   oldBandE = (opus_val16*)(st->in_mem+CC*(2*st->overlap+COMBFILTER_MAXPERIOD));
   oldLogE = oldBandE + CC*st->mode->nbEBands;
   oldLogE2 = oldLogE + CC*st->mode->nbEBands;

   if (enc==NULL)
   {
      tell=1;
      nbFilledBytes=0;
   } else {
      tell=ec_tell(enc);
      nbFilledBytes=(tell+4)>>3;
   }

   if (st->signalling && enc==NULL)
   {
      int tmp = (st->mode->effEBands-st->end)>>1;
      st->end = IMAX(1, st->mode->effEBands-tmp);
      compressed[0] = tmp<<5;
      compressed[0] |= LM<<3;
      compressed[0] |= (C==2)<<2;
      /* Convert "standard mode" to Opus header */
      if (st->mode->Fs==48000 && st->mode->shortMdctSize==120)
      {
         int c0 = toOpus(compressed[0]);
         if (c0<0)
            return CELT_BAD_ARG;
         compressed[0] = c0;
      }
      compressed++;
      nbCompressedBytes--;
   }

   /* Can't produce more than 1275 output bytes */
   nbCompressedBytes = IMIN(nbCompressedBytes,1275);
   nbAvailableBytes = nbCompressedBytes - nbFilledBytes;