Some work towards being able to encode a 48 kHz stream from 32 kHz audio (incomplete)
[opus.git] / libcelt / quant_bands.c
index 0e2b065..af2b7e1 100644 (file)
@@ -1,5 +1,6 @@
-/* (C) 2007-2008 Jean-Marc Valin, CSIRO
-*/
+/* Copyright (c) 2007-2008 CSIRO
+   Copyright (c) 2007-2009 Xiph.Org Foundation
+   Written by Jean-Marc Valin */
 /*
    Redistribution and use in source and binary forms, with or without
    modification, are permitted provided that the following conditions
 #include "mathops.h"
 #include "stack_alloc.h"
 
-#ifdef FIXED_POINT
-const celt_word16_t eMeans[24] = {11520, -2048, -3072, -640, 256, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
-#else
-const celt_word16_t eMeans[24] = {45.f, -8.f, -12.f, -2.5f, 1.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f};
-#endif
+#define E_MEANS_SIZE (5)
 
-/*const int frac[24] = {4, 3, 3, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2};*/
-/*const int frac[24] = {8, 6, 5, 4, 3, 3, 3, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2};*/
+const celt_word16 eMeans[E_MEANS_SIZE] = {QCONST16(7.5f,DB_SHIFT), -QCONST16(1.33f,DB_SHIFT), -QCONST16(2.f,DB_SHIFT), -QCONST16(0.42f,DB_SHIFT), QCONST16(0.17f,DB_SHIFT)};
 
-static void compute_fine_allocation(const CELTMode *m, celt_int16_t *bits, int budget)
+/* FIXME: Implement for stereo */
+int intra_decision(celt_word16 *eBands, celt_word16 *oldEBands, int len)
 {
-   int i,j;
-   int len;
-   len = m->nbEBands;
-   for (i=0;i<m->nbAllocVectors;i++)
-   {
-      if (m->energy_alloc[i*(len+1)+len] > budget)
-         break;
-   }
-   if (i==0)
-   {
-      for (j=0;j<len;j++)
-         bits[j] = 0;
-   } else {
-      for (j=0;j<len;j++)
-         bits[j] = m->energy_alloc[(i-1)*(len+1)+j];
-      budget -= m->energy_alloc[(i-1)*(len+1)+len];
-   }
-   if (i<m->nbAllocVectors)
+   int i;
+   celt_word32 dist = 0;
+   for (i=0;i<len;i++)
    {
-      j=0;
-      while (budget>0)
-      {
-         if (m->energy_alloc[i*(len+1)+j]>bits[j])
-         {
-            bits[j]++;
-            budget--;
-         }
-         j++;
-         if (j>=len)
-            j=0;
-      }
+      celt_word16 d = SUB16(eBands[i], oldEBands[i]);
+      dist = MAC16_16(dist, d,d);
    }
-   
-   /*for (j=0;j<len;j++)
-      printf ("%d ", bits[j]);
-   printf ("\n");*/
-}
-
-#ifdef FIXED_POINT
-static inline celt_ener_t dB2Amp(celt_ener_t dB)
-{
-   celt_ener_t amp;
-   amp = PSHR32(celt_exp2(MULT16_16_Q14(21771,dB)),2)-QCONST16(.3f, 14);
-   if (amp < 0)
-      amp = 0;
-   return amp;
-}
-
-#define DBofTWO 24661
-static inline celt_word16_t amp2dB(celt_ener_t amp)
-{
-   /* equivalent to return 6.0207*log2(.3+amp) */
-   return ROUND16(MULT16_16(24661,celt_log2(ADD32(QCONST32(.3f,14),amp))),12);
-   /* return DB_SCALING*20*log10(.3+ENER_SCALING_1*amp); */
-}
-#else
-static inline celt_ener_t dB2Amp(celt_ener_t dB)
-{
-   celt_ener_t amp;
-   amp = pow(10, .05*dB)-.3;
-   if (amp < 0)
-      amp = 0;
-   return amp;
-}
-static inline celt_word16_t amp2dB(celt_ener_t amp)
-{
-   return 20*log10(.3+amp);
+   return SHR32(dist,2*DB_SHIFT) > 2*len;
 }
