Squashed commit of the following:
[opus.git] / libcelt / rate.c
index 6d98dd8..9763b6f 100644 (file)
@@ -1,5 +1,6 @@
-/* (C) 2007 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
    SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */
 
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
 #include <math.h>
 #include "modes.h"
 #include "cwrs.h"
 #include "os_support.h"
 
 #include "entcode.h"
+#include "rate.h"
 
-#define BITRES 4
-#define BITOVERFLOW 1000
 
-int log2_frac(ec_uint32 val, int frac)
-{
-   int i;
-   /* EC_ILOG() actually returns log2()+1, go figure */
-   int L = EC_ILOG(val)-1;
-   //printf ("in: %d %d ", val, L);
-   if (L>14)
-      val >>= L-14;
-   else if (L<14)
-      val <<= 14-L;
-   L <<= frac;
-   //printf ("%d\n", val);
-   for (i=0;i<frac;i++)
-   {
-      val = (val*val) >> 15;
-      //printf ("%d\n", val);
-      if (val > 16384)
-         L |= (1<<(frac-i-1));
-      else   
-         val <<= 1;
-   }
-   return L;
-}
+#ifndef STATIC_MODES
 
-int log2_frac64(ec_uint64 val, int frac)
+/*Determines if V(N,K) fits in a 32-bit unsigned integer.
+  N and K are themselves limited to 15 bits.*/
+static int fits_in32(int _n, int _k)
 {
-   int i;
-   /* EC_ILOG64() actually returns log2()+1, go figure */
-   int L = EC_ILOG64(val)-1;
-   //printf ("in: %d %d ", val, L);
-   if (L>14)
-      val >>= L-14;
-   else if (L<14)
-      val <<= 14-L;
-   L <<= frac;
-   //printf ("%d\n", val);
-   for (i=0;i<frac;i++)
+   static const celt_int16 maxN[15] = {
+      32767, 32767, 32767, 1476, 283, 109,  60,  40,
+       29,  24,  20,  18,  16,  14,  13};
+   static const celt_int16 maxK[15] = {
+      32767, 32767, 32767, 32767, 1172, 238,  95,  53,
+       36,  27,  22,  18,  16,  15,  13};
+   if (_n>=14)
    {
-      val = (val*val) >> 15;
-      //printf ("%d\n", val);
-      if (val > 16384)
-         L |= (1<<(frac-i-1));
-      else   
-         val <<= 1;
+      if (_k>=14)
+         return 0;
+      else
+         return _n <= maxN[_k];
+   } else {
+      return _k <= maxK[_n];
    }
-   return L;
 }
 
-int bits2pulses(int bits, int N)
+void compute_pulse_cache(CELTMode *m, int LM)
 {
-   int i, b, prev;
-   /* FIXME: This is terribly inefficient. Do a bisection instead
-   but be careful about overflows */
-   prev = 0;
-   i=1;
-   b = log2_frac64(ncwrs64(N, i),0);
-   while (b<bits)
-   {
-      prev=b;
-      i++;
-      b = log2_frac64(ncwrs64(N, i),0);
-   }
-   if (bits-prev < b-bits)
-      i--;
-   return i;
-}
-
+   int i;
+   int curr=0;
+   int nbEntries=0;
+   int entryN[100], entryK[100], entryI[100];
+   const celt_int16 *eBands = m->eBands;
+   PulseCache *cache = &m->cache;
+   celt_int16 *cindex;
+   unsigned char *bits;
 
-struct alloc_data {
-   int len;
-   int **bits;
-   int **rev_bits;
-};
+   cindex = celt_alloc(sizeof(cache->index[0])*m->nbEBands*(LM+2));
+   cache->index = cindex;
 
