More allocation function updates.
[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 (unforced_skips == -1)
231          {
232             /*This if() block is the only part of the allocation function that
233                is not a mandatory part of the bitstream: any bands we choose to
234                skip here must be explicitly signaled.*/
235             /*If we have enough for the fine energy, but not more than a full
236                bit beyond that, or no more than one bit total, then don't bother
237                skipping this band: there's no extra bits to redistribute.*/
238             if ((band_bits>=alloc_floor && band_bits<=alloc_floor+(1<<BITRES))
239                   || band_bits<(1<<BITRES))
240                break;
241             /*Never skip the first band: we'd be coding a bit to signal that
242                we're going to waste all of the other bits.*/
243             if (j==start)break;
244             /*Choose a threshold with some hysteresis to keep bands from
245                fluctuating in and out.*/
246             if (band_bits > ((j<prev?7:9)*band_width<<LM<<BITRES)>>4)
247                break;
248          } else if(unforced_skips--<=0)
249             break;
250          (*skip)++;
251          /*Use a bit to skip this band.*/
252          psum += 1<<BITRES;
253          band_bits -= 1<<BITRES;
254       }
255       /*Reclaim the bits originally allocated to this band.*/
256       psum -= bits[j];
257       if (band_bits >= alloc_floor)
258       {
259          /*If we have enough for a fine energy bit per channel, use it.*/
260          psum += alloc_floor;
261          bits[j] = alloc_floor;
262       } else {
263          /*Otherwise this band gets nothing at all.*/
264          bits[j] = 0;
265       }
266       codedBands--;
267    }
268
269    /* Allocate the remaining bits */
270    if (codedBands>start) {
271       for (j=start;j<codedBands;j++)
272          bits[j] += percoeff*(m->eBands[j+1]-m->eBands[j]);
273       for (j=start;j<codedBands;j++)
274       {
275          int tmp = IMIN(left, m->eBands[j+1]-m->eBands[j]);
276          bits[j] += tmp;
277          left -= tmp;
278       }
279    }
280    /*for (j=0;j<end;j++)printf("%d ", bits[j]);printf("\n");*/
281    for (j=start;j<end;j++)
282    {
283       int N0, N, den;
284       int offset;
285       int NClogN;
286
287       celt_assert(bits[j] >= 0);
288       N0 = m->eBands[j+1]-m->eBands[j];
289       N=N0<<LM;
290       NClogN = N*C*(m->logN[j] + logM);
291
292       /* Compensate for the extra DoF in stereo */
293       den=(C*N+ ((C==2 && N>2) ? 1 : 0));
294
295       /* Offset for the number of fine bits by log2(N)/2 + FINE_OFFSET
296          compared to their "fair share" of total/N */
297       offset = (NClogN>>1)-N*C*FINE_OFFSET;
298
299       /* N=2 is the only point that doesn't match the curve */
300       if (N==2)
301          offset += N*C<<BITRES>>2;
302
303       /* Changing the offset for allocating the second and third fine energy bit */
304       if (bits[j] + offset < den*2<<BITRES)
305          offset += NClogN>>2;
306       else if (bits[j] + offset < den*3<<BITRES)
307          offset += NClogN>>3;
308
309       /* Divide with rounding */
310       ebits[j] = IMAX(0, (bits[j] + offset + (den<<(BITRES-1))) / (den<<BITRES));
311
312       /* If we rounded down, make it a candidate for final fine energy pass */
313       fine_priority[j] = ebits[j]*(den<<BITRES) >= bits[j]+offset;
314
315       /* For N=1, all bits go to fine energy except for a single sign bit */
316       if (N==1)
317       {
318          ebits[j] = IMAX(0,(bits[j]/C >> BITRES)-1);
319          fine_priority[j] = (ebits[j]+1)*C<<BITRES >= bits[j];
320       }
321       /* Make sure not to bust */
322       if (C*ebits[j] > (bits[j]>>BITRES))
323          ebits[j] = bits[j]/C >> BITRES;
324
325       /* More than that is useless because that's about as far as PVQ can go */
326       if (ebits[j]>7)
327          ebits[j]=7;
328
329       /* The other bits are assigned to PVQ */
330       bits[j] -= C*ebits[j]<<BITRES;
331       celt_assert(bits[j] >= 0);
332       celt_assert(ebits[j] >= 0);
333    }
334    RESTORE_STACK;
335    return codedBands;
336 }
337
338 int compute_allocation(const CELTMode *m, int start, int end, int *offsets, int alloc_trim,
339       int total, int *pulses, int *ebits, int *fine_priority, int _C, int LM, int *skip, int prev)
340 {
341    int lo, hi, len, j;
342    const int C = CHANNELS(_C);
343    int codedBands;
344    VARDECL(int, bits1);
345    VARDECL(int, bits2);
346    VARDECL(int, thresh);
347    VARDECL(int, trim_offset);
348    SAVE_STACK;
349    
350    total = IMAX(total, 0);
351    len = m->nbEBands;
352    ALLOC(bits1, len, int);
353    ALLOC(bits2, len, int);
354    ALLOC(thresh, len, int);
355    ALLOC(trim_offset, len, int);
356
357    /* Below this threshold, we're sure not to allocate any PVQ bits */
358    for (j=start;j<end;j++)
359       thresh[j] = IMAX((C)<<BITRES, (3*(m->eBands[j+1]-m->eBands[j])<<LM<<BITRES)>>4);
360    /* Tilt of the allocation curve */
361    for (j=start;j<end;j++)
362       trim_offset[j] = C*(m->eBands[j+1]-m->eBands[j])*(alloc_trim-5-LM)*(m->nbEBands-j-1)
363             <<(LM+BITRES)>>6;
364
365    lo = 0;
366    hi = m->nbAllocVectors - 1;
367    while (hi-lo != 1)
368    {
369       int psum = 0;
370       int mid = (lo+hi) >> 1;
371       for (j=start;j<end;j++)
372       {
373          int N = m->eBands[j+1]-m->eBands[j];
374          bits1[j] = C*N*m->allocVectors[mid*len+j]<<LM>>2;
375          if (bits1[j] > 0)
376             bits1[j] += trim_offset[j];
377          if (bits1[j] < 0)
378             bits1[j] = 0;
379          bits1[j] += offsets[j];
380          if (bits1[j] >= thresh[j])
381             psum += bits1[j];
382          else if (bits1[j] >= C<<BITRES)
383             psum += C<<BITRES;
384
385          /*printf ("%d ", bits[j]);*/
386       }
387       /*printf ("\n");*/
388       if (psum > (total<<BITRES))
389          hi = mid;
390       else
391          lo = mid;
392       /*printf ("lo = %d, hi = %d\n", lo, hi);*/
393    }
394    /*printf ("interp between %d and %d\n", lo, hi);*/
395    for (j=start;j<end;j++)
396    {
397       int N = m->eBands[j+1]-m->eBands[j];
398       bits1[j] = (C*N*m->allocVectors[lo*len+j]<<LM>>2);
399       bits2[j] = (C*N*m->allocVectors[hi*len+j]<<LM>>2) - bits1[j];
400       if (bits1[j] > 0)
401          bits1[j] += trim_offset[j];
402       if (bits1[j] < 0)
403          bits1[j] = 0;
404       bits1[j] += offsets[j];
405    }
406    codedBands = interp_bits2pulses(m, start, end, bits1, bits2, thresh,
407          total, pulses, ebits, fine_priority, len, C, LM, skip, prev);
408    RESTORE_STACK;
409    return codedBands;
410 }
411