-#endif
-
-static const celt_word16_t base_resolution = QCONST16(6.f,8);
 
 int *quant_prob_alloc(const CELTMode *m)
 {
    int i;
    int *prob;
-   prob = celt_alloc(2*m->nbEBands*sizeof(int));
+   prob = celt_alloc(4*m->nbEBands*sizeof(int));
+   if (prob==NULL)
+     return NULL;
    for (i=0;i<m->nbEBands;i++)
    {
       prob[2*i] = 6000-i*200;
       prob[2*i+1] = ec_laplace_get_start_freq(prob[2*i]);
    }
+   for (i=0;i<m->nbEBands;i++)
+   {
+      prob[2*m->nbEBands+2*i] = 9000-i*240;
+      prob[2*m->nbEBands+2*i+1] = ec_laplace_get_start_freq(prob[2*m->nbEBands+2*i]);
+   }
    return prob;
 }
 
@@ -142,226 +84,242 @@ void quant_prob_free(int *freq)
    celt_free(freq);
 }
 
-static void quant_energy_mono(const CELTMode *m, celt_ener_t *eBands, celt_word16_t *oldEBands, unsigned budget, int *prob, ec_enc *enc)
+unsigned quant_coarse_energy(const CELTMode *m, int start, int end, const celt_word16 *eBands, celt_word16 *oldEBands, int budget, int intra, int *prob, celt_word16 *error, ec_enc *enc, int _C, celt_word16 max_decay)
 {
-   int i;
-   unsigned bits;
-   celt_word16_t prev = 0;
-   celt_word16_t coef = m->ePredCoef;
-   celt_word16_t beta;
-   VARDECL(celt_word16_t, error);
-   VARDECL(celt_int16_t, fine_quant);
-   SAVE_STACK;
-   /* The .7 is a heuristic */
-   beta = MULT16_16_Q15(QCONST16(.8f,15),coef);
-   
-   ALLOC(error, m->nbEBands, celt_word16_t);
-   ALLOC(fine_quant, m->nbEBands, celt_int16_t);
-   bits = ec_enc_tell(enc, 0);
+   int i, c;
+   unsigned bits_used = 0;
+   celt_word32 prev[2] = {0,0};
+   celt_word16 coef = m->ePredCoef;
+   celt_word16 beta;
+   const int C = CHANNELS(_C);
+
+   if (intra)
+   {
+      coef = 0;
+      prob += 2*m->nbEBands;
+   }
+   /* The .8 is a heuristic */
+   beta = MULT16_16_P15(QCONST16(.8f,15),coef);
+
    /* Encode at a fixed coarse resolution */
-   for (i=0;i<m->nbEBands;i++)
+   for (i=start;i<end;i++)
    {
-      int qi;
-      celt_word16_t q;   /* dB */
-      celt_word16_t x;   /* dB */
-      celt_word16_t f;   /* Q8 */
-      celt_word16_t mean = MULT16_16_Q15(Q15ONE-coef,eMeans[i]);
-      x = amp2dB(eBands[i]);
-      f = EXTRACT16(celt_div(SHL32(EXTEND32(x-mean-MULT16_16_Q15(coef,oldEBands[i])-prev),8),base_resolution));
+      c=0;
+      do {
+         int qi;
+         celt_word16 q;
+         celt_word16 x;
+         celt_word32 f;
+         celt_word32 mean =  (i-start < E_MEANS_SIZE) ? SUB32(SHL32(EXTEND32(eMeans[i-start]),15), MULT16_16(coef,eMeans[i-start])) : 0;
+         x = eBands[i+c*m->nbEBands];
 #ifdef FIXED_POINT
-      /* Rounding to nearest integer here is really important! */
-      qi = (f+128)>>8;
+         f = SHL32(EXTEND32(x),15)-mean -MULT16_16(coef,oldEBands[i+c*m->nbEBands])-prev[c];
+         /* Rounding to nearest integer here is really important! */
+         qi = (f+QCONST32(.5,DB_SHIFT+15))>>(DB_SHIFT+15);
 #else
-      qi = (int)floor(.5+f);
+         f = x-mean-coef*oldEBands[i+c*m->nbEBands]-prev[c];
+         /* Rounding to nearest integer here is really important! */
+         qi = (int)floor(.5f+f);
 #endif
-      /* If we don't have enough bits to encode all the energy, just assume something safe.
-         We allow slightly busting the budget here */
-      if (ec_enc_tell(enc, 0) - bits > budget+16)
-         qi = -1;
-      else
-         ec_laplace_encode_start(enc, &qi, prob[2*i], prob[2*i+1]);
-      q = qi*base_resolution;
-      error[i] = f - SHL16(qi,8);
-      
-      oldEBands[i] = mean+MULT16_16_Q15(coef,oldEBands[i])+prev+q;
-      if (oldEBands[i] < -QCONST16(12.f,8))
-         oldEBands[i] = -QCONST16(12.f,8);
-      prev = mean+prev+MULT16_16_Q15(Q15ONE-beta,q);
+         if (qi < 0 && x < oldEBands[i+c*m->nbEBands]-max_decay)
+         {
+            qi += SHR16(oldEBands[i+c*m->nbEBands]-max_decay-x, DB_SHIFT);
+            if (qi > 0)
+               qi = 0;
+         }
+         /* If we don't have enough bits to encode all the energy, just assume something safe.
+            We allow slightly busting the budget here */
+         bits_used=ec_enc_tell(enc, 0);
+         if (bits_used > budget)
+         {
+            qi = -1;
+            error[i+c*m->nbEBands] = QCONST16(.5f,DB_SHIFT);
+         } else {
+            ec_laplace_encode_start(enc, &qi, prob[2*i], prob[2*i+1]);
+            error[i+c*m->nbEBands] = PSHR32(f,15) - SHL16(qi,DB_SHIFT);
+         }
+         q = SHL16(qi,DB_SHIFT);
+         
+         oldEBands[i+c*m->nbEBands] = PSHR32(MULT16_16(coef,oldEBands[i+c*m->nbEBands]) + mean + prev[c] + SHL32(EXTEND32(q),15), 15);
+         prev[c] = mean + prev[c] + SHL32(EXTEND32(q),15) - MULT16_16(beta,q);
+      } while (++c < C);
    }
-   
-   compute_fine_allocation(m, fine_quant, budget-(ec_enc_tell(enc, 0)-bits));
+   return bits_used;
+}
+
+void quant_fine_energy(const CELTMode *m, int start, int end, celt_ener *eBands, celt_word16 *oldEBands, celt_word16 *error, int *fine_quant, ec_enc *enc, int _C)
+{
+   int i, c;
+   const int C = CHANNELS(_C);
 
    /* Encode finer resolution */
-   for (i=0;i<m->nbEBands;i++)
+   for (i=start;i<end;i++)
    {
-      int q2;
-      celt_int16_t frac = 1<<fine_quant[i];
-      celt_word16_t offset = (error[i]+QCONST16(.5f,8))*frac;
+      celt_int16 frac = 1<<fine_quant[i];
       if (fine_quant[i] <= 0)
          continue;
+      c=0;
+      do {
+         int q2;
+         celt_word16 offset;
+#ifdef FIXED_POINT
+         /* Has to be without rounding */
+         q2 = (error[i+c*m->nbEBands]+QCONST16(.5f,DB_SHIFT))>>(DB_SHIFT-fine_quant[i]);
+#else
+         q2 = (int)floor((error[i+c*m->nbEBands]+.5f)*frac);
+#endif
+         if (q2 > frac-1)
+            q2 = frac-1;
+         if (q2<0)
+            q2 = 0;
+         ec_enc_bits(enc, q2, fine_quant[i]);
 #ifdef FIXED_POINT
-      /* Has to be without rounding */
-      q2 = offset>>8;
+         offset = SUB16(SHR16(SHL16(q2,DB_SHIFT)+QCONST16(.5,DB_SHIFT),fine_quant[i]),QCONST16(.5f,DB_SHIFT));
 #else
-      q2 = (int)floor(offset);
+         offset = (q2+.5f)*(1<<(14-fine_quant[i]))*(1.f/16384) - .5f;
 #endif
-      if (q2 > frac-1)
-         q2 = frac-1;
-      ec_enc_bits(enc, q2, fine_quant[i]);
-      offset = EXTRACT16(celt_div(SHL16(q2,8)+QCONST16(.5,8),frac)-QCONST16(.5f,8));
-      oldEBands[i] += PSHR32(MULT16_16(DB_SCALING*6,offset),8);
-      /*printf ("%f ", error[i] - offset);*/
+         oldEBands[i+c*m->nbEBands] += offset;
+         error[i+c*m->nbEBands] -= offset;
+         /*printf ("%f ", error[i] - offset);*/
+      } while (++c < C);
    }
-   for (i=0;i<m->nbEBands;i++)
+}
+
+void quant_energy_finalise(const CELTMode *m, int start, int end, celt_ener *eBands, celt_word16 *oldEBands, celt_word16 *error, int *fine_quant, int *fine_priority, int bits_left, ec_enc *enc, int _C)
+{
+   int i, prio, c;
+   const int C = CHANNELS(_C);
+
+   /* Use up the remaining bits */
+   for (prio=0;prio<2;prio++)
    {
-      eBands[i] = dB2Amp(oldEBands[i]);
+      for (i=start;i<end && bits_left>=C ;i++)
+      {
+         if (fine_quant[i] >= 7 || fine_priority[i]!=prio)
+            continue;
+         c=0;
+         do {
+            int q2;
+            celt_word16 offset;
+            q2 = error[i+c*m->nbEBands]<0 ? 0 : 1;
+            ec_enc_bits(enc, q2, 1);
+#ifdef FIXED_POINT
+            offset = SHR16(SHL16(q2,DB_SHIFT)-QCONST16(.5,DB_SHIFT),fine_quant[i]+1);
+#else
+            offset = (q2-.5f)*(1<<(14-fine_quant[i]-1))*(1.f/16384);
+#endif
+            oldEBands[i+c*m->nbEBands] += offset;
+            bits_left--;
+         } while (++c < C);
+      }
    }
-   /*printf ("%d\n", ec_enc_tell(enc, 0)-9);*/
-
-   /*printf ("\n");*/
-   RESTORE_STACK;
+   c=0;
+   do {
+      for (i=start;i<m->nbEBands;i++)
+      {
+         eBands[i+c*m->nbEBands] = log2Amp(oldEBands[i+c*m->nbEBands]);
+         if (oldEBands[i+c*m->nbEBands] < -QCONST16(7.f,DB_SHIFT))
+            oldEBands[i+c*m->nbEBands] = -QCONST16(7.f,DB_SHIFT);
+      }
+   } while (++c < C);
 }
 
