Max delta: +/- 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 len, 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 8 is useless because that's about as far as PVQ can go */
370          if (ebits[j]>8)
371             ebits[j]=8;
372
373          /* If we rounded down or capped this band, make it a candidate for the
374              final fine energy pass */
375          fine_priority[j] = ebits[j]*(den<<BITRES) >= bits[j]+offset;
376
377       } else {
378          /* For N=1, all bits go to fine energy except for a single sign bit */
379          ebits[j] = IMIN(IMAX(0,(bits[j] >> stereo >> BITRES)-1),8);
380          fine_priority[j] = (ebits[j]+1)*C<<BITRES >= (bits[j]-balance);
381          /* N=1 bands can't take advantage of the re-balancing in
382              quant_all_bands() because they don't have shape, only fine energy.
383             Instead, do the re-balancing here.*/
384          balance = IMAX(0,bits[j] - ((ebits[j]+1)*C<<BITRES));
385          if (j+1<codedBands)
386          {
387             bits[j] -= balance;
388             bits[j+1] += balance;
389          }
390       }
391
392       /* The other bits are assigned to PVQ */
393       bits[j] -= C*ebits[j]<<BITRES;
394       celt_assert(bits[j] >= 0);
395       celt_assert(ebits[j] >= 0);
396    }
397    /* The skipped bands use all their bits for fine energy. */
398    for (;j<end;j++)
399    {
400       ebits[j] = bits[j] >> stereo >> BITRES;
401       celt_assert(C*ebits[j]<<BITRES == bits[j]);
402       bits[j] = 0;
403       fine_priority[j] = ebits[j]<1;
404    }
405    RESTORE_STACK;
406    return codedBands;
407 }
408
409 int compute_allocation(const CELTMode *m, int start, int end, const int *offsets, int alloc_trim, int *intensity, int *dual_stereo,
410       int total, int *pulses, int *ebits, int *fine_priority, int _C, int LM, void *ec, int encode, int prev)
411 {
412    int lo, hi, len, j;
413    const int C = CHANNELS(_C);
414    int codedBands;
415    int skip_start;
416    int skip_rsv;
417    int intensity_rsv;
418    int dual_stereo_rsv;
419    VARDECL(int, bits1);
420    VARDECL(int, bits2);
421    VARDECL(int, thresh);
422    VARDECL(int, trim_offset);
423    SAVE_STACK;
424    
425    total = IMAX(total, 0);
426    len = m->nbEBands;
427    skip_start = start;
428    /* Reserve a bit to signal the end of manually skipped bands. */
429    skip_rsv = total >= 1<<BITRES ? 1<<BITRES : 0;
430    total -= skip_rsv;
431    /* Reserve bits for the intensity and dual stereo parameters. */
432    intensity_rsv = dual_stereo_rsv = 0;
433    if (C==2)
434    {
435       intensity_rsv = LOG2_FRAC_TABLE[end-start];
436       if (intensity_rsv>total)
437          intensity_rsv = 0;
438       else
439       {
440          total -= intensity_rsv;
441          dual_stereo_rsv = total>=1<<BITRES ? 1<<BITRES : 0;
442          total -= dual_stereo_rsv;
443       }
444    }
445    ALLOC(bits1, len, int);
446    ALLOC(bits2, len, int);
447    ALLOC(thresh, len, int);
448    ALLOC(trim_offset, len, int);
449
450    for (j=start;j<end;j++)
451    {
452       /* Below this threshold, we're sure not to allocate any PVQ bits */
453       thresh[j] = IMAX((C)<<BITRES, (3*(m->eBands[j+1]-m->eBands[j])<<LM<<BITRES)>>4);
454       /* Tilt of the allocation curve */
455       trim_offset[j] = C*(m->eBands[j+1]-m->eBands[j])*(alloc_trim-5-LM)*(m->nbEBands-j-1)
456             <<(LM+BITRES)>>6;
457       /* Giving less resolution to single-coefficient bands because they get
458          more benefit from having one coarse value per coefficient*/
459       if ((m->eBands[j+1]-m->eBands[j])<<LM==1)
460          trim_offset[j] -= C<<BITRES;
461    }
462    lo = 1;
463    hi = m->nbAllocVectors - 2;
464    do
465    {
466       int done = 0;
467       int psum = 0;
468       int mid = (lo+hi) >> 1;
469       for (j=end;j-->start;)
470       {
471          int N = m->eBands[j+1]-m->eBands[j];
472          bits1[j] = C*N*m->allocVectors[mid*len+j]<<LM>>2;
473          if (bits1[j] > 0)
474             bits1[j] = IMAX(0, bits1[j] + trim_offset[j]);
475          bits1[j] += offsets[j];
476          if (bits1[j] >= thresh[j] || done)
477          {
478             done = 1;
479             /* Don't allocate more than we can actually use */
480             psum += IMIN(bits1[j], 64*C<<BITRES<<LM);
481          } else {
482             if (bits1[j] >= C<<BITRES)
483                psum += C<<BITRES;
484          }
485       }
486       if (psum > total)
487          hi = mid - 1;
488       else
489          lo = mid + 1;
490       /*printf ("lo = %d, hi = %d\n", lo, hi);*/
491    }
492    while (lo <= hi);
493    hi = lo--;
494    /*printf ("interp between %d and %d\n", lo, hi);*/
495    for (j=start;j<end;j++)
496    {
497       int N = m->eBands[j+1]-m->eBands[j];
498       bits1[j] = C*N*m->allocVectors[lo*len+j]<<LM>>2;
499       bits2[j] = C*N*m->allocVectors[hi*len+j]<<LM>>2;
500       if (bits1[j] > 0)
501          bits1[j] = IMAX(0, bits1[j] + trim_offset[j]);
502       if (bits2[j] > 0)
503          bits2[j] = IMAX(0, bits2[j] + trim_offset[j]);
504       if (lo > 0)
505          bits1[j] += offsets[j];
506       bits2[j] += offsets[j];
507       if (offsets[j]>0)
508          skip_start = j;
509       bits2[j] -= bits1[j];
510    }
511    codedBands = interp_bits2pulses(m, start, end, skip_start, bits1, bits2, thresh,
512          total, skip_rsv, intensity, intensity_rsv, dual_stereo, dual_stereo_rsv,
513          pulses, ebits, fine_priority, len, C, LM, ec, encode, prev);
514    RESTORE_STACK;
515    return codedBands;
516 }
517