Give the bit we reserved to end skipping back when we don't use it.
[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, int skip_start,
144       const int *bits1, const int *bits2, const int *thresh, int total, int skip_rsv,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    int balance;
157    SAVE_STACK;
158
159    alloc_floor = C<<BITRES;
160
161    logM = LM<<BITRES;
162    lo = 0;
163    hi = 1<<ALLOC_STEPS;
164    for (i=0;i<ALLOC_STEPS;i++)
165    {
166       int mid = (lo+hi)>>1;
167       psum = 0;
168       done = 0;
169       for (j=end;j-->start;)
170       {
171          int tmp = bits1[j] + (mid*bits2[j]>>ALLOC_STEPS);
172          if (tmp >= thresh[j] || done)
173          {
174             done = 1;
175             /* Don't allocate more than we can actually use */
176             psum += IMIN(tmp, 64*C<<BITRES<<LM);
177          } else {
178             if (tmp >= alloc_floor)
179                psum += alloc_floor;
180          }
181       }
182       if (psum > total)
183          hi = mid;
184       else
185          lo = mid;
186    }
187    psum = 0;
188    /*printf ("interp bisection gave %d\n", lo);*/
189    done = 0;
190    for (j=end;j-->start;)
191    {
192       int tmp = bits1[j] + (lo*bits2[j]>>ALLOC_STEPS);
193       if (tmp < thresh[j] && !done)
194       {
195          if (tmp >= alloc_floor)
196             tmp = alloc_floor;
197          else
198             tmp = 0;
199       } else
200          done = 1;
201       /* Don't allocate more than we can actually use */
202       tmp = IMIN(tmp, 64*C<<BITRES<<LM);
203       bits[j] = tmp;
204       psum += tmp;
205    }
206
207    /* Decide which bands to skip, working backwards from the end. */
208    for (codedBands=end;;codedBands--)
209    {
210       int band_width;
211       int band_bits;
212       int rem;
213       j = codedBands-1;
214       /*Figure out how many left-over bits we would be adding to this band.
215         This can include bits we've stolen back from higher, skipped bands.*/
216       left = total-psum;
217       percoeff = left/(m->eBands[codedBands]-m->eBands[start]);
218       left -= (m->eBands[codedBands]-m->eBands[start])*percoeff;
219       /* Never skip the first band, nor a band that has been boosted by
220           dynalloc.
221          In the first case, we'd be coding a bit to signal we're going to waste
222           all the other bits.
223          In the second case, we'd be coding a bit to redistribute all the bits
224           we just signaled should be cocentrated in this band. */
225       if (j<=skip_start)
226       {
227          /* Give the bit we reserved to end skipping back to this band. */
228          bits[j] += skip_rsv;
229          break;
230       }
231       rem = IMAX(left-(m->eBands[j]-m->eBands[start]),0);
232       band_width = m->eBands[codedBands]-m->eBands[j];
233       band_bits = bits[j] + percoeff*band_width + rem;
234       /*Only code a skip decision if we're above the threshold for this band.
235         Otherwise it is force-skipped.
236         This ensures that we have enough bits to code the skip flag.*/
237       if (band_bits >= IMAX(thresh[j], alloc_floor+(1<<BITRES)))
238       {
239          if (encode)
240          {
241             /*This if() block is the only part of the allocation function that
242                is not a mandatory part of the bitstream: any bands we choose to
243                skip here must be explicitly signaled.*/
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             {
248                ec_enc_bit_prob((ec_enc *)ec, 1, 32768);
249                break;
250             }
251             ec_enc_bit_prob((ec_enc *)ec, 0, 32768);
252          } else if (ec_dec_bit_prob((ec_dec *)ec, 32768)) {
253             break;
254          }
255          /*We used a bit to skip this band.*/
256          psum += 1<<BITRES;
257          band_bits -= 1<<BITRES;
258       }
259       /*Reclaim the bits originally allocated to this band.*/
260       psum -= bits[j];
261       if (band_bits >= alloc_floor)
262       {
263          /*If we have enough for a fine energy bit per channel, use it.*/
264          psum += alloc_floor;
265          bits[j] = alloc_floor;
266       } else {
267          /*Otherwise this band gets nothing at all.*/
268          bits[j] = 0;
269       }
270    }
271
272    /* Allocate the remaining bits */
273    if (codedBands>start) {
274       for (j=start;j<codedBands;j++)
275          bits[j] += percoeff*(m->eBands[j+1]-m->eBands[j]);
276       for (j=start;j<codedBands;j++)
277       {
278          int tmp = IMIN(left, m->eBands[j+1]-m->eBands[j]);
279          bits[j] += tmp;
280          left -= tmp;
281       }
282    }
283    /*for (j=0;j<end;j++)printf("%d ", bits[j]);printf("\n");*/
284
285    balance = 0;
286    for (j=start;j<codedBands;j++)
287    {
288       int N0, N, den;
289       int offset;
290       int NClogN;
291
292       celt_assert(bits[j] >= 0);
293       N0 = m->eBands[j+1]-m->eBands[j];
294       N=N0<<LM;
295
296       if (N>1)
297       {
298          NClogN = N*C*(m->logN[j] + logM);
299
300          /* Compensate for the extra DoF in stereo */
301          den=(C*N+ ((C==2 && N>2) ? 1 : 0));
302
303          /* Offset for the number of fine bits by log2(N)/2 + FINE_OFFSET
304             compared to their "fair share" of total/N */
305          offset = (NClogN>>1)-N*C*FINE_OFFSET;
306
307          /* N=2 is the only point that doesn't match the curve */
308          if (N==2)
309             offset += N*C<<BITRES>>2;
310
311          /* Changing the offset for allocating the second and third
312              fine energy bit */
313          if (bits[j] + offset < den*2<<BITRES)
314             offset += NClogN>>2;
315          else if (bits[j] + offset < den*3<<BITRES)
316             offset += NClogN>>3;
317
318          /* Divide with rounding */
319          ebits[j] = IMAX(0, (bits[j] + offset + (den<<(BITRES-1))) / (den<<BITRES));
320
321          /* If we rounded down, make it a candidate for final
322              fine energy pass */
323          fine_priority[j] = ebits[j]*(den<<BITRES) >= bits[j]+offset;
324
325          /* Make sure not to bust */
326          if (C*ebits[j] > (bits[j]>>BITRES))
327             ebits[j] = bits[j]/C >> BITRES;
328
329          /* More than 7 is useless because that's about as far as PVQ can go */
330          if (ebits[j]>7)
331             ebits[j]=7;
332
333       } else {
334          /* For N=1, all bits go to fine energy except for a single sign bit */
335          ebits[j] = IMIN(IMAX(0,(bits[j]/C >> BITRES)-1),7);
336          fine_priority[j] = (ebits[j]+1)*C<<BITRES >= (bits[j]-balance);
337          /* N=1 bands can't take advantage of the re-balancing in
338              quant_all_bands() because they don't have shape, only fine energy.
339             Instead, do the re-balancing here.*/
340          balance = IMAX(0,bits[j] - ((ebits[j]+1)*C<<BITRES));
341          if (j+1<codedBands)
342          {
343             bits[j] -= balance;
344             bits[j+1] += balance;
345          }
346       }
347
348       /* The other bits are assigned to PVQ */
349       bits[j] -= C*ebits[j]<<BITRES;
350       celt_assert(bits[j] >= 0);
351       celt_assert(ebits[j] >= 0);
352    }
353    /* The skipped bands use all their bits for fine energy. */
354    for (;j<end;j++)
355    {
356       ebits[j] = bits[j]/C >> BITRES;
357       celt_assert(C*ebits[j]<<BITRES == bits[j]);
358       bits[j] = 0;
359       fine_priority[j] = ebits[j]<1;
360    }
361    RESTORE_STACK;
362    return codedBands;
363 }
364
365 int compute_allocation(const CELTMode *m, int start, int end, const int *offsets, int alloc_trim,
366       int total, int *pulses, int *ebits, int *fine_priority, int _C, int LM, void *ec, int encode, int prev)
367 {
368    int lo, hi, len, j;
369    const int C = CHANNELS(_C);
370    int codedBands;
371    int skip_start;
372    int skip_rsv;
373    VARDECL(int, bits1);
374    VARDECL(int, bits2);
375    VARDECL(int, thresh);
376    VARDECL(int, trim_offset);
377    SAVE_STACK;
378    
379    total = IMAX(total, 0);
380    len = m->nbEBands;
381    skip_start = start;
382    /* Reserve a bit to signal the end of manually skipped bands. */
383    skip_rsv = total >= 1<<BITRES ? 1<<BITRES : 0;
384    total -= skip_rsv;
385    ALLOC(bits1, len, int);
386    ALLOC(bits2, len, int);
387    ALLOC(thresh, len, int);
388    ALLOC(trim_offset, len, int);
389
390    /* Below this threshold, we're sure not to allocate any PVQ bits */
391    for (j=start;j<end;j++)
392       thresh[j] = IMAX((C)<<BITRES, (3*(m->eBands[j+1]-m->eBands[j])<<LM<<BITRES)>>4);
393    /* Tilt of the allocation curve */
394    for (j=start;j<end;j++)
395       trim_offset[j] = C*(m->eBands[j+1]-m->eBands[j])*(alloc_trim-5-LM)*(m->nbEBands-j-1)
396             <<(LM+BITRES)>>6;
397
398    lo = 0;
399    hi = m->nbAllocVectors - 1;
400    while (hi-lo != 1)
401    {
402       int psum = 0;
403       int mid = (lo+hi) >> 1;
404       for (j=start;j<end;j++)
405       {
406          int N = m->eBands[j+1]-m->eBands[j];
407          bits1[j] = C*N*m->allocVectors[mid*len+j]<<LM>>2;
408          if (bits1[j] > 0)
409             bits1[j] += trim_offset[j];
410          if (bits1[j] < 0)
411             bits1[j] = 0;
412          bits1[j] += offsets[j];
413          if (bits1[j] >= thresh[j])
414             psum += bits1[j];
415          else if (bits1[j] >= C<<BITRES)
416             psum += C<<BITRES;
417
418          /*printf ("%d ", bits[j]);*/
419       }
420       /*printf ("\n");*/
421       if (psum > total)
422          hi = mid;
423       else
424          lo = mid;
425       /*printf ("lo = %d, hi = %d\n", lo, hi);*/
426    }
427    /*printf ("interp between %d and %d\n", lo, hi);*/
428    for (j=start;j<end;j++)
429    {
430       int N = m->eBands[j+1]-m->eBands[j];
431       bits1[j] = (C*N*m->allocVectors[lo*len+j]<<LM>>2);
432       bits2[j] = (C*N*m->allocVectors[hi*len+j]<<LM>>2) - bits1[j];
433       if (bits1[j] > 0)
434          bits1[j] += trim_offset[j];
435       if (bits1[j] < 0)
436          bits1[j] = 0;
437       bits1[j] += offsets[j];
438       if (offsets[j]>0)
439          skip_start = j;
440    }
441    codedBands = interp_bits2pulses(m, start, end, skip_start, bits1, bits2, thresh,
442          total, skip_rsv, pulses, ebits, fine_priority, len, C, LM, ec, encode, prev);
443    RESTORE_STACK;
444    return codedBands;
445 }
446