-static void unquant_energy_mono(const CELTMode *m, celt_ener_t *eBands, celt_word16_t *oldEBands, unsigned budget, int *prob, ec_dec *dec)
+void unquant_coarse_energy(const CELTMode *m, int start, int end, celt_ener *eBands, celt_word16 *oldEBands, int budget, int intra, int *prob, ec_dec *dec, int _C)
 {
-   int i;
-   unsigned bits;
-   celt_word16_t prev = 0;
-   celt_word16_t coef = m->ePredCoef;
-   /* The .7 is a heuristic */
-   VARDECL(celt_int16_t, fine_quant);
-   celt_word16_t beta = MULT16_16_Q15(QCONST16(.8f,15),coef);
-   SAVE_STACK;
-   
-   ALLOC(fine_quant, m->nbEBands, celt_int16_t);
-   bits = ec_dec_tell(dec, 0);
-   
+   int i, c;
+   celt_word32 prev[2] = {0, 0};
+   celt_word16 coef = m->ePredCoef;
+   celt_word16 beta;
+   const int C = CHANNELS(_C);
+
+   if (intra)
+   {
+      coef = 0;
+      prob += 2*m->nbEBands;
+   }
+   /* The .8 is a heuristic */
+   beta = MULT16_16_P15(QCONST16(.8f,15),coef);
+
    /* Decode at a fixed coarse resolution */
-   for (i=0;i<m->nbEBands;i++)
+   for (i=start;i<end;i++)
    {
-      int qi;
-      celt_word16_t q;
-      celt_word16_t mean = MULT16_16_Q15(Q15ONE-coef,eMeans[i]);
-      /* If we didn't have enough bits to encode all the energy, just assume something safe.
-         We allow slightly busting the budget here */
-      if (ec_dec_tell(dec, 0) - bits > budget+16)
-         qi = -1;
-      else
-         qi = ec_laplace_decode_start(dec, prob[2*i], prob[2*i+1]);
-      q = qi*base_resolution;
-      
-      oldEBands[i] = mean+MULT16_16_Q15(coef,oldEBands[i])+prev+q;
-      if (oldEBands[i] < -QCONST16(12.f,8))
-         oldEBands[i] = -QCONST16(12.f,8);
-      
-      prev = mean+prev+MULT16_16_Q15(Q15ONE-beta,q);
+      c=0;
+      do {
+         int qi;
+         celt_word16 q;
+         celt_word32 mean =  (i-start < E_MEANS_SIZE) ? SUB32(SHL32(EXTEND32(eMeans[i-start]),15), MULT16_16(coef,eMeans[i-start])) : 0;
+         /* If we didn't have enough bits to encode all the energy, just assume something safe.
+            We allow slightly busting the budget here */
+         if (ec_dec_tell(dec, 0) > budget)
+            qi = -1;
+         else
+            qi = ec_laplace_decode_start(dec, prob[2*i], prob[2*i+1]);
+         q = SHL16(qi,DB_SHIFT);
+
+         oldEBands[i+c*m->nbEBands] = PSHR32(MULT16_16(coef,oldEBands[i+c*m->nbEBands]) + mean + prev[c] + SHL32(EXTEND32(q),15), 15);
+         prev[c] = mean + prev[c] + SHL32(EXTEND32(q),15) - MULT16_16(beta,q);
+      } while (++c < C);
    }
-   
-   compute_fine_allocation(m, fine_quant, budget-(ec_dec_tell(dec, 0)-bits));
+}
 
