Refactor ec_enc_patch_initial_bits().
[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    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
17    ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
18    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
19    A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE FOUNDATION OR
20    CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
21    EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
22    PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
23    PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
24    LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
25    NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
26    SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
27 */
28
29 #ifdef HAVE_CONFIG_H
30 #include "config.h"
31 #endif
32
33 #include <math.h>
34 #include "modes.h"
35 #include "cwrs.h"
36 #include "arch.h"
37 #include "os_support.h"
38
39 #include "entcode.h"
40 #include "rate.h"
41
42
43 static const unsigned char LOG2_FRAC_TABLE[24]={
44    0,
45    8,13,
46   16,19,21,23,
47   24,26,27,28,29,30,31,32,
48   32,33,34,34,35,36,36,37,37
49 };
50
51 #ifdef CUSTOM_MODES
52
53 /*Determines if V(N,K) fits in a 32-bit unsigned integer.
54   N and K are themselves limited to 15 bits.*/
55 static int fits_in32(int _n, int _k)
56 {
57    static const celt_int16 maxN[15] = {
58       32767, 32767, 32767, 1476, 283, 109,  60,  40,
59        29,  24,  20,  18,  16,  14,  13};
60    static const celt_int16 maxK[15] = {
61       32767, 32767, 32767, 32767, 1172, 238,  95,  53,
62        36,  27,  22,  18,  16,  15,  13};
63    if (_n>=14)
64    {
65       if (_k>=14)
66          return 0;
67       else
68          return _n <= maxN[_k];
69    } else {
70       return _k <= maxK[_n];
71    }
72 }
73
74 void compute_pulse_cache(CELTMode *m, int LM)
75 {
76    int C;
77    int i;
78    int j;
79    int curr=0;
80    int nbEntries=0;
81    int entryN[100], entryK[100], entryI[100];
82    const celt_int16 *eBands = m->eBands;
83    PulseCache *cache = &m->cache;
84    celt_int16 *cindex;
85    unsigned char *bits;
86    unsigned char *cap;
87
88    cindex = celt_alloc(sizeof(cache->index[0])*m->nbEBands*(LM+2));
89    cache->index = cindex;
90
91    /* Scan for all unique band sizes */
92    for (i=0;i<=LM+1;i++)
93    {
94       for (j=0;j<m->nbEBands;j++)
95       {
96          int k;
97          int N = (eBands[j+1]-eBands[j])<<i>>1;
98          cindex[i*m->nbEBands+j] = -1;
99          /* Find other bands that have the same size */
100          for (k=0;k<=i;k++)
101          {
102             int n;
103             for (n=0;n<m->nbEBands && (k!=i || n<j);n++)
104             {
105                if (N == (eBands[n+1]-eBands[n])<<k>>1)
106                {
107                   cindex[i*m->nbEBands+j] = cindex[k*m->nbEBands+n];
108                   break;
109                }
110             }
111          }
112          if (cache->index[i*m->nbEBands+j] == -1 && N!=0)
113          {
114             int K;
115             entryN[nbEntries] = N;
116             K = 0;
117             while (fits_in32(N,get_pulses(K+1)) && K<MAX_PSEUDO)
118                K++;
119             entryK[nbEntries] = K;
120             cindex[i*m->nbEBands+j] = curr;
121             entryI[nbEntries] = curr;
122
123             curr += K+1;
124             nbEntries++;
125          }
126       }
127    }
128    bits = celt_alloc(sizeof(unsigned char)*curr);
129    cache->bits = bits;
130    cache->size = curr;
131    /* Compute the cache for all unique sizes */
132    for (i=0;i<nbEntries;i++)
133    {
134       unsigned char *ptr = bits+entryI[i];
135       celt_int16 tmp[MAX_PULSES+1];
136       get_required_bits(tmp, entryN[i], get_pulses(entryK[i]), BITRES);
137       for (j=1;j<=entryK[i];j++)
138          ptr[j] = tmp[get_pulses(j)]-1;
139       ptr[0] = entryK[i];
140    }
141
142    /* Compute the maximum rate for each band at which we'll reliably use as
143        many bits as we ask for. */
144    cache->caps = cap = celt_alloc(sizeof(cache->caps[0])*(LM+1)*2*m->nbEBands);
145    for (i=0;i<=LM;i++)
146    {
147       for (C=1;C<=2;C++)
148       {
149          for (j=0;j<m->nbEBands;j++)
150          {
151             int N0;
152             int max_bits;
153             N0 = m->eBands[j+1]-m->eBands[j];
154             /* N=1 bands only have a sign bit and fine bits. */
155             if (N0<<i == 1)
156                max_bits = C*(1+MAX_FINE_BITS)<<BITRES;
157             else
158             {
159                const unsigned char *pcache;
160                celt_int32           num;
161                celt_int32           den;
162                int                  LM0;
163                int                  N;
164                int                  offset;
165                int                  ndof;
166                int                  qb;
167                int                  k;
168                LM0 = 0;
169                /* Even-sized bands bigger than N=2 can be split one more
170                    time. */
171                if (N0 > 2 && !(N0&1))
172                {
173                   N0>>=1;
174                   LM0--;
175                }
176                /* N0=1 bands can't be split down to N<2. */
177                else if (N0 <= 1)
178                {
179                   LM0=IMIN(i,1);
180                   N0<<=LM0;
181                }
182                /* Compute the cost for the lowest-level PVQ of a fully split
183                    band. */
184                pcache = bits + cindex[(LM0+1)*m->nbEBands+j];
185                max_bits = pcache[pcache[0]]+1;
186                /* Add in the cost of coding regular splits. */
187                N = N0;
188                for(k=0;k<i-LM0;k++){
189                   max_bits <<= 1;
190                   /* Offset the number of qtheta bits by log2(N)/2
191                       + QTHETA_OFFSET compared to their "fair share" of
192                       total/N */
193                   offset = (m->logN[j]+(LM0+k<<BITRES)>>1)-QTHETA_OFFSET;
194                   /* The number of qtheta bits we'll allocate if the remainder
195                       is to be max_bits.
196                      The average measured cost for theta is 0.89701 times qb,
197                       approximated here as 459/512. */
198                   num=459*(celt_int32)((2*N-1)*offset+max_bits);
199                   den=((celt_int32)(2*N-1)<<9)-459;
200                   qb = IMIN((num+(den>>1))/den, 57);
201                   celt_assert(qb >= 0);
202                   max_bits += qb;
203                   N <<= 1;
204                }
205                /* Add in the cost of a stereo split, if necessary. */
206                if (C==2)
207                {
208                   max_bits <<= 1;
209                   offset = (m->logN[j]+(i<<BITRES)>>1)-(N==2?QTHETA_OFFSET_TWOPHASE:QTHETA_OFFSET);
210                   ndof = 2*N-1-(N==2);
211                   /* The average measured cost for theta with the step PDF is
212                       0.95164 times qb, approximated here as 487/512. */
213                   num = (N==2?512:487)*(celt_int32)(max_bits+ndof*offset);
214                   den = ((celt_int32)ndof<<9)-(N==2?512:487);
215                   qb = IMIN((num+(den>>1))/den, (N==2?64:61));
216                   celt_assert(qb >= 0);
217                   max_bits += qb;
218                }
219                /* Add the fine bits we'll use. */
220                /* Compensate for the extra DoF in stereo */
221                ndof = C*N + ((C==2 && N>2) ? 1 : 0);
222                /* Offset the number of fine bits by log2(N)/2 + FINE_OFFSET
223                    compared to their "fair share" of total/N */
224                offset = (m->logN[j] + (i<<BITRES)>>1)-FINE_OFFSET;
225                /* N=2 is the only point that doesn't match the curve */
226                if (N==2)
227                   offset += 1<<BITRES>>2;
228                /* The number of fine bits we'll allocate if the remainder is
229                    to be max_bits. */
230                num = max_bits+ndof*offset;
231                den = ndof-1<<BITRES;
232                qb = IMIN((num+(den>>1))/den, MAX_FINE_BITS);
233                celt_assert(qb >= 0);
234                max_bits += C*qb<<BITRES;
235             }
236             max_bits = (4*max_bits/(C*(m->eBands[j+1]-m->eBands[j]<<i)))-64;
237             celt_assert(max_bits >= 0);
238             celt_assert(max_bits < 256);
239             *cap++ = (unsigned char)max_bits;
240          }
241       }
242    }
243 }
244
245 #endif /* CUSTOM_MODES */
246
247
248 #define ALLOC_STEPS 6
249
250 static inline int interp_bits2pulses(const CELTMode *m, int start, int end, int skip_start,
251       const int *bits1, const int *bits2, const int *thresh, const int *cap, celt_int32 total, celt_int32 *_balance,
252       int skip_rsv, int *intensity, int intensity_rsv, int *dual_stereo, int dual_stereo_rsv, int *bits,
253       int *ebits, int *fine_priority, int _C, int LM, ec_ctx *ec, int encode, int prev)
254 {
255    celt_int32 psum;
256    int lo, hi;
257    int i, j;
258    int logM;
259    const int C = CHANNELS(_C);
260    int stereo;
261    int codedBands=-1;
262    int alloc_floor;
263    celt_int32 left, percoeff;
264    int done;
265    int balance;
266    SAVE_STACK;
267
268    alloc_floor = C<<BITRES;
269    stereo = C>1;
270
271    logM = LM<<BITRES;
272    lo = 0;
273    hi = 1<<ALLOC_STEPS;
274    for (i=0;i<ALLOC_STEPS;i++)
275    {
276       int mid = (lo+hi)>>1;
277       psum = 0;
278       done = 0;
279       for (j=end;j-->start;)
280       {
281          int tmp = bits1[j] + (mid*(celt_int32)bits2[j]>>ALLOC_STEPS);
282          if (tmp >= thresh[j] || done)
283          {
284             done = 1;
285             /* Don't allocate more than we can actually use */
286             psum += IMIN(tmp, cap[j]);
287          } else {
288             if (tmp >= alloc_floor)
289                psum += alloc_floor;
290          }
291       }
292       if (psum > total)
293          hi = mid;
294       else
295          lo = mid;
296    }
297    psum = 0;
298    /*printf ("interp bisection gave %d\n", lo);*/
299    done = 0;
300    for (j=end;j-->start;)
301    {
302       int tmp = bits1[j] + (lo*bits2[j]>>ALLOC_STEPS);
303       if (tmp < thresh[j] && !done)
304       {
305          if (tmp >= alloc_floor)
306             tmp = alloc_floor;
307          else
308             tmp = 0;
309       } else
310          done = 1;
311       /* Don't allocate more than we can actually use */
312       tmp = IMIN(tmp, cap[j]);
313       bits[j] = tmp;
314       psum += tmp;
315    }
316
317    /* Decide which bands to skip, working backwards from the end. */
318    for (codedBands=end;;codedBands--)
319    {
320       int band_width;
321       int band_bits;
322       int rem;
323       j = codedBands-1;
324       /* Never skip the first band, nor a band that has been boosted by
325           dynalloc.
326          In the first case, we'd be coding a bit to signal we're going to waste
327           all the other bits.
328          In the second case, we'd be coding a bit to redistribute all the bits
329           we just signaled should be cocentrated in this band. */
330       if (j<=skip_start)
331       {
332          /* Give the bit we reserved to end skipping back. */
333          total += skip_rsv;
334          break;
335       }
336       /*Figure out how many left-over bits we would be adding to this band.
337         This can include bits we've stolen back from higher, skipped bands.*/
338       left = total-psum;
339       percoeff = left/(m->eBands[codedBands]-m->eBands[start]);
340       left -= (m->eBands[codedBands]-m->eBands[start])*percoeff;
341       rem = IMAX(left-(m->eBands[j]-m->eBands[start]),0);
342       band_width = m->eBands[codedBands]-m->eBands[j];
343       band_bits = (int)(bits[j] + percoeff*band_width + rem);
344       /*Only code a skip decision if we're above the threshold for this band.
345         Otherwise it is force-skipped.
346         This ensures that we have enough bits to code the skip flag.*/
347       if (band_bits >= IMAX(thresh[j], alloc_floor+(1<<BITRES)))
348       {
349          if (encode)
350          {
351             /*This if() block is the only part of the allocation function that
352                is not a mandatory part of the bitstream: any bands we choose to
353                skip here must be explicitly signaled.*/
354             /*Choose a threshold with some hysteresis to keep bands from
355                fluctuating in and out.*/
356             if (band_bits > ((j<prev?7:9)*band_width<<LM<<BITRES)>>4)
357             {
358                ec_enc_bit_logp(ec, 1, 1);
359                break;
360             }
361             ec_enc_bit_logp(ec, 0, 1);
362          } else if (ec_dec_bit_logp(ec, 1)) {
363             break;
364          }
365          /*We used a bit to skip this band.*/
366          psum += 1<<BITRES;
367          band_bits -= 1<<BITRES;
368       }
369       /*Reclaim the bits originally allocated to this band.*/
370       psum -= bits[j]+intensity_rsv;
371       if (intensity_rsv > 0)
372          intensity_rsv = LOG2_FRAC_TABLE[j-start];
373       psum += intensity_rsv;
374       if (band_bits >= alloc_floor)
375       {
376          /*If we have enough for a fine energy bit per channel, use it.*/
377          psum += alloc_floor;
378          bits[j] = alloc_floor;
379       } else {
380          /*Otherwise this band gets nothing at all.*/
381          bits[j] = 0;
382       }
383    }
384
385    celt_assert(codedBands > start);
386    /* Code the intensity and dual stereo parameters. */
387    if (intensity_rsv > 0)
388    {
389       if (encode)
390       {
391          *intensity = IMIN(*intensity, codedBands);
392          ec_enc_uint(ec, *intensity-start, codedBands+1-start);
393       }
394       else
395          *intensity = start+ec_dec_uint(ec, codedBands+1-start);
396    }
397    else
398       *intensity = 0;
399    if (*intensity <= start)
400    {
401       total += dual_stereo_rsv;
402       dual_stereo_rsv = 0;
403    }
404    if (dual_stereo_rsv > 0)
405    {
406       if (encode)
407          ec_enc_bit_logp(ec, *dual_stereo, 1);
408       else
409          *dual_stereo = ec_dec_bit_logp(ec, 1);
410    }
411    else
412       *dual_stereo = 0;
413
414    /* Allocate the remaining bits */
415    left = total-psum;
416    percoeff = left/(m->eBands[codedBands]-m->eBands[start]);
417    left -= (m->eBands[codedBands]-m->eBands[start])*percoeff;
418    for (j=start;j<codedBands;j++)
419       bits[j] += ((int)percoeff*(m->eBands[j+1]-m->eBands[j]));
420    for (j=start;j<codedBands;j++)
421    {
422       int tmp = (int)IMIN(left, m->eBands[j+1]-m->eBands[j]);
423       bits[j] += tmp;
424       left -= tmp;
425    }
426    /*for (j=0;j<end;j++)printf("%d ", bits[j]);printf("\n");*/
427
428    balance = 0;
429    for (j=start;j<codedBands;j++)
430    {
431       int N0, N, den;
432       int offset;
433       int NClogN;
434       int excess;
435
436       celt_assert(bits[j] >= 0);
437       N0 = m->eBands[j+1]-m->eBands[j];
438       N=N0<<LM;
439       bits[j] += balance;
440
441       if (N>1)
442       {
443          excess = IMAX(bits[j]-cap[j],0);
444          bits[j] -= excess;
445
446          /* Compensate for the extra DoF in stereo */
447          den=(C*N+ ((C==2 && N>2 && !*dual_stereo && j<*intensity) ? 1 : 0));
448
449          NClogN = den*(m->logN[j] + logM);
450
451          /* Offset for the number of fine bits by log2(N)/2 + FINE_OFFSET
452             compared to their "fair share" of total/N */
453          offset = (NClogN>>1)-den*FINE_OFFSET;
454
455          /* N=2 is the only point that doesn't match the curve */
456          if (N==2)
457             offset += den<<BITRES>>2;
458
459          /* Changing the offset for allocating the second and third
460              fine energy bit */
461          if (bits[j] + offset < den*2<<BITRES)
462             offset += NClogN>>2;
463          else if (bits[j] + offset < den*3<<BITRES)
464             offset += NClogN>>3;
465
466          /* Divide with rounding */
467          ebits[j] = IMAX(0, (bits[j] + offset + (den<<(BITRES-1))) / (den<<BITRES));
468
469          /* Make sure not to bust */
470          if (C*ebits[j] > (bits[j]>>BITRES))
471             ebits[j] = bits[j] >> stereo >> BITRES;
472
473          /* More than that is useless because that's about as far as PVQ can go */
474          ebits[j] = IMIN(ebits[j], MAX_FINE_BITS);
475
476          /* If we rounded down or capped this band, make it a candidate for the
477              final fine energy pass */
478          fine_priority[j] = ebits[j]*(den<<BITRES) >= bits[j]+offset;
479
480          /* Remove the allocated fine bits; the rest are assigned to PVQ */
481          bits[j] -= C*ebits[j]<<BITRES;
482
483       } else {
484          /* For N=1, all bits go to fine energy except for a single sign bit */
485          excess = IMAX(0,bits[j]-(C<<BITRES));
486          bits[j] -= excess;
487          ebits[j] = 0;
488          fine_priority[j] = 1;
489       }
490
491       /* Fine energy can't take advantage of the re-balancing in
492           quant_all_bands().
493          Instead, do the re-balancing here.*/
494       if(excess > 0)
495       {
496          int extra_fine;
497          int extra_bits;
498          extra_fine = IMIN(excess >> stereo+BITRES, MAX_FINE_BITS-ebits[j]);
499          ebits[j] += extra_fine;
500          extra_bits = extra_fine*C<<BITRES;
501          fine_priority[j] = extra_bits >= excess-balance;
502          excess -= extra_bits;
503       }
504       balance = excess;
505
506       celt_assert(bits[j] >= 0);
507       celt_assert(ebits[j] >= 0);
508    }
509    /* Save any remaining bits over the cap for the rebalancing in
510        quant_all_bands(). */
511    *_balance = balance;
512
513    /* The skipped bands use all their bits for fine energy. */
514    for (;j<end;j++)
515    {
516       ebits[j] = bits[j] >> stereo >> BITRES;
517       celt_assert(C*ebits[j]<<BITRES == bits[j]);
518       bits[j] = 0;
519       fine_priority[j] = ebits[j]<1;
520    }
521    RESTORE_STACK;
522    return codedBands;
523 }
524
525 int compute_allocation(const CELTMode *m, int start, int end, const int *offsets, const int *cap, int alloc_trim, int *intensity, int *dual_stereo,
526       int total, celt_int32 *balance, int *pulses, int *ebits, int *fine_priority, int _C, int LM, ec_ctx *ec, int encode, int prev)
527 {
528    int lo, hi, len, j;
529    const int C = CHANNELS(_C);
530    int codedBands;
531    int skip_start;
532    int skip_rsv;
533    int intensity_rsv;
534    int dual_stereo_rsv;
535    VARDECL(int, bits1);
536    VARDECL(int, bits2);
537    VARDECL(int, thresh);
538    VARDECL(int, trim_offset);
539    SAVE_STACK;
540    
541    total = IMAX(total, 0);
542    len = m->nbEBands;
543    skip_start = start;
544    /* Reserve a bit to signal the end of manually skipped bands. */
545    skip_rsv = total >= 1<<BITRES ? 1<<BITRES : 0;
546    total -= skip_rsv;
547    /* Reserve bits for the intensity and dual stereo parameters. */
548    intensity_rsv = dual_stereo_rsv = 0;
549    if (C==2)
550    {
551       intensity_rsv = LOG2_FRAC_TABLE[end-start];
552       if (intensity_rsv>total)
553          intensity_rsv = 0;
554       else
555       {
556          total -= intensity_rsv;
557          dual_stereo_rsv = total>=1<<BITRES ? 1<<BITRES : 0;
558          total -= dual_stereo_rsv;
559       }
560    }
561    ALLOC(bits1, len, int);
562    ALLOC(bits2, len, int);
563    ALLOC(thresh, len, int);
564    ALLOC(trim_offset, len, int);
565
566    for (j=start;j<end;j++)
567    {
568       /* Below this threshold, we're sure not to allocate any PVQ bits */
569       thresh[j] = IMAX((C)<<BITRES, (3*(m->eBands[j+1]-m->eBands[j])<<LM<<BITRES)>>4);
570       /* Tilt of the allocation curve */
571       trim_offset[j] = C*(m->eBands[j+1]-m->eBands[j])*(alloc_trim-5-LM)*(end-j-1)
572             <<(LM+BITRES)>>6;
573       /* Giving less resolution to single-coefficient bands because they get
574          more benefit from having one coarse value per coefficient*/
575       if ((m->eBands[j+1]-m->eBands[j])<<LM==1)
576          trim_offset[j] -= C<<BITRES;
577    }
578    lo = 1;
579    hi = m->nbAllocVectors - 1;
580    do
581    {
582       int done = 0;
583       int psum = 0;
584       int mid = (lo+hi) >> 1;
585       for (j=end;j-->start;)
586       {
587          int N = m->eBands[j+1]-m->eBands[j];
588          bits1[j] = C*N*m->allocVectors[mid*len+j]<<LM>>2;
589          if (bits1[j] > 0)
590             bits1[j] = IMAX(0, bits1[j] + trim_offset[j]);
591          bits1[j] += offsets[j];
592          if (bits1[j] >= thresh[j] || done)
593          {
594             done = 1;
595             /* Don't allocate more than we can actually use */
596             psum += IMIN(bits1[j], cap[j]);
597          } else {
598             if (bits1[j] >= C<<BITRES)
599                psum += C<<BITRES;
600          }
601       }
602       if (psum > total)
603          hi = mid - 1;
604       else
605          lo = mid + 1;
606       /*printf ("lo = %d, hi = %d\n", lo, hi);*/
607    }
608    while (lo <= hi);
609    hi = lo--;
610    /*printf ("interp between %d and %d\n", lo, hi);*/
611    for (j=start;j<end;j++)
612    {
613       int N = m->eBands[j+1]-m->eBands[j];
614       bits1[j] = C*N*m->allocVectors[lo*len+j]<<LM>>2;
615       bits2[j] = hi>=m->nbAllocVectors ?
616             cap[j] : C*N*m->allocVectors[hi*len+j]<<LM>>2;
617       if (bits1[j] > 0)
618          bits1[j] = IMAX(0, bits1[j] + trim_offset[j]);
619       if (bits2[j] > 0)
620          bits2[j] = IMAX(0, bits2[j] + trim_offset[j]);
621       if (lo > 0)
622          bits1[j] += offsets[j];
623       bits2[j] += offsets[j];
624       if (offsets[j]>0)
625          skip_start = j;
626       bits2[j] = IMAX(0,bits2[j]-bits1[j]);
627    }
628    codedBands = interp_bits2pulses(m, start, end, skip_start, bits1, bits2, thresh, cap,
629          total, balance, skip_rsv, intensity, intensity_rsv, dual_stereo, dual_stereo_rsv,
630          pulses, ebits, fine_priority, C, LM, ec, encode, prev);
631    RESTORE_STACK;
632    return codedBands;
633 }
634