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