Tim just rewrote half of the bit allocator -- hope it works now
[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    int alloc_floor;
154    int left, percoeff;
155    int force_skipping;
156    int unforced_skips;
157    int done;
158    SAVE_STACK;
159
160    alloc_floor = C<<BITRES;
161
162    logM = LM<<BITRES;
163    lo = 0;
164    hi = 1<<ALLOC_STEPS;
165    for (i=0;i<ALLOC_STEPS;i++)
166    {
167       int mid = (lo+hi)>>1;
168       psum = 0;
169       done = 0;
170       for (j=start;j<end;j++)
171       {
172          int tmp = bits1[j] + (mid*bits2[j]>>ALLOC_STEPS);
173          /* Don't allocate more than we can actually use */
174          if (tmp >= thresh[j] && !done)
175          {
176             psum += tmp;
177          } else {
178             done = 1;
179             if (tmp >= alloc_floor)
180                psum += alloc_floor;
181          }
182       }
183       if (psum > (total<<BITRES))
184          hi = mid;
185       else
186          lo = mid;
187    }
188    psum = 0;
189    /*printf ("interp bisection gave %d\n", lo);*/
190    done = 0;
191    for (j=start;j<end;j++)
192    {
193       int tmp = bits1[j] + (lo*bits2[j]>>ALLOC_STEPS);
194       if (tmp < thresh[j] || done)
195       {
196          done = 1;
197          if (tmp >= alloc_floor)
198             tmp = alloc_floor;
199          else
200             tmp = 0;
201       }
202       /* Don't allocate more than we can actually use */
203       tmp = IMIN(tmp, 64*C<<BITRES<<LM);
204       bits[j] = tmp;
205       psum += tmp;
206    }
207
208    force_skipping = 1;
209    unforced_skips = *skip;
210    *skip = 0;
211    codedBands=end;
212    for (j=end;j-->start;)
213    {
214       int band_width;
215       int band_bits;
216       int rem;
217       /*Figure out how many left-over bits we would be adding to this band.
218         This can include bits we've stolen back from higher, skipped bands.*/
219       left = (total<<BITRES)-psum;
220       percoeff = left/(m->eBands[codedBands]-m->eBands[start]);
221       left -= (m->eBands[codedBands]-m->eBands[start])*percoeff;
222       rem = IMAX(left-m->eBands[j],0);
223       band_width = m->eBands[codedBands]-m->eBands[j];
224       band_bits = bits[j] + percoeff*band_width + rem;
225       /*As long as, even after adding these bits, we're below the threshold for
226          this band, it is force-skipped.*/
227       force_skipping = force_skipping && band_bits < thresh[j];
228       if (!force_skipping)
229       {
230          /*If we have enough for the fine energy, but not more than a full bit
231             beyond that, or no more than one bit total, then don't bother
232             skipping this band: there's no extra bits to redistribute.*/
233          if ((band_bits >= alloc_floor && band_bits <= alloc_floor + (1<<BITRES))
234                || band_bits < (1<<BITRES))
235             break;
236          /*Never skip the first band: we'd be coding a bit to signal that we're
237             going to waste all of the other bits.*/
238          if (j==start)break;
239          if (unforced_skips == -1)
240          {
241             /*Choose a threshold with some hysteresis to keep bands from
242                fluctuating in and out.*/
243             if (band_bits > ((j<prev?7:9)*band_width<<LM<<BITRES)>>4)
244                break;
245          } else if(unforced_skips--<=0)
246             break;
247          (*skip)++;
248          /*Use a bit to skip this band.*/
249          psum += 1<<BITRES;
250          band_bits -= 1<<BITRES;
251       }
252       /*Reclaim the bits originally allocated to this band.*/
253       psum -= bits[j];
254       if (band_bits >= alloc_floor + (1<<BITRES))
255       {
256          /*If we have enough for a fine energy bit per channel, use it.*/
257          psum += alloc_floor;
258          bits[codedBands-1] = alloc_floor;
259       } else {
260          /*Otherwise this band gets nothing at all.*/
261          bits[codedBands-1] = 0;
262       }
263       codedBands--;
264    }
265
266    /* Allocate the remaining bits */
267    if (codedBands>start) {
268       for (j=start;j<codedBands;j++)
269          bits[j] += percoeff*(m->eBands[j+1]-m->eBands[j]);
270       for (j=start;j<codedBands;j++)
271       {
272          int tmp = IMIN(left, m->eBands[j+1]-m->eBands[j]);
273          bits[j] += tmp;
274          left -= tmp;
275       }
276    }
277    /*for (j=0;j<end;j++)printf("%d ", bits[j]);printf("\n");*/
278    for (j=start;j<end;j++)
279    {
280       int N0, N, den;
281       int offset;
282       int NClogN;
283
284       celt_assert(bits[j] >= 0);
285       N0 = m->eBands[j+1]-m->eBands[j];
286       N=N0<<LM;
287       NClogN = N*C*(m->logN[j] + logM);
288
289       /* Compensate for the extra DoF in stereo */
290       den=(C*N+ ((C==2 && N>2) ? 1 : 0));
291
292       /* Offset for the number of fine bits by log2(N)/2 + FINE_OFFSET
293          compared to their "fair share" of total/N */
294       offset = (NClogN>>1)-N*C*FINE_OFFSET;
295
296       /* N=2 is the only point that doesn't match the curve */
297       if (N==2)
298          offset += N*C<<BITRES>>2;
299
300       /* Changing the offset for allocating the second and third fine energy bit */
301       if (bits[j] + offset < den*2<<BITRES)
302          offset += NClogN>>2;
303       else if (bits[j] + offset < den*3<<BITRES)
304          offset += NClogN>>3;
305
306       /* Divide with rounding */
307       ebits[j] = IMAX(0, (bits[j] + offset + (den<<(BITRES-1))) / (den<<BITRES));
308
309       /* If we rounded down, make it a candidate for final fine energy pass */
310       fine_priority[j] = ebits[j]*(den<<BITRES) >= bits[j]+offset;
311
312       /* For N=1, all bits go to fine energy except for a single sign bit */
313       if (N==1)
314       {
315          ebits[j] = IMAX(0,(bits[j]/C >> BITRES)-1);
316          fine_priority[j] = (ebits[j]+1)*C<<BITRES >= bits[j];
317       }
318       /* Make sure not to bust */
319       if (C*ebits[j] > (bits[j]>>BITRES))
320          ebits[j] = bits[j]/C >> BITRES;
321
322       /* More than that is useless because that's about as far as PVQ can go */
323       if (ebits[j]>7)
324          ebits[j]=7;
325
326       /* The other bits are assigned to PVQ */
327       bits[j] -= C*ebits[j]<<BITRES;
328       celt_assert(bits[j] >= 0);
329       celt_assert(ebits[j] >= 0);
330    }
331    RESTORE_STACK;
332    return codedBands;
333 }
334
335 int compute_allocation(const CELTMode *m, int start, int end, int *offsets, int alloc_trim,
336       int total, int *pulses, int *ebits, int *fine_priority, int _C, int LM, int *skip, int prev)
337 {
338    int lo, hi, len, j;
339    const int C = CHANNELS(_C);
340    int codedBands;
341    VARDECL(int, bits1);
342    VARDECL(int, bits2);
343    VARDECL(int, thresh);
344    VARDECL(int, trim_offset);
345    SAVE_STACK;
346    
347    total = IMAX(total, 0);
348    len = m->nbEBands;
349    ALLOC(bits1, len, int);
350    ALLOC(bits2, len, int);
351    ALLOC(thresh, len, int);
352    ALLOC(trim_offset, len, int);
353
354    /* Below this threshold, we're sure not to allocate any PVQ bits */
355    for (j=start;j<end;j++)
356       thresh[j] = IMAX((C)<<BITRES, (3*(m->eBands[j+1]-m->eBands[j])<<LM<<BITRES)>>4);
357    /* Tilt of the allocation curve */
358    for (j=start;j<end;j++)
359       trim_offset[j] = C*(m->eBands[j+1]-m->eBands[j])*(alloc_trim-5-LM)*(m->nbEBands-j-1)
360             <<(LM+BITRES)>>6;
361
362    lo = 0;
363    hi = m->nbAllocVectors - 1;
364    while (hi-lo != 1)
365    {
366       int psum = 0;
367       int mid = (lo+hi) >> 1;
368       for (j=start;j<end;j++)
369       {
370          int N = m->eBands[j+1]-m->eBands[j];
371          bits1[j] = C*N*m->allocVectors[mid*len+j]<<LM>>2;
372          if (bits1[j] > 0)
373             bits1[j] += trim_offset[j];
374          if (bits1[j] < 0)
375             bits1[j] = 0;
376          bits1[j] += offsets[j];
377          if (bits1[j] >= thresh[j])
378             psum += bits1[j];
379          else if (bits1[j] >= C<<BITRES)
380             psum += C<<BITRES;
381
382          /*printf ("%d ", bits[j]);*/
383       }
384       /*printf ("\n");*/
385       if (psum > (total<<BITRES))
386          hi = mid;
387       else
388          lo = mid;
389       /*printf ("lo = %d, hi = %d\n", lo, hi);*/
390    }
391    /*printf ("interp between %d and %d\n", lo, hi);*/
392    for (j=start;j<end;j++)
393    {
394       int N = m->eBands[j+1]-m->eBands[j];
395       bits1[j] = (C*N*m->allocVectors[lo*len+j]<<LM>>2);
396       bits2[j] = (C*N*m->allocVectors[hi*len+j]<<LM>>2) - bits1[j];
397       if (bits1[j] > 0)
398          bits1[j] += trim_offset[j];
399       if (bits1[j] < 0)
400          bits1[j] = 0;
401       bits1[j] += offsets[j];
402    }
403    codedBands = interp_bits2pulses(m, start, end, bits1, bits2, thresh,
404          total, pulses, ebits, fine_priority, len, C, LM, skip, prev);
405    RESTORE_STACK;
406    return codedBands;
407 }
408