+void unquant_fine_energy(const CELTMode *m, int start, int end, celt_ener *eBands, celt_word16 *oldEBands, int *fine_quant, ec_dec *dec, int _C)
+{
+   int i, c;
+   const int C = CHANNELS(_C);
    /* Decode finer resolution */
-   for (i=0;i<m->nbEBands;i++)
+   for (i=start;i<end;i++)
    {
-      int q2;
-      celt_int16_t frac = 1<<fine_quant[i];
-      celt_word16_t offset;
       if (fine_quant[i] <= 0)
          continue;
-      q2 = ec_dec_bits(dec, fine_quant[i]);
-      offset = EXTRACT16(celt_div(SHL16(q2,8)+QCONST16(.5,8),frac)-QCONST16(.5f,8));
-      oldEBands[i] += PSHR32(MULT16_16(DB_SCALING*6,offset),8);
-   }
-   for (i=0;i<m->nbEBands;i++)
-   {
-      eBands[i] = dB2Amp(oldEBands[i]);
+      c=0; 
+      do {
+         int q2;
+         celt_word16 offset;
+         q2 = ec_dec_bits(dec, fine_quant[i]);
+#ifdef FIXED_POINT
+         offset = SUB16(SHR16(SHL16(q2,DB_SHIFT)+QCONST16(.5,DB_SHIFT),fine_quant[i]),QCONST16(.5f,DB_SHIFT));
+#else
+         offset = (q2+.5f)*(1<<(14-fine_quant[i]))*(1.f/16384) - .5f;
+#endif
+         oldEBands[i+c*m->nbEBands] += offset;
+      } while (++c < C);
    }
-   RESTORE_STACK;
-   /*printf ("\n");*/
 }
 
