Don't rebalance bits for itheta=0 or 16384
[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 static const unsigned char LOG2_FRAC_TABLE[24]={
48    0,
49    8,13,
50   16,19,21,23,
51   24,26,27,28,29,30,31,32,
52   32,33,34,34,35,36,36,37,37
53 };
54
55 #ifndef STATIC_MODES
56
57 /*Determines if V(N,K) fits in a 32-bit unsigned integer.
58   N and K are themselves limited to 15 bits.*/
59 static int fits_in32(int _n, int _k)
60 {
61    static const celt_int16 maxN[15] = {
62       32767, 32767, 32767, 1476, 283, 109,  60,  40,
63        29,  24,  20,  18,  16,  14,  13};
64    static const celt_int16 maxK[15] = {
65       32767, 32767, 32767, 32767, 1172, 238,  95,  53,
66        36,  27,  22,  18,  16,  15,  13};
67    if (_n>=14)
68    {
69       if (_k>=14)
70          return 0;
71       else
72          return _n <= maxN[_k];
73    } else {
74       return _k <= maxK[_n];
75    }
76 }
77
78 void compute_pulse_cache(CELTMode *m, int LM)
79 {
80    int i;
81    int curr=0;
82    int nbEntries=0;
83    int entryN[100], entryK[100], entryI[100];
84    const celt_int16 *eBands = m->eBands;
85    PulseCache *cache = &m->cache;
86    celt_int16 *cindex;
87    unsigned char *bits;
88
89    cindex = celt_alloc(sizeof(cache->index[0])*m->nbEBands*(LM+2));
90    cache->index = cindex;
91
92    /* Scan for all unique band sizes */
93    for (i=0;i<=LM+1;i++)
94    {
95       int j;
96       for (j=0;j<m->nbEBands;j++)
97       {
98          int k;
99          int N = (eBands[j+1]-eBands[j])<<i>>1;
100          cindex[i*m->nbEBands+j] = -1;
101          /* Find other bands that have the same size */
102          for (k=0;k<=i;k++)
103          {
104             int n;
105             for (n=0;n<m->nbEBands && (k!=i || n<j);n++)
106             {
107                if (N == (eBands[n+1]-eBands[n])<<k>>1)
108                {
109                   cindex[i*m->nbEBands+j] = cindex[k*m->nbEBands+n];
110                   break;
111                }
112             }
113          }
114          if (cache->index[i*m->nbEBands+j] == -1 && N!=0)
115          {
116             int K;
117             entryN[nbEntries] = N;
118             K = 0;
119             while (fits_in32(N,get_pulses(K+1)) && K<MAX_PSEUDO)
120                K++;
121             entryK[nbEntries] = K;
122             cindex[i*m->nbEBands+j] = curr;
123             entryI[nbEntries] = curr;
124
125             curr += K+1;
126             nbEntries++;
127          }
128       }
129    }
130    bits = celt_alloc(sizeof(unsigned char)*curr);
131    cache->bits = bits;
132    cache->size = curr;
133    /* Compute the cache for all unique sizes */
134    for (i=0;i<nbEntries;i++)
135    {
136       int j;
137       unsigned char *ptr = bits+entryI[i];
138       celt_int16 tmp[MAX_PULSES+1];
139       get_required_bits(tmp, entryN[i], get_pulses(entryK[i]), BITRES);
140       for (j=1;j<=entryK[i];j++)
141          ptr[j] = tmp[get_pulses(j)]-1;
142       ptr[0] = entryK[i];
143    }
144 }
145
146 #endif /* !STATIC_MODES */
147
148
149 #define ALLOC_STEPS 6
150
151 static inline int interp_bits2pulses(const CELTMode *m, int start, int end, int skip_start,
152       const int *bits1, const int *bits2, const int *thresh, int total, int skip_rsv,
153       int *intensity, int intensity_rsv, int *dual_stereo, int dual_stereo_rsv, int *bits,
154       int *ebits, int *fine_priority, int _C, int LM, void *ec, int encode, int prev)
155 {
156    int psum;
157    int lo, hi;
158    int i, j;
159    int logM;
160    const int C = CHANNELS(_C);
161    int stereo;
162    int codedBands=-1;
163    int alloc_floor;
164    int left, percoeff;
165    int done;
166    int balance;
167    SAVE_STACK;
168
169    alloc_floor = C<<BITRES;
170    stereo = C>1;
171
172    logM = LM<<BITRES;
173    lo = 0;
174    hi = 1<<ALLOC_STEPS;
175    for (i=0;i<ALLOC_STEPS;i++)
176    {
177       int mid = (lo+hi)>>1;
178       psum = 0;
179       done = 0;
180       for (j=end;j-->start;)
181       {
182          int tmp = bits1[j] + (mid*bits2[j]>>ALLOC_STEPS);
183          if (tmp >= thresh[j] || done)
184          {
185             done = 1;
186             /* Don't allocate more than we can actually use */
187             psum += IMIN(tmp, 64*C<<BITRES<<LM);
188          } else {
189             if (tmp >= alloc_floor)
190                psum += alloc_floor;
191          }
192       }
193       if (psum > total)
194          hi = mid;
195       else
196          lo = mid;
197    }
198    psum = 0;
199    /*printf ("interp bisection gave %d\n", lo);*/
200    done = 0;
201    for (j=end;j-->start;)
202    {
203       int tmp = bits1[j] + (lo*bits2[j]>>ALLOC_STEPS);
204       if (tmp < thresh[j] && !done)
205       {
206          if (tmp >= alloc_floor)
207             tmp = alloc_floor;
208          else
209             tmp = 0;
210       } else
211          done = 1;
212       /* Don't allocate more than we can actually use */
213       tmp = IMIN(tmp, 64*C<<BITRES<<LM);
214       bits[j] = tmp;
215       psum += tmp;
216    }
217
218    /* Decide which bands to skip, working backwards from the end. */
219    for (codedBands=end;;codedBands--)
220    {
221       int band_width;
222       int band_bits;
223       int rem;
224       j = codedBands-1;
225       /* Never skip the first band, nor a band that has been boosted by
226           dynalloc.
227          In the first case, we'd be coding a bit to signal we're going to waste
228           all the other bits.
229          In the second case, we'd be coding a bit to redistribute all the bits
230           we just signaled should be cocentrated in this band. */
231       if (j<=skip_start)
232       {
233          /* Give the bit we reserved to end skipping back. */
234          total += skip_rsv;
235          break;
236       }
237       /*Figure out how many left-over bits we would be adding to this band.
238         This can include bits we've stolen back from higher, skipped bands.*/
239       left = total-psum;
240       percoeff = left/(m->eBands[codedBands]-m->eBands[start]);
241       left -= (m->eBands[codedBands]-m->eBands[start])*percoeff;
242       rem = IMAX(left-(m->eBands[j]-m->eBands[start]),0);
243       band_width = m->eBands[codedBands]-m->eBands[j];
244       band_bits = bits[j] + percoeff*band_width + rem;
245       /*Only code a skip decision if we're above the threshold for this band.
246         Otherwise it is force-skipped.
247         This ensures that we have enough bits to code the skip flag.*/
248       if (band_bits >= IMAX(thresh[j], alloc_floor+(1<<BITRES)))
249       {
250          if (encode)
251          {
252             /*This if() block is the only part of the allocation function that
253                is not a mandatory part of the bitstream: any bands we choose to
254                skip here must be explicitly signaled.*/
255             /*Choose a threshold with some hysteresis to keep bands from
256                fluctuating in and out.*/
257             if (band_bits > ((j<prev?7:9)*band_width<<LM<<BITRES)>>4)
258             {
259                ec_enc_bit_logp((ec_enc *)ec, 1, 1);
260                break;
261             }
262             ec_enc_bit_logp((ec_enc *)ec, 0, 1);
263          } else if (ec_dec_bit_logp((ec_dec *)ec, 1)) {
264             break;
265          }
266          /*We used a bit to skip this band.*/
267          psum += 1<<BITRES;
268          band_bits -= 1<<BITRES;
269       }
270       /*Reclaim the bits originally allocated to this band.*/
271       psum -= bits[j]+intensity_rsv;
272       if (intensity_rsv > 0)
273          intensity_rsv = LOG2_FRAC_TABLE[j-start];
274       psum += intensity_rsv;
275       if (band_bits >= alloc_floor)
276       {
277          /*If we have enough for a fine energy bit per channel, use it.*/
278          psum += alloc_floor;
279          bits[j] = alloc_floor;
280       } else {
281          /*Otherwise this band gets nothing at all.*/
282          bits[j] = 0;
283       }
284    }
285
286    celt_assert(codedBands > start);
287    /* Code the intensity and dual stereo parameters. */
288    if (intensity_rsv > 0)
289    {
290       if (encode)
291       {
292          *intensity = IMIN(*intensity, codedBands);
293          ec_enc_uint((ec_enc *)ec, *intensity-start, codedBands+1-start);
294       }
295       else
296          *intensity = start+ec_dec_uint((ec_dec *)ec, codedBands+1-start);
297    }
298    else
299       *intensity = 0;
300    if (*intensity <= start)
301    {
302       total += dual_stereo_rsv;
303       dual_stereo_rsv = 0;
304    }
305    if (dual_stereo_rsv > 0)
306    {
307       if (encode)
308          ec_enc_bit_logp((ec_enc *)ec, *dual_stereo, 1);
309       else
310          *dual_stereo = ec_dec_bit_logp((ec_dec *)ec, 1);
311    }
312    else
313       *dual_stereo = 0;
314
315    /* Allocate the remaining bits */
316    left = total-psum;
317    percoeff = left/(m->eBands[codedBands]-m->eBands[start]);
318    left -= (m->eBands[codedBands]-m->eBands[start])*percoeff;
319    for (j=start;j<codedBands;j++)
320       bits[j] += percoeff*(m->eBands[j+1]-m->eBands[j]);
321    for (j=start;j<codedBands;j++)
322    {
323       int tmp = IMIN(left, m->eBands[j+1]-m->eBands[j]);
324       bits[j] += tmp;
325       left -= tmp;
326    }
327    /*for (j=0;j<end;j++)printf("%d ", bits[j]);printf("\n");*/
328
329    balance = 0;
330    for (j=start;j<codedBands;j++)
331    {
332       int N0, N, den;
333       int offset;
334       int NClogN;
335
336       celt_assert(bits[j] >= 0);
337       N0 = m->eBands[j+1]-m->eBands[j];
338       N=N0<<LM;
339
340       if (N>1)
341       {
342          NClogN = N*C*(m->logN[j] + logM);
343
344          /* Compensate for the extra DoF in stereo */
345          den=(C*N+ ((C==2 && N>2) ? 1 : 0));
346
347          /* Offset for the number of fine bits by log2(N)/2 + FINE_OFFSET
348             compared to their "fair share" of total/N */
349          offset = (NClogN>>1)-N*C*FINE_OFFSET;
350
351          /* N=2 is the only point that doesn't match the curve */
352          if (N==2)
353             offset += N*C<<BITRES>>2;
354
355          /* Changing the offset for allocating the second and third
356              fine energy bit */
357          if (bits[j] + offset < den*2<<BITRES)
358             offset += NClogN>>2;
359          else if (bits[j] + offset < den*3<<BITRES)
360             offset += NClogN>>3;
361
362          /* Divide with rounding */
363          ebits[j] = IMAX(0, (bits[j] + offset + (den<<(BITRES-1))) / (den<<BITRES));
364
365          /* Make sure not to bust */
366          if (C*ebits[j] > (bits[j]>>BITRES))
367             ebits[j] = bits[j] >> stereo >> BITRES;
368
369          /* More than that is useless because that's about as far as PVQ can go */
370          ebits[j] = IMIN(ebits[j], MAX_FINE_BITS);
371
372          /* If we rounded down or capped this band, make it a candidate for the
373              final fine energy pass */
374          fine_priority[j] = ebits[j]*(den<<BITRES) >= bits[j]+offset;
375
376       } else {
377          /* For N=1, all bits go to fine energy except for a single sign bit */
378          ebits[j] = IMIN(IMAX(0,(bits[j] >> stereo >> BITRES)-1),MAX_FINE_BITS);
379          fine_priority[j] = (ebits[j]+1)*C<<BITRES >= (bits[j]-balance);
380          /* N=1 bands can't take advantage of the re-balancing in
381              quant_all_bands() because they don't have shape, only fine energy.
382             Instead, do the re-balancing here.*/
383          balance = IMAX(0,bits[j] - ((ebits[j]+1)*C<<BITRES));
384          if (j+1<codedBands)
385          {
386             bits[j] -= balance;
387             bits[j+1] += balance;
388          }
389       }
390
391       /* The other bits are assigned to PVQ */
392       bits[j] -= C*ebits[j]<<BITRES;
393       celt_assert(bits[j] >= 0);
394       celt_assert(ebits[j] >= 0);
395    }
396    /* The skipped bands use all their bits for fine energy. */
397    for (;j<end;j++)
398    {
399       ebits[j] = bits[j] >> stereo >> BITRES;
400       celt_assert(C*ebits[j]<<BITRES == bits[j]);
401       bits[j] = 0;
402       fine_priority[j] = ebits[j]<1;
403    }
404    RESTORE_STACK;
405    return codedBands;
406 }
407
408 int compute_allocation(const CELTMode *m, int start, int end, const int *offsets, int alloc_trim, int *intensity, int *dual_stereo,
409       int total, int *pulses, int *ebits, int *fine_priority, int _C, int LM, void *ec, int encode, int prev)
410 {
411    int lo, hi, len, j;
412    const int C = CHANNELS(_C);
413    int codedBands;
414    int skip_start;
415    int skip_rsv;
416    int intensity_rsv;
417    int dual_stereo_rsv;
418    VARDECL(int, bits1);
419    VARDECL(int, bits2);
420    VARDECL(int, thresh);
421    VARDECL(int, trim_offset);
422    SAVE_STACK;
423    
424    total = IMAX(total, 0);
425    len = m->nbEBands;
426    skip_start = start;
427    /* Reserve a bit to signal the end of manually skipped bands. */
428    skip_rsv = total >= 1<<BITRES ? 1<<BITRES : 0;
429    total -= skip_rsv;
430    /* Reserve bits for the intensity and dual stereo parameters. */
431    intensity_rsv = dual_stereo_rsv = 0;
432    if (C==2)
433    {
434       intensity_rsv = LOG2_FRAC_TABLE[end-start];
435       if (intensity_rsv>total)
436          intensity_rsv = 0;
437       else
438       {
439          total -= intensity_rsv;
440          dual_stereo_rsv = total>=1<<BITRES ? 1<<BITRES : 0;
441          total -= dual_stereo_rsv;
442       }
443    }
444    ALLOC(bits1, len, int);
445    ALLOC(bits2, len, int);
446    ALLOC(thresh, len, int);
447    ALLOC(trim_offset, len, int);
448
449    for (j=start;j<end;j++)
450    {
451       /* Below this threshold, we're sure not to allocate any PVQ bits */
452       thresh[j] = IMAX((C)<<BITRES, (3*(m->eBands[j+1]-m->eBands[j])<<LM<<BITRES)>>4);
453       /* Tilt of the allocation curve */
454       trim_offset[j] = C*(m->eBands[j+1]-m->eBands[j])*(alloc_trim-5-LM)*(m->nbEBands-j-1)
455             <<(LM+BITRES)>>6;
456       /* Giving less resolution to single-coefficient bands because they get
457          more benefit from having one coarse value per coefficient*/
458       if ((m->eBands[j+1]-m->eBands[j])<<LM==1)
459          trim_offset[j] -= C<<BITRES;
460    }
461    lo = 1;
462    hi = m->nbAllocVectors - 2;
463    do
464    {
465       int done = 0;
466       int psum = 0;
467       int mid = (lo+hi) >> 1;
468       for (j=end;j-->start;)
469       {
470          int N = m->eBands[j+1]-m->eBands[j];
471          bits1[j] = C*N*m->allocVectors[mid*len+j]<<LM>>2;
472          if (bits1[j] > 0)
473             bits1[j] = IMAX(0, bits1[j] + trim_offset[j]);
474          bits1[j] += offsets[j];
475          if (bits1[j] >= thresh[j] || done)
476          {
477             done = 1;
478             /* Don't allocate more than we can actually use */
479             psum += IMIN(bits1[j], 64*C<<BITRES<<LM);
480          } else {
481             if (bits1[j] >= C<<BITRES)
482                psum += C<<BITRES;
483          }
484       }
485       if (psum > total)
486          hi = mid - 1;
487       else
488          lo = mid + 1;
489       /*printf ("lo = %d, hi = %d\n", lo, hi);*/
490    }
491    while (lo <= hi);
492    hi = lo--;
493    /*printf ("interp between %d and %d\n", lo, hi);*/
494    for (j=start;j<end;j++)
495    {
496       int N = m->eBands[j+1]-m->eBands[j];
497       bits1[j] = C*N*m->allocVectors[lo*len+j]<<LM>>2;
498       bits2[j] = C*N*m->allocVectors[hi*len+j]<<LM>>2;
499       if (bits1[j] > 0)
500          bits1[j] = IMAX(0, bits1[j] + trim_offset[j]);
501       if (bits2[j] > 0)
502          bits2[j] = IMAX(0, bits2[j] + trim_offset[j]);
503       if (lo > 0)
504          bits1[j] += offsets[j];
505       bits2[j] += offsets[j];
506       if (offsets[j]>0)
507          skip_start = j;
508       bits2[j] -= bits1[j];
509    }
510    codedBands = interp_bits2pulses(m, start, end, skip_start, bits1, bits2, thresh,
511          total, skip_rsv, intensity, intensity_rsv, dual_stereo, dual_stereo_rsv,
512          pulses, ebits, fine_priority, C, LM, ec, encode, prev);
513    RESTORE_STACK;
514    return codedBands;
515 }
516