Speed up the bit allocator interpolation
[opus.git] / libcelt / rate.c
1 /* Copyright (c) 2007-2008 CSIRO
2    Copyright (c) 2007-2009 Xiph.Org Foundation
3    Written by Jean-Marc Valin */
4 /*
5    Redistribution and use in source and binary forms, with or without
6    modification, are permitted provided that the following conditions
7    are met:
8    
9    - Redistributions of source code must retain the above copyright
10    notice, this list of conditions and the following disclaimer.
11    
12    - Redistributions in binary form must reproduce the above copyright
13    notice, this list of conditions and the following disclaimer in the
14    documentation and/or other materials provided with the distribution.
15    
16    - Neither the name of the Xiph.org Foundation nor the names of its
17    contributors may be used to endorse or promote products derived from
18    this software without specific prior written permission.
19    
20    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
21    ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
22    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
23    A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE FOUNDATION OR
24    CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
25    EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
26    PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
27    PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
28    LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
29    NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
30    SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
31 */
32
33 #ifdef HAVE_CONFIG_H
34 #include "config.h"
35 #endif
36
37 #include <math.h>
38 #include "modes.h"
39 #include "cwrs.h"
40 #include "arch.h"
41 #include "os_support.h"
42
43 #include "entcode.h"
44 #include "rate.h"
45
46
47 #ifndef STATIC_MODES
48
49 /*Determines if V(N,K) fits in a 32-bit unsigned integer.
50   N and K are themselves limited to 15 bits.*/
51 static int fits_in32(int _n, int _k)
52 {
53    static const celt_int16 maxN[15] = {
54       32767, 32767, 32767, 1476, 283, 109,  60,  40,
55        29,  24,  20,  18,  16,  14,  13};
56    static const celt_int16 maxK[15] = {
57       32767, 32767, 32767, 32767, 1172, 238,  95,  53,
58        36,  27,  22,  18,  16,  15,  13};
59    if (_n>=14)
60    {
61       if (_k>=14)
62          return 0;
63       else
64          return _n <= maxN[_k];
65    } else {
66       return _k <= maxK[_n];
67    }
68 }
69
70 void compute_pulse_cache(CELTMode *m, int LM)
71 {
72    int i;
73    int curr=0;
74    int nbEntries=0;
75    int entryN[100], entryK[100], entryI[100];
76    const celt_int16 *eBands = m->eBands;
77    PulseCache *cache = &m->cache;
78    celt_int16 *cindex;
79    unsigned char *bits;
80
81    cindex = celt_alloc(sizeof(cache->index[0])*m->nbEBands*(LM+2));
82    cache->index = cindex;
83
84    for (i=0;i<=LM+1;i++)
85    {
86       int j;
87       for (j=0;j<m->nbEBands;j++)
88       {
89          int k;
90          int N = (eBands[j+1]-eBands[j])<<i>>1;
91          cindex[i*m->nbEBands+j] = -1;
92          for (k=0;k<=i;k++)
93          {
94             int n;
95             for (n=0;n<m->nbEBands && (k!=i || n<j);n++)
96             {
97                if (N == (eBands[n+1]-eBands[n])<<k>>1)
98                {
99                   cindex[i*m->nbEBands+j] = cindex[k*m->nbEBands+n];
100                   break;
101                }
102             }
103          }
104          if (cache->index[i*m->nbEBands+j] == -1 && N!=0)
105          {
106             int K;
107             entryN[nbEntries] = N;
108             K = 0;
109             while (fits_in32(N,get_pulses(K+1)) && K<MAX_PSEUDO-1)
110                K++;
111             entryK[nbEntries] = K;
112             cindex[i*m->nbEBands+j] = curr;
113             entryI[nbEntries] = curr;
114
115             curr += K+1;
116             nbEntries++;
117          }
118       }
119    }
120    bits = celt_alloc(sizeof(unsigned char)*curr);
121    cache->bits = bits;
122    cache->size = curr;
123    for (i=0;i<nbEntries;i++)
124    {
125       int j;
126       unsigned char *ptr = bits+entryI[i];
127       celt_int16 tmp[MAX_PULSES];
128       get_required_bits(tmp, entryN[i], get_pulses(entryK[i]), BITRES);
129       for (j=1;j<=entryK[i];j++)
130          ptr[j] = tmp[get_pulses(j)]-1;
131       ptr[0] = entryK[i];
132    }
133 }
134
135 #endif /* !STATIC_MODES */
136
137
138 #define ALLOC_STEPS 6
139
140 static inline void 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)
141 {
142    int psum;
143    int lo, hi;
144    int i, j;
145    int logM;
146    int psum0;
147    const int C = CHANNELS(_C);
148    SAVE_STACK;
149
150    psum0=0;
151    for (j=start;j<end;j++)
152       psum0 += bits1[j];
153
154    logM = LM<<BITRES;
155    lo = 0;
156    hi = 1<<ALLOC_STEPS;
157    for (i=0;i<ALLOC_STEPS;i++)
158    {
159       int mid = (lo+hi)>>1;
160       psum = psum0;
161       for (j=start;j<end;j++)
162          psum += (mid*bits2[j]>>ALLOC_STEPS);
163       if (psum > (total<<BITRES))
164          hi = mid;
165       else
166          lo = mid;
167    }
168    psum = 0;
169    /*printf ("interp bisection gave %d\n", lo);*/
170    for (j=start;j<end;j++)
171    {
172       bits[j] = bits1[j] + (lo*bits2[j]>>ALLOC_STEPS);
173       psum += bits[j];
174    }
175    /* Allocate the remaining bits */
176    {
177       int left, perband;
178       left = (total<<BITRES)-psum;
179       perband = left/(end-start);
180       for (j=start;j<end;j++)
181          bits[j] += perband;
182       left = left-end*perband;
183       for (j=start;j<start+left;j++)
184          bits[j]++;
185    }
186    for (j=start;j<end;j++)
187    {
188       int N0, N, den;
189       int offset;
190       int NClogN;
191
192       N0 = m->eBands[j+1]-m->eBands[j];
193       N=N0<<LM;
194       NClogN = N*C*(m->logN[j] + logM);
195
196       /* Compensate for the extra DoF in stereo */
197       den=(C*N+ ((C==2 && N>2) ? 1 : 0));
198
199       /* Offset for the number of fine bits by log2(N)/2 + FINE_OFFSET
200          compared to their "fair share" of total/N */
201       offset = (NClogN>>1)-N*C*FINE_OFFSET;
202
203       /* N=2 is the only point that doesn't match the curve */
204       if (N==2)
205          offset += N*C<<BITRES>>2;
206
207       /* Changing the offset for allocating the second and third fine energy bit */
208       if (bits[j] + offset < den*2<<BITRES)
209          offset += NClogN>>2;
210       else if (bits[j] + offset < den*3<<BITRES)
211          offset += NClogN>>3;
212
213       /* Divide with rounding */
214       ebits[j] = (bits[j] + offset + (den<<(BITRES-1))) / (den<<BITRES);
215
216       /* If we rounded down, make it a candidate for final fine energy pass */
217       fine_priority[j] = ebits[j]*(den<<BITRES) >= bits[j]+offset;
218
219       /* For N=1, all bits go to fine energy except for a single sign bit */
220       if (N==1)
221          ebits[j] = (bits[j]/C >> BITRES)-1;
222       /* Make sure the first bit is spent on fine energy */
223       if (ebits[j] < 1)
224          ebits[j] = 1;
225
226       /* Make sure not to bust */
227       if (C*ebits[j] > (bits[j]>>BITRES))
228          ebits[j] = bits[j]/C >> BITRES;
229
230       /* More than that is useless because that's about as far as PVQ can go */
231       if (ebits[j]>7)
232          ebits[j]=7;
233
234       /* The other bits are assigned to PVQ */
235       bits[j] -= C*ebits[j]<<BITRES;
236       if (bits[j] < 0)
237          bits[j] = 0;
238    }
239    RESTORE_STACK;
240 }
241
242 void compute_allocation(const CELTMode *m, int start, int end, int *offsets, int total, int *pulses, int *ebits, int *fine_priority, int _C, int LM)
243 {
244    int lo, hi, len, j;
245    const int C = CHANNELS(_C);
246    VARDECL(int, bits1);
247    VARDECL(int, bits2);
248    SAVE_STACK;
249    
250    len = m->nbEBands;
251    ALLOC(bits1, len, int);
252    ALLOC(bits2, len, int);
253
254    lo = 0;
255    hi = m->nbAllocVectors - 1;
256    while (hi-lo != 1)
257    {
258       int psum = 0;
259       int mid = (lo+hi) >> 1;
260       for (j=start;j<end;j++)
261       {
262          int N = m->eBands[j+1]-m->eBands[j];
263          bits1[j] = ((C*N*m->allocVectors[mid*len+j]<<LM>>2) + offsets[j]);
264          psum += bits1[j];
265          /*printf ("%d ", bits[j]);*/
266       }
267       /*printf ("\n");*/
268       if (psum > (total<<BITRES))
269          hi = mid;
270       else
271          lo = mid;
272       /*printf ("lo = %d, hi = %d\n", lo, hi);*/
273    }
274    /*printf ("interp between %d and %d\n", lo, hi);*/
275    for (j=start;j<end;j++)
276    {
277       int N = m->eBands[j+1]-m->eBands[j];
278       bits1[j] = (C*N*m->allocVectors[lo*len+j]<<LM>>2);
279       bits2[j] = (C*N*m->allocVectors[hi*len+j]<<LM>>2) - bits1[j];
280       bits1[j] += offsets[j];
281    }
282    interp_bits2pulses(m, start, end, bits1, bits2, total, pulses, ebits, fine_priority, len, C, LM);
283    RESTORE_STACK;
284 }
285