-void alloc_init(struct alloc_data *alloc, const CELTMode *m)
-{
-   int i, prevN, BC;
-   const int *eBands = m->eBands;
-   alloc->bits = celt_alloc(m->nbEBands*sizeof(int*));
-   alloc->rev_bits = celt_alloc(m->nbEBands*sizeof(int*));
-   alloc->len = m->nbEBands;
-   BC = m->nbMdctBlocks*m->nbChannels;
-   prevN = -1;
-   for (i=0;i<alloc->len;i++)
+   /* Scan for all unique band sizes */
+   for (i=0;i<=LM+1;i++)
    {
-      int N = BC*(eBands[i+1]-eBands[i]);
-      if (N == prevN)
+      int j;
+      for (j=0;j<m->nbEBands;j++)
       {
-         alloc->bits[i] = alloc->bits[i-1];
-         alloc->rev_bits[i] = alloc->rev_bits[i-1];
-      } else {
-         int j;
-         /* FIXME: We could save memory here */
-         alloc->bits[i] = celt_alloc(64*sizeof(int));
-         alloc->rev_bits[i] = celt_alloc(64*sizeof(int));
-         for (j=0;j<64;j++)
+         int k;
+         int N = (eBands[j+1]-eBands[j])<<i>>1;
+         cindex[i*m->nbEBands+j] = -1;
+         /* Find other bands that have the same size */
+         for (k=0;k<=i;k++)
+         {
+            int n;
+            for (n=0;n<m->nbEBands && (k!=i || n<j);n++)
+            {
+               if (N == (eBands[n+1]-eBands[n])<<k>>1)
+               {
+                  cindex[i*m->nbEBands+j] = cindex[k*m->nbEBands+n];
+                  break;
+               }
+            }
+         }
+         if (cache->index[i*m->nbEBands+j] == -1 && N!=0)
          {
-            alloc->bits[i][j] = log2_frac64(ncwrs64(N, j),BITRES);
-            /* We could just update rev_bits here */
-            if (alloc->bits[i][j] > (60>>BITRES))
-               break;
+            int K;
+            entryN[nbEntries] = N;
+            K = 0;
+            while (fits_in32(N,get_pulses(K+1)) && K<MAX_PSEUDO)
+               K++;
+            entryK[nbEntries] = K;
+            cindex[i*m->nbEBands+j] = curr;
+            entryI[nbEntries] = curr;
+
+            curr += K+1;
+            nbEntries++;
          }
-         for (;j<64;j++)
-            alloc->bits[i][j] = BITOVERFLOW;
-         for (j=0;j<32;j++)
-            alloc->rev_bits[i][j] = bits2pulses(j, N);
-         prevN = N;
       }
    }
-}
-
-void alloc_clear(struct alloc_data *alloc)
-{
-   int i;
-   int *prevPtr = NULL;
-   for (i=0;i<alloc->len;i++)
+   bits = celt_alloc(sizeof(unsigned char)*curr);
+   cache->bits = bits;
+   cache->size = curr;
+   /* Compute the cache for all unique sizes */
+   for (i=0;i<nbEntries;i++)
    {
-      if (alloc->bits[i] != prevPtr)
-      {
-         prevPtr = alloc->bits[i];
-         celt_free(alloc->bits[i]);
-         celt_free(alloc->rev_bits[i]);
-      }
+      int j;
+      unsigned char *ptr = bits+entryI[i];
+      celt_int16 tmp[MAX_PULSES+1];
+      get_required_bits(tmp, entryN[i], get_pulses(entryK[i]), BITRES);
+      for (j=1;j<=entryK[i];j++)
+         ptr[j] = tmp[get_pulses(j)]-1;
+      ptr[0] = entryK[i];
    }
-   celt_free(alloc->bits);
-   celt_free(alloc->rev_bits);
 }
 
+#endif /* !STATIC_MODES */
+
 