-
-
-void quant_energy(const CELTMode *m, celt_ener_t *eBands, celt_word16_t *oldEBands, int budget, int *prob, ec_enc *enc)
+void unquant_energy_finalise(const CELTMode *m, int start, int end, celt_ener *eBands, celt_word16 *oldEBands, int *fine_quant,  int *fine_priority, int bits_left, ec_dec *dec, int _C)
 {
-   int C;
-   SAVE_STACK;
-   
-   C = m->nbChannels;
+   int i, prio, c;
+   const int C = CHANNELS(_C);
 
-   if (C==1)
-      quant_energy_mono(m, eBands, oldEBands, budget, prob, enc);
-   else 
-#if 1
+   /* Use up the remaining bits */
+   for (prio=0;prio<2;prio++)
    {
-      int c;
-      VARDECL(celt_ener_t, E);
-      ALLOC(E, m->nbEBands, celt_ener_t);
-      for (c=0;c<C;c++)
+      for (i=start;i<end && bits_left>=C ;i++)
       {
-         int i;
-         for (i=0;i<m->nbEBands;i++)
-            E[i] = eBands[C*i+c];
-         quant_energy_mono(m, E, oldEBands+c*m->nbEBands, budget/C, prob, enc);
-         for (i=0;i<m->nbEBands;i++)
-            eBands[C*i+c] = E[i];
-      }
-   }
+         if (fine_quant[i] >= 7 || fine_priority[i]!=prio)
+            continue;
+         c=0;
+         do {
+            int q2;
+            celt_word16 offset;
+            q2 = ec_dec_bits(dec, 1);
+#ifdef FIXED_POINT
+            offset = SHR16(SHL16(q2,DB_SHIFT)-QCONST16(.5,DB_SHIFT),fine_quant[i]+1);
 #else
-      if (C==2)
-   {
-      int i;
-      int NB = m->nbEBands;
-      celt_ener_t mid[NB];
-      celt_ener_t side[NB];
-      for (i=0;i<NB;i++)
-      {
-         //left = eBands[C*i];
-         //right = eBands[C*i+1];
-         mid[i] = ENER_SCALING_1*sqrt(eBands[C*i]*eBands[C*i] + eBands[C*i+1]*eBands[C*i+1]);
-         side[i] = 20*log10((ENER_SCALING_1*eBands[2*i]+.3)/(ENER_SCALING_1*eBands[2*i+1]+.3));
-         //printf ("%f %f ", mid[i], side[i]);
-      }
-      //printf ("\n");
-      quant_energy_mono(m, mid, oldEBands, enc);
-      for (i=0;i<NB;i++)
-         side[i] = pow(10.f,floor(.5f+side[i])/10.f);
-         
-      //quant_energy_side(m, side, oldEBands+NB, enc);
-      for (i=0;i<NB;i++)
-      {
-         eBands[C*i] = ENER_SCALING*mid[i]*sqrt(side[i]/(1.f+side[i]));
-         eBands[C*i+1] = ENER_SCALING*mid[i]*sqrt(1.f/(1.f+side[i]));
-         //printf ("%f %f ", mid[i], side[i]);
+            offset = (q2-.5f)*(1<<(14-fine_quant[i]-1))*(1.f/16384);
+#endif
+            oldEBands[i+c*m->nbEBands] += offset;
+            bits_left--;
+         } while (++c < C);
       }
-
-   } else {
-      celt_fatal("more than 2 channels not supported");
    }
-#endif
-   RESTORE_STACK;
-}
-
-
-
-void unquant_energy(const CELTMode *m, celt_ener_t *eBands, celt_word16_t *oldEBands, int budget, int *prob, ec_dec *dec)
-{
-   int C;   
-   SAVE_STACK;
-   C = m->nbChannels;
-
-   if (C==1)
-      unquant_energy_mono(m, eBands, oldEBands, budget, prob, dec);
-   else {
-      int c;
-      VARDECL(celt_ener_t, E);
-      ALLOC(E, m->nbEBands, celt_ener_t);
-      for (c=0;c<C;c++)
+   c=0;
+   do {
+      for (i=start;i<m->nbEBands;i++)
       {
-         int i;
-         unquant_energy_mono(m, E, oldEBands+c*m->nbEBands, budget/C, prob, dec);
-         for (i=0;i<m->nbEBands;i++)
-            eBands[C*i+c] = E[i];
+         eBands[i+c*m->nbEBands] = log2Amp(oldEBands[i+c*m->nbEBands]);
+         if (oldEBands[i+c*m->nbEBands] < -QCONST16(7.f,DB_SHIFT))
+            oldEBands[i+c*m->nbEBands] = -QCONST16(7.f,DB_SHIFT);
       }
-   }
-   RESTORE_STACK;
+   } while (++c < C);
 }