Change strategies for allocation hole prevention.
[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, void *ec, int encode, 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 done;
156    SAVE_STACK;
157
158    alloc_floor = C<<BITRES;
159
160    logM = LM<<BITRES;
161    lo = 0;
162    hi = 1<<ALLOC_STEPS;
163    for (i=0;i<ALLOC_STEPS;i++)
164    {
165       int mid = (lo+hi)>>1;
166       psum = 0;
167       done = 0;
168       for (j=end;j-->start;)
169       {
170          int tmp = bits1[j] + (mid*bits2[j]>>ALLOC_STEPS);
171          if (tmp >= thresh[j] || done)
172          {
173             done = 1;
174             /* Don't allocate more than we can actually use */
175             psum += IMIN(tmp, 64*C<<BITRES<<LM);
176          } else {
177             if (tmp >= alloc_floor)
178                psum += alloc_floor;
179          }
180       }
181       if (psum > total)
182          hi = mid;
183       else
184          lo = mid;
185    }
186    psum = 0;
187    /*printf ("interp bisection gave %d\n", lo);*/
188    done = 0;
189    for (j=end;j-->start;)
190    {
191       int tmp = bits1[j] + (lo*bits2[j]>>ALLOC_STEPS);
192       if (tmp < thresh[j] && !done)
193       {
194          if (tmp >= alloc_floor)
195             tmp = alloc_floor;
196          else
197             tmp = 0;
198       } else
199          done = 1;
200       /* Don't allocate more than we can actually use */
201       tmp = IMIN(tmp, 64*C<<BITRES<<LM);
202       bits[j] = tmp;
203       psum += tmp;
204    }
205
206    /* Decide which bands to skip, working backwards from the end. */
207    for (codedBands=end;;codedBands--)
208    {
209       int band_width;
210       int band_bits;
211       int rem;
212       j = codedBands-1;
213       /*Figure out how many left-over bits we would be adding to this band.
214         This can include bits we've stolen back from higher, skipped bands.*/
215       left = total-psum;
216       percoeff = left/(m->eBands[codedBands]-m->eBands[start]);
217       left -= (m->eBands[codedBands]-m->eBands[start])*percoeff;
218       /* Never skip the first band: we'd be coding a bit to signal that we're
219           going to waste all the other bits.
220          This means we won't be using the extra bit we reserved to signal the
221           end of manual skipping, but that will get added back in by
222           quant_all_bands().*/
223       if (j<=start)
224          break;
225       rem = IMAX(left-(m->eBands[j]-m->eBands[start]),0);
226       band_width = m->eBands[codedBands]-m->eBands[j];
227       band_bits = bits[j] + percoeff*band_width + rem;
228       /*Only code a skip decision if we're above the threshold for this band.
229         Otherwise it is force-skipped.
230         This ensures that we have enough bits to code the skip flag.*/
231       if (band_bits >= IMAX(thresh[j], alloc_floor+(1<<BITRES)))
232       {
233          if (encode)
234          {
235             /*This if() block is the only part of the allocation function that
236                is not a mandatory part of the bitstream: any bands we choose to
237                skip here must be explicitly signaled.*/
238             /*Choose a threshold with some hysteresis to keep bands from
239                fluctuating in and out.*/
240             if (band_bits > ((j<prev?7:9)*band_width<<LM<<BITRES)>>4)
241             {
242                ec_enc_bit_prob((ec_enc *)ec, 1, 32768);
243                break;
244             }
245             ec_enc_bit_prob((ec_enc *)ec, 0, 32768);
246          } else if (ec_dec_bit_prob((ec_dec *)ec, 32768)) {
247             break;
248          }
249          /*We used a bit to skip this band.*/
250          psum += 1<<BITRES;
251          band_bits -= 1<<BITRES;
252       }
253       /*Reclaim the bits originally allocated to this band.*/
254       psum -= bits[j];
255       if (band_bits >= alloc_floor)
256       {
257          /*If we have enough for a fine energy bit per channel, use it.*/
258          psum += alloc_floor;
259          bits[j] = alloc_floor;
260       } else {
261          /*Otherwise this band gets nothing at all.*/
262          bits[j] = 0;
263       }
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<codedBands;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    /* The skipped bands use all their bits for fine energy. */
332    for (;j<end;j++)
333    {
334       ebits[j] = bits[j]/C >> BITRES;
335       celt_assert(C*ebits[j]<<BITRES == bits[j]);
336       bits[j] = 0;
337       fine_priority[j] = 0;
338    }
339    RESTORE_STACK;
340    return codedBands;
341 }
342
343 int compute_allocation(const CELTMode *m, int start, int end, int *offsets, int alloc_trim,
344       int total, int *pulses, int *ebits, int *fine_priority, int _C, int LM, void *ec, int encode, int prev)
345 {
346    int lo, hi, len, j;
347    const int C = CHANNELS(_C);
348    int codedBands;
349    VARDECL(int, bits1);
350    VARDECL(int, bits2);
351    VARDECL(int, thresh);
352    VARDECL(int, trim_offset);
353    SAVE_STACK;
354    
355    total = IMAX(total, 0);
356    len = m->nbEBands;
357    ALLOC(bits1, len, int);
358    ALLOC(bits2, len, int);
359    ALLOC(thresh, len, int);
360    ALLOC(trim_offset, len, int);
361
362    /* Below this threshold, we're sure not to allocate any PVQ bits */
363    for (j=start;j<end;j++)
364       thresh[j] = IMAX((C)<<BITRES, (3*(m->eBands[j+1]-m->eBands[j])<<LM<<BITRES)>>4);
365    /* Tilt of the allocation curve */
366    for (j=start;j<end;j++)
367       trim_offset[j] = C*(m->eBands[j+1]-m->eBands[j])*(alloc_trim-5-LM)*(m->nbEBands-j-1)
368             <<(LM+BITRES)>>6;
369
370    lo = 0;
371    hi = m->nbAllocVectors - 1;
372    while (hi-lo != 1)
373    {
374       int psum = 0;
375       int mid = (lo+hi) >> 1;
376       for (j=start;j<end;j++)
377       {
378          int N = m->eBands[j+1]-m->eBands[j];
379          bits1[j] = C*N*m->allocVectors[mid*len+j]<<LM>>2;
380          if (bits1[j] > 0)
381             bits1[j] += trim_offset[j];
382          if (bits1[j] < 0)
383             bits1[j] = 0;
384          bits1[j] += offsets[j];
385          if (bits1[j] >= thresh[j])
386             psum += bits1[j];
387          else if (bits1[j] >= C<<BITRES)
388             psum += C<<BITRES;
389
390          /*printf ("%d ", bits[j]);*/
391       }
392       /*printf ("\n");*/
393       if (psum > total)
394          hi = mid;
395       else
396          lo = mid;
397       /*printf ("lo = %d, hi = %d\n", lo, hi);*/
398    }
399    /*printf ("interp between %d and %d\n", lo, hi);*/
400    for (j=start;j<end;j++)
401    {
402       int N = m->eBands[j+1]-m->eBands[j];
403       bits1[j] = (C*N*m->allocVectors[lo*len+j]<<LM>>2);
404       bits2[j] = (C*N*m->allocVectors[hi*len+j]<<LM>>2) - bits1[j];
405       if (bits1[j] > 0)
406          bits1[j] += trim_offset[j];
407       if (bits1[j] < 0)
408          bits1[j] = 0;
409       bits1[j] += offsets[j];
410    }
411    codedBands = interp_bits2pulses(m, start, end, bits1, bits2, thresh,
412          total, pulses, ebits, fine_priority, len, C, LM, ec, encode, prev);
413    RESTORE_STACK;
414    return codedBands;
415 }
416