Move skip coding into interp_bits2pulses().
[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,
144       int *bits1, int *bits2, const int *thresh, int total, 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 force_skipping;
156    int done;
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=start;j<end;j++)
170       {
171          int tmp = bits1[j] + (mid*bits2[j]>>ALLOC_STEPS);
172          /* Don't allocate more than we can actually use */
173          if (tmp >= thresh[j] && !done)
174          {
175             psum += tmp;
176          } else {
177             done = 1;
178             if (tmp >= alloc_floor)
179                psum += alloc_floor;
180          }
181       }
182       if (psum > (total<<BITRES))
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=start;j<end;j++)
191    {
192       int tmp = bits1[j] + (lo*bits2[j]>>ALLOC_STEPS);
193       if (tmp < thresh[j] || done)
194       {
195          done = 1;
196          if (tmp >= alloc_floor)
197             tmp = alloc_floor;
198          else
199             tmp = 0;
200       }
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    force_skipping = 1;
208    codedBands=end;
209    for (j=end;j-->start;)
210    {
211       int band_width;
212       int band_bits;
213       int rem;
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<<BITRES)-psum;
217       percoeff = left/(m->eBands[codedBands]-m->eBands[start]);
218       left -= (m->eBands[codedBands]-m->eBands[start])*percoeff;
219       rem = IMAX(left-m->eBands[j],0);
220       band_width = m->eBands[codedBands]-m->eBands[j];
221       band_bits = bits[j] + percoeff*band_width + rem;
222       /*As long as, even after adding these bits, we're below the threshold for
223          this band, it is force-skipped.*/
224       force_skipping = force_skipping && band_bits < thresh[j];
225       if (!force_skipping)
226       {
227          /*Never skip the first band: we'd be coding a bit to signal that we're
228             going to waste all of the other bits.*/
229          if (j<=start)break;
230          /*If we have enough for the fine energy, but not more than a full
231             bit beyond that, or no more than one bit total, then don't bother
232             skipping this band: there's no extra bits to redistribute.*/
233          if ((band_bits>=alloc_floor && band_bits<=alloc_floor+(1<<BITRES))
234                || band_bits<(1<<BITRES))
235             break;
236          if (encode)
237          {
238             /*This if() block is the only part of the allocation function that
239                is not a mandatory part of the bitstream: any bands we choose to
240                skip here must be explicitly signaled.*/
241             /*Choose a threshold with some hysteresis to keep bands from
242                fluctuating in and out.*/
243             if (band_bits > ((j<prev?7:9)*band_width<<LM<<BITRES)>>4)
244             {
245                ec_enc_bit_prob((ec_enc *)ec, 1, 32768);
246                break;
247             }
248             ec_enc_bit_prob((ec_enc *)ec, 0, 32768);
249          } else if (ec_dec_bit_prob((ec_dec *)ec, 32768)) {
250             break;
251          }
252          /*We used a bit to skip this band.*/
253          psum += 1<<BITRES;
254          band_bits -= 1<<BITRES;
255       }
256       /*Reclaim the bits originally allocated to this band.*/
257       psum -= bits[j];
258       if (band_bits >= alloc_floor)
259       {
260          /*If we have enough for a fine energy bit per channel, use it.*/
261          psum += alloc_floor;
262          bits[j] = alloc_floor;
263       } else {
264          /*Otherwise this band gets nothing at all.*/
265          bits[j] = 0;
266       }
267       codedBands--;
268    }
269
270    /* Allocate the remaining bits */
271    if (codedBands>start) {
272       for (j=start;j<codedBands;j++)
273          bits[j] += percoeff*(m->eBands[j+1]-m->eBands[j]);
274       for (j=start;j<codedBands;j++)
275       {
276          int tmp = IMIN(left, m->eBands[j+1]-m->eBands[j]);
277          bits[j] += tmp;
278          left -= tmp;
279       }
280    }
281    /*for (j=0;j<end;j++)printf("%d ", bits[j]);printf("\n");*/
282    for (j=start;j<end;j++)
283    {
284       int N0, N, den;
285       int offset;
286       int NClogN;
287
288       celt_assert(bits[j] >= 0);
289       N0 = m->eBands[j+1]-m->eBands[j];
290       N=N0<<LM;
291       NClogN = N*C*(m->logN[j] + logM);
292
293       /* Compensate for the extra DoF in stereo */
294       den=(C*N+ ((C==2 && N>2) ? 1 : 0));
295
296       /* Offset for the number of fine bits by log2(N)/2 + FINE_OFFSET
297          compared to their "fair share" of total/N */
298       offset = (NClogN>>1)-N*C*FINE_OFFSET;
299
300       /* N=2 is the only point that doesn't match the curve */
301       if (N==2)
302          offset += N*C<<BITRES>>2;
303
304       /* Changing the offset for allocating the second and third fine energy bit */
305       if (bits[j] + offset < den*2<<BITRES)
306          offset += NClogN>>2;
307       else if (bits[j] + offset < den*3<<BITRES)
308          offset += NClogN>>3;
309
310       /* Divide with rounding */
311       ebits[j] = IMAX(0, (bits[j] + offset + (den<<(BITRES-1))) / (den<<BITRES));
312
313       /* If we rounded down, make it a candidate for final fine energy pass */
314       fine_priority[j] = ebits[j]*(den<<BITRES) >= bits[j]+offset;
315
316       /* For N=1, all bits go to fine energy except for a single sign bit */
317       if (N==1)
318       {
319          ebits[j] = IMAX(0,(bits[j]/C >> BITRES)-1);
320          fine_priority[j] = (ebits[j]+1)*C<<BITRES >= bits[j];
321       }
322       /* Make sure not to bust */
323       if (C*ebits[j] > (bits[j]>>BITRES))
324          ebits[j] = bits[j]/C >> BITRES;
325
326       /* More than that is useless because that's about as far as PVQ can go */
327       if (ebits[j]>7)
328          ebits[j]=7;
329
330       /* The other bits are assigned to PVQ */
331       bits[j] -= C*ebits[j]<<BITRES;
332       celt_assert(bits[j] >= 0);
333       celt_assert(ebits[j] >= 0);
334    }
335    RESTORE_STACK;
336    return codedBands;
337 }
338
339 int compute_allocation(const CELTMode *m, int start, int end, int *offsets, int alloc_trim,
340       int total, int *pulses, int *ebits, int *fine_priority, int _C, int LM, void *ec, int encode, int prev)
341 {
342    int lo, hi, len, j;
343    const int C = CHANNELS(_C);
344    int codedBands;
345    VARDECL(int, bits1);
346    VARDECL(int, bits2);
347    VARDECL(int, thresh);
348    VARDECL(int, trim_offset);
349    SAVE_STACK;
350    
351    total = IMAX(total, 0);
352    len = m->nbEBands;
353    ALLOC(bits1, len, int);
354    ALLOC(bits2, len, int);
355    ALLOC(thresh, len, int);
356    ALLOC(trim_offset, len, int);
357
358    /* Below this threshold, we're sure not to allocate any PVQ bits */
359    for (j=start;j<end;j++)
360       thresh[j] = IMAX((C)<<BITRES, (3*(m->eBands[j+1]-m->eBands[j])<<LM<<BITRES)>>4);
361    /* Tilt of the allocation curve */
362    for (j=start;j<end;j++)
363       trim_offset[j] = C*(m->eBands[j+1]-m->eBands[j])*(alloc_trim-5-LM)*(m->nbEBands-j-1)
364             <<(LM+BITRES)>>6;
365
366    lo = 0;
367    hi = m->nbAllocVectors - 1;
368    while (hi-lo != 1)
369    {
370       int psum = 0;
371       int mid = (lo+hi) >> 1;
372       for (j=start;j<end;j++)
373       {
374          int N = m->eBands[j+1]-m->eBands[j];
375          bits1[j] = C*N*m->allocVectors[mid*len+j]<<LM>>2;
376          if (bits1[j] > 0)
377             bits1[j] += trim_offset[j];
378          if (bits1[j] < 0)
379             bits1[j] = 0;
380          bits1[j] += offsets[j];
381          if (bits1[j] >= thresh[j])
382             psum += bits1[j];
383          else if (bits1[j] >= C<<BITRES)
384             psum += C<<BITRES;
385
386          /*printf ("%d ", bits[j]);*/
387       }
388       /*printf ("\n");*/
389       if (psum > (total<<BITRES))
390          hi = mid;
391       else
392          lo = mid;
393       /*printf ("lo = %d, hi = %d\n", lo, hi);*/
394    }
395    /*printf ("interp between %d and %d\n", lo, hi);*/
396    for (j=start;j<end;j++)
397    {
398       int N = m->eBands[j+1]-m->eBands[j];
399       bits1[j] = (C*N*m->allocVectors[lo*len+j]<<LM>>2);
400       bits2[j] = (C*N*m->allocVectors[hi*len+j]<<LM>>2) - bits1[j];
401       if (bits1[j] > 0)
402          bits1[j] += trim_offset[j];
403       if (bits1[j] < 0)
404          bits1[j] = 0;
405       bits1[j] += offsets[j];
406    }
407    codedBands = interp_bits2pulses(m, start, end, bits1, bits2, thresh,
408          total, pulses, ebits, fine_priority, len, C, LM, ec, encode, prev);
409    RESTORE_STACK;
410    return codedBands;
411 }
412