-int compute_allocation(const CELTMode *m, int *pulses)
+#define ALLOC_STEPS 6
+
+static inline int interp_bits2pulses(const CELTMode *m, int start, int end, int *bits1, int *bits2, int total, int *bits, int *ebits, int *fine_priority, int len, int _C, int LM)
 {
-   int i, N, BC, bits;
-   const int *eBands = m->eBands;
-   BC = m->nbMdctBlocks*m->nbChannels;
-   bits = 0;
-   for (i=0;i<m->nbEBands;i++)
+   int psum;
+   int lo, hi;
+   int i, j;
+   int logM;
+   const int C = CHANNELS(_C);
+   int codedBands=-1;
+   VARDECL(int, thresh);
+   SAVE_STACK;
+
+   ALLOC(thresh, len, int);
+
+   /* Threshold: don't allow any band to go below 3/8 bit/sample */
+   for (j=start;j<end;j++)
+      thresh[j] = 3*(C*(m->eBands[j+1]-m->eBands[j])<<LM<<BITRES)>>3;
+   logM = LM<<BITRES;
+   lo = 0;
+   hi = 1<<ALLOC_STEPS;
+   for (i=0;i<ALLOC_STEPS;i++)
    {
-      int q;
-      N = BC*(eBands[i+1]-eBands[i]);
-      q = pulses[i];
-      if (q<=0)
+      int mid = (lo+hi)>>1;
+      psum = 0;
+      for (j=start;j<end;j++)
       {
-         bits += log2_frac64(eBands[i] - (eBands[i+1]-eBands[i]), 8) + (1<<8);
-         q = -q;
+         int tmp = bits1[j] + (mid*bits2[j]>>ALLOC_STEPS);
+         if (tmp >= thresh[j])
+            psum += tmp;
+         else if (tmp >= 1<<BITRES)
+            psum += 1<<BITRES;
       }
-      if (q != 0)
-         bits += log2_frac64(ncwrs64(N, pulses[i]), 8);
+      if (psum > (total<<BITRES))
+         hi = mid;
+      else
+         lo = mid;
    }
-   return (bits+255)>>8;
-}
-
-
-int vec_bits2pulses(int *bands, int *bits, int *pulses, int len, int B)
-{
-   int i;
-   int sum=0;
-   for (i=0;i<len;i++)
+   psum = 0;
+   /*printf ("interp bisection gave %d\n", lo);*/
+   for (j=start;j<end;j++)
    {
-      int N = (bands[i+1]-bands[i])*B;
-      pulses[i] = bits2pulses(bits[i], N);
-      sum += log2_frac64(ncwrs(N, pulses[i]),8);
+      int tmp = bits1[j] + (lo*bits2[j]>>ALLOC_STEPS);
+      if (tmp >= thresh[j])
+      {
+         bits[j] = tmp;
+         codedBands = j;
+      } else if (tmp >= 1<<BITRES)
+         bits[j] = 1<<BITRES;
+      else
+         bits[j] = 0;
+      psum += bits[j];
    }
-   return (sum+255)>>8;
-}
-
-int interp_bits2pulses(int *bands, int *bits1, int *bits2, int total, int *pulses, int len, int B)
-{
-   int i;
-   /* FIXME: This too is terribly inefficient. We should do a bisection instead */
-   for (i=0;i<16;i++)
-   {
-      int j;
-      int bits[len];
-      for (j=0;j<len;j++)
-         bits[j] = ((16-i)*bits1[j] + i*bits2[j]) >> 4;
-      if (vec_bits2pulses(bands, bits, pulses, len, B) > total)
-         break;
+   codedBands++;
+   /* Allocate the remaining bits */
+   if (codedBands) {
+      int left, perband;
+      left = (total<<BITRES)-psum;
+      perband = left/(codedBands-start);
+      for (j=start;j<codedBands;j++)
+         bits[j] += perband;
+      left = left-codedBands*perband;
+      for (j=start;j<start+left;j++)
+         bits[j]++;
    }
-   if (i==0)
-      return -1;
-   else {
-      int j;
-      int bits[len];
-      /* Get the previous one (that didn't bust). Should rewrite that anyway */
-      i--;
-      for (j=0;j<len;j++)
-         bits[j] = ((16-i)*bits1[j] + i*bits2[j]) >> 4;      
-      return vec_bits2pulses(bands, bits, pulses, len, B);
+   for (j=start;j<end;j++)
+   {
+      int N0, N, den;
+      int offset;
+      int NClogN;
+
+      N0 = m->eBands[j+1]-m->eBands[j];
+      N=N0<<LM;
+      NClogN = N*C*(m->logN[j] + logM);
+
+      /* Compensate for the extra DoF in stereo */
+      den=(C*N+ ((C==2 && N>2) ? 1 : 0));
+
+      /* Offset for the number of fine bits by log2(N)/2 + FINE_OFFSET
+         compared to their "fair share" of total/N */
+      offset = (NClogN>>1)-N*C*FINE_OFFSET;
+
+      /* N=2 is the only point that doesn't match the curve */
+      if (N==2)
+         offset += N*C<<BITRES>>2;
+
+      /* Changing the offset for allocating the second and third fine energy bit */
+      if (bits[j] + offset < den*2<<BITRES)
+         offset += NClogN>>2;
+      else if (bits[j] + offset < den*3<<BITRES)
+         offset += NClogN>>3;
+
+      /* Divide with rounding */
+      ebits[j] = (bits[j] + offset + (den<<(BITRES-1))) / (den<<BITRES);
+
+      /* If we rounded down, make it a candidate for final fine energy pass */
+      fine_priority[j] = ebits[j]*(den<<BITRES) >= bits[j]+offset;
+
+      /* For N=1, all bits go to fine energy except for a single sign bit */
+      if (N==1)
+         ebits[j] = (bits[j]/C >> BITRES)-1;
+      /* Make sure the first bit is spent on fine energy */
+      if (ebits[j] < 1)
+         ebits[j] = 1;
+
+      /* Make sure not to bust */
+      if (C*ebits[j] > (bits[j]>>BITRES))
+         ebits[j] = bits[j]/C >> BITRES;
+
+      /* More than that is useless because that's about as far as PVQ can go */
+      if (ebits[j]>7)
+         ebits[j]=7;
+
+      /* The other bits are assigned to PVQ */
+      bits[j] -= C*ebits[j]<<BITRES;
+      if (bits[j] < 0)
+         bits[j] = 0;
    }
+   RESTORE_STACK;
+   return codedBands;
 }
 
-#if 0
-int main()
+int compute_allocation(const CELTMode *m, int start, int end, int *offsets, int alloc_trim,
+      int total, int *pulses, int *ebits, int *fine_priority, int _C, int LM)
 {
-   int i;
-   /*for(i=1;i<2000000000;i+=1738)
-   {
-      printf ("%d %d\n", i, frac_log2(i, 10));
-   }*/
-   for (i=4;i<=32;i*=2)
+   int lo, hi, len, j;
+   const int C = CHANNELS(_C);
+   int codedBands;
+   VARDECL(int, bits1);
+   VARDECL(int, bits2);
+   VARDECL(int, thresh);
+   VARDECL(int, trim_offset);
+   SAVE_STACK;
+   
+   len = m->nbEBands;
+   ALLOC(bits1, len, int);
+   ALLOC(bits2, len, int);
+   ALLOC(thresh, len, int);
+   ALLOC(trim_offset, len, int);
+
+   /* Below this threshold, we don't allocate any PVQ bits */
+   for (j=start;j<end;j++)
+      thresh[j] = 3*(C*(m->eBands[j+1]-m->eBands[j])<<LM<<BITRES)>>3;
+   /* Tilt of the allocation curve */
+   for (j=start;j<end;j++)
+      trim_offset[j] = C*(m->eBands[j+1]-m->eBands[j])*(2*alloc_trim-7)*(m->nbEBands-j-1)
+            <<(LM+BITRES)>>6;
+
+   lo = 0;
+   hi = m->nbAllocVectors - 1;
+   while (hi-lo != 1)
    {
-      int j;
-      for (j=0;j<30;j++)
+      int psum = 0;
+      int mid = (lo+hi) >> 1;
+      for (j=start;j<end;j++)
       {
-         printf ("%d %d %d\n", i, j, bits2pulses(j,i));
+         int N = m->eBands[j+1]-m->eBands[j];
+         bits1[j] = C*N*m->allocVectors[mid*len+j]<<LM>>2;
+         if (bits1[j] > 0)
+            bits1[j] += trim_offset[j];
+         if (bits1[j] < 0)
+            bits1[j] = 0;
+         bits1[j] += offsets[j];
+         if (bits1[j] >= thresh[j])
+            psum += bits1[j];
+         else if (bits1[j] >= 1<<BITRES)
+            psum += 1<<BITRES;
+
+         /*printf ("%d ", bits[j]);*/
       }
+      /*printf ("\n");*/
+      if (psum > (total<<BITRES))
+         hi = mid;
+      else
+         lo = mid;
+      /*printf ("lo = %d, hi = %d\n", lo, hi);*/
    }
-   return 0;
-}
-#endif
-#if 0
-int main()
-{
-   int i;
-   int bits[18] = {10, 9, 9, 8, 8, 8, 8, 8, 8, 8, 9, 10, 8, 9, 10, 11, 6, 7};
-   int bits1[18] = {8, 7, 7, 6, 6, 6, 5, 4, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5};
-   int bits2[18] = {15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15};
-   int bank[20] = {0,  4,  8, 12, 16, 20, 24, 28, 32, 38, 44, 52, 62, 74, 90,112,142,182, 232,256};
-   int pulses[18];
-   struct alloc_data alloc;
-   
-   alloc_init(&alloc, celt_mode0);
-   int b = vec_bits2pulses(bank, bits, pulses, 18, 1);
-   printf ("total: %d bits\n", b);
-   for (i=0;i<18;i++)
-      printf ("%d ", pulses[i]);
-   printf ("\n");
-   b = interp_bits2pulses(bank, bits1, bits2, 160, pulses, 18, 1);
-   printf ("total: %d bits\n", b);
-   for (i=0;i<18;i++)
-      printf ("%d ", pulses[i]);
-   printf ("\n");
-
-   alloc_clear(&alloc);
-   return 0;
+   /*printf ("interp between %d and %d\n", lo, hi);*/
+   for (j=start;j<end;j++)
+   {
+      int N = m->eBands[j+1]-m->eBands[j];
+      bits1[j] = (C*N*m->allocVectors[lo*len+j]<<LM>>2);
+      bits2[j] = (C*N*m->allocVectors[hi*len+j]<<LM>>2) - bits1[j];
+      if (bits1[j] > 0)
+         bits1[j] += trim_offset[j];
+      if (bits1[j] < 0)
+         bits1[j] = 0;
+      bits1[j] += offsets[j];
+   }
+   codedBands = interp_bits2pulses(m, start, end, bits1, bits2, total, pulses, ebits, fine_priority, len, C, LM);
+   RESTORE_STACK;
+   return codedBands;
 }
-#endif
+