Preventing negative bit allocation
[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    /* Scan for all unique band sizes */
85    for (i=0;i<=LM+1;i++)
86    {
87       int j;
88       for (j=0;j<m->nbEBands;j++)
89       {
90          int k;
91          int N = (eBands[j+1]-eBands[j])<<i>>1;
92          cindex[i*m->nbEBands+j] = -1;
93          /* Find other bands that have the same size */
94          for (k=0;k<=i;k++)
95          {
96             int n;
97             for (n=0;n<m->nbEBands && (k!=i || n<j);n++)
98             {
99                if (N == (eBands[n+1]-eBands[n])<<k>>1)
100                {
101                   cindex[i*m->nbEBands+j] = cindex[k*m->nbEBands+n];
102                   break;
103                }
104             }
105          }
106          if (cache->index[i*m->nbEBands+j] == -1 && N!=0)
107          {
108             int K;
109             entryN[nbEntries] = N;
110             K = 0;
111             while (fits_in32(N,get_pulses(K+1)) && K<MAX_PSEUDO)
112                K++;
113             entryK[nbEntries] = K;
114             cindex[i*m->nbEBands+j] = curr;
115             entryI[nbEntries] = curr;
116
117             curr += K+1;
118             nbEntries++;
119          }
120       }
121    }
122    bits = celt_alloc(sizeof(unsigned char)*curr);
123    cache->bits = bits;
124    cache->size = curr;
125    /* Compute the cache for all unique sizes */
126    for (i=0;i<nbEntries;i++)
127    {
128       int j;
129       unsigned char *ptr = bits+entryI[i];
130       celt_int16 tmp[MAX_PULSES+1];
131       get_required_bits(tmp, entryN[i], get_pulses(entryK[i]), BITRES);
132       for (j=1;j<=entryK[i];j++)
133          ptr[j] = tmp[get_pulses(j)]-1;
134       ptr[0] = entryK[i];
135    }
136 }
137
138 #endif /* !STATIC_MODES */
139
140
141 #define ALLOC_STEPS 6
142
143 static inline int interp_bits2pulses(const CELTMode *m, int start, int end,
144       int *bits1, int *bits2, const int *thresh, int total, int *bits,
145       int *ebits, int *fine_priority, int len, int _C, int LM, int *skip, int prev)
146 {
147    int psum;
148    int lo, hi;
149    int i, j;
150    int logM;
151    const int C = CHANNELS(_C);
152    int codedBands=-1;
153    SAVE_STACK;
154
155    logM = LM<<BITRES;
156    lo = 0;
157    hi = 1<<ALLOC_STEPS;
158    for (i=0;i<ALLOC_STEPS;i++)
159    {
160       int mid = (lo+hi)>>1;
161       psum = 0;
162       for (j=start;j<end;j++)
163       {
164          int tmp = bits1[j] + (mid*bits2[j]>>ALLOC_STEPS);
165          /* Don't allocate more than we can actually use */
166          if (tmp >= thresh[j])
167             psum += tmp;
168          else if (tmp >= C<<BITRES)
169             psum += C<<BITRES;
170       }
171       if (psum > (total<<BITRES))
172          hi = mid;
173       else
174          lo = mid;
175    }
176    psum = 0;
177    /*printf ("interp bisection gave %d\n", lo);*/
178    for (j=start;j<end;j++)
179    {
180       int tmp = bits1[j] + (lo*bits2[j]>>ALLOC_STEPS);
181       if (tmp >= thresh[j])
182       {
183          bits[j] = tmp;
184          codedBands = j;
185       } else if (tmp >= C<<BITRES)
186          bits[j] = C<<BITRES;
187       else
188          bits[j] = 0;
189       /* Don't allocate more than we can actually use */
190       if (bits[j] > 64*C<<BITRES<<LM)
191          bits[j] = 64*C<<BITRES<<LM;
192       psum += bits[j];
193    }
194    codedBands++;
195    if (*skip==-1)
196    {
197       *skip=0;
198       for (j=codedBands-1;j>=0;j--)
199       {
200          if ((bits[j] > (7*(m->eBands[j+1]-m->eBands[j])<<LM<<BITRES)>>4 && j<prev)
201                || (bits[j] > (9*(m->eBands[j+1]-m->eBands[j])<<LM<<BITRES)>>4))
202             break;
203          else
204             (*skip)++;
205       }
206       *skip = IMIN(*skip, codedBands-start-1);
207    }
208    for (i=0;i<*skip;i++)
209    {
210       if (bits[codedBands-1] >= C<<BITRES)
211       {
212          psum = psum - bits[codedBands-1] + ((C+1)<<BITRES);
213          bits[codedBands-1] = C<<BITRES;
214       } else {
215          psum = psum - bits[codedBands-1];
216          bits[codedBands-1] = 0;
217       }
218       codedBands--;
219    }
220    /* Allocate the remaining bits */
221    if (codedBands) {
222       int left, perband;
223       left = (total<<BITRES)-psum;
224       perband = left/(m->eBands[codedBands]-m->eBands[start]);
225       for (j=start;j<codedBands;j++)
226          bits[j] += perband*(m->eBands[j+1]-m->eBands[j]);
227       left = left-(m->eBands[codedBands]-m->eBands[start])*perband;
228       for (j=start;j<codedBands;j++)
229       {
230          int tmp = IMIN(left, m->eBands[j+1]-m->eBands[j]);
231          bits[j] += tmp;
232          left -= tmp;
233       }
234    }
235    /*for (j=0;j<end;j++)printf("%d ", bits[j]);printf("\n");*/
236    for (j=start;j<end;j++)
237    {
238       int N0, N, den;
239       int offset;
240       int NClogN;
241
242       N0 = m->eBands[j+1]-m->eBands[j];
243       N=N0<<LM;
244       NClogN = N*C*(m->logN[j] + logM);
245
246       /* Compensate for the extra DoF in stereo */
247       den=(C*N+ ((C==2 && N>2) ? 1 : 0));
248
249       /* Offset for the number of fine bits by log2(N)/2 + FINE_OFFSET
250          compared to their "fair share" of total/N */
251       offset = (NClogN>>1)-N*C*FINE_OFFSET;
252
253       /* N=2 is the only point that doesn't match the curve */
254       if (N==2)
255          offset += N*C<<BITRES>>2;
256
257       /* Changing the offset for allocating the second and third fine energy bit */
258       if (bits[j] + offset < den*2<<BITRES)
259          offset += NClogN>>2;
260       else if (bits[j] + offset < den*3<<BITRES)
261          offset += NClogN>>3;
262
263       /* Divide with rounding */
264       ebits[j] = (bits[j] + offset + (den<<(BITRES-1))) / (den<<BITRES);
265
266       /* If we rounded down, make it a candidate for final fine energy pass */
267       fine_priority[j] = ebits[j]*(den<<BITRES) >= bits[j]+offset;
268
269       /* For N=1, all bits go to fine energy except for a single sign bit */
270       if (N==1)
271       {
272          ebits[j] = IMAX(0,(bits[j]/C >> BITRES)-1);
273          fine_priority[j] = (ebits[j]+1)*C<<BITRES >= bits[j];
274       }
275       /* Make sure not to bust */
276       if (C*ebits[j] > (bits[j]>>BITRES))
277          ebits[j] = bits[j]/C >> BITRES;
278
279       /* More than that is useless because that's about as far as PVQ can go */
280       if (ebits[j]>7)
281          ebits[j]=7;
282
283       /* The other bits are assigned to PVQ */
284       bits[j] -= C*ebits[j]<<BITRES;
285       celt_assert(bits[j] >= 0);
286    }
287    RESTORE_STACK;
288    return codedBands;
289 }
290
291 int compute_allocation(const CELTMode *m, int start, int end, int *offsets, int alloc_trim,
292       int total, int *pulses, int *ebits, int *fine_priority, int _C, int LM, int *skip, int prev)
293 {
294    int lo, hi, len, j;
295    const int C = CHANNELS(_C);
296    int codedBands;
297    VARDECL(int, bits1);
298    VARDECL(int, bits2);
299    VARDECL(int, thresh);
300    VARDECL(int, trim_offset);
301    SAVE_STACK;
302    
303    len = m->nbEBands;
304    ALLOC(bits1, len, int);
305    ALLOC(bits2, len, int);
306    ALLOC(thresh, len, int);
307    ALLOC(trim_offset, len, int);
308
309    /* Below this threshold, we don't allocate any PVQ bits */
310    for (j=start;j<end;j++)
311       thresh[j] = ((C+1)<<BITRES) + (3*(m->eBands[j+1]-m->eBands[j])<<LM<<BITRES)>>4;
312    /* Tilt of the allocation curve */
313    for (j=start;j<end;j++)
314       trim_offset[j] = C*(m->eBands[j+1]-m->eBands[j])*(alloc_trim-5-LM)*(m->nbEBands-j-1)
315             <<(LM+BITRES)>>6;
316
317    lo = 0;
318    hi = m->nbAllocVectors - 1;
319    while (hi-lo != 1)
320    {
321       int psum = 0;
322       int mid = (lo+hi) >> 1;
323       for (j=start;j<end;j++)
324       {
325          int N = m->eBands[j+1]-m->eBands[j];
326          bits1[j] = C*N*m->allocVectors[mid*len+j]<<LM>>2;
327          if (bits1[j] > 0)
328             bits1[j] += trim_offset[j];
329          if (bits1[j] < 0)
330             bits1[j] = 0;
331          bits1[j] += offsets[j];
332          if (bits1[j] >= thresh[j])
333             psum += bits1[j];
334          else if (bits1[j] >= C<<BITRES)
335             psum += C<<BITRES;
336
337          /*printf ("%d ", bits[j]);*/
338       }
339       /*printf ("\n");*/
340       if (psum > (total<<BITRES))
341          hi = mid;
342       else
343          lo = mid;
344       /*printf ("lo = %d, hi = %d\n", lo, hi);*/
345    }
346    /*printf ("interp between %d and %d\n", lo, hi);*/
347    for (j=start;j<end;j++)
348    {
349       int N = m->eBands[j+1]-m->eBands[j];
350       bits1[j] = (C*N*m->allocVectors[lo*len+j]<<LM>>2);
351       bits2[j] = (C*N*m->allocVectors[hi*len+j]<<LM>>2) - bits1[j];
352       if (bits1[j] > 0)
353          bits1[j] += trim_offset[j];
354       if (bits1[j] < 0)
355          bits1[j] = 0;
356       bits1[j] += offsets[j];
357    }
358    codedBands = interp_bits2pulses(m, start, end, bits1, bits2, thresh,
359          total, pulses, ebits, fine_priority, len, C, LM, skip, prev);
360    RESTORE_STACK;
361    return codedBands;
362 }
363