Two-pass algorithm for filling the remaining bits. Still not perfect, but close
[opus.git] / libcelt / rate.c
1 /* (C) 2007 Jean-Marc Valin, CSIRO
2 */
3 /*
4    Redistribution and use in source and binary forms, with or without
5    modification, are permitted provided that the following conditions
6    are met:
7    
8    - Redistributions of source code must retain the above copyright
9    notice, this list of conditions and the following disclaimer.
10    
11    - Redistributions in binary form must reproduce the above copyright
12    notice, this list of conditions and the following disclaimer in the
13    documentation and/or other materials provided with the distribution.
14    
15    - Neither the name of the Xiph.org Foundation nor the names of its
16    contributors may be used to endorse or promote products derived from
17    this software without specific prior written permission.
18    
19    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
20    ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
21    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
22    A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE FOUNDATION OR
23    CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
24    EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
25    PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
26    PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
27    LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
28    NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
29    SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
30 */
31
32 #include <math.h>
33 #include "modes.h"
34 #include "cwrs.h"
35 #include "arch.h"
36 #include "os_support.h"
37
38 #include "entcode.h"
39 #include "rate.h"
40
41 #define BITRES 4
42 #define BITROUND 8
43 #define BITOVERFLOW 10000
44
45 #define MAX_PULSES 64
46
47 int log2_frac(ec_uint32 val, int frac)
48 {
49    int i;
50    /* EC_ILOG() actually returns log2()+1, go figure */
51    int L = EC_ILOG(val)-1;
52    //printf ("in: %d %d ", val, L);
53    if (L>14)
54       val >>= L-14;
55    else if (L<14)
56       val <<= 14-L;
57    L <<= frac;
58    //printf ("%d\n", val);
59    for (i=0;i<frac;i++)
60    {
61       val = (val*val) >> 15;
62       //printf ("%d\n", val);
63       if (val > 16384)
64          L |= (1<<(frac-i-1));
65       else   
66          val <<= 1;
67    }
68    return L;
69 }
70
71 int log2_frac64(ec_uint64 val, int frac)
72 {
73    int i;
74    /* EC_ILOG64() actually returns log2()+1, go figure */
75    int L = EC_ILOG64(val)-1;
76    //printf ("in: %d %d ", val, L);
77    if (L>14)
78       val >>= L-14;
79    else if (L<14)
80       val <<= 14-L;
81    L <<= frac;
82    //printf ("%d\n", val);
83    for (i=0;i<frac;i++)
84    {
85       val = (val*val) >> 15;
86       //printf ("%d\n", val);
87       if (val > 16384)
88          L |= (1<<(frac-i-1));
89       else   
90          val <<= 1;
91    }
92    return L;
93 }
94
95
96 void alloc_init(struct alloc_data *alloc, const CELTMode *m)
97 {
98    int i, prevN, BC;
99    const int *eBands = m->eBands;
100    
101    alloc->mode = m;
102    alloc->len = m->nbEBands;
103    alloc->bands = m->eBands;
104    alloc->bits = celt_alloc(m->nbEBands*sizeof(int*));
105    
106    BC = m->nbMdctBlocks*m->nbChannels;
107    prevN = -1;
108    for (i=0;i<alloc->len;i++)
109    {
110       int N = BC*(eBands[i+1]-eBands[i]);
111       if (N == prevN && eBands[i] < m->pitchEnd)
112       {
113          alloc->bits[i] = alloc->bits[i-1];
114       } else {
115          int j;
116          /* FIXME: We could save memory here */
117          alloc->bits[i] = celt_alloc(MAX_PULSES*sizeof(int));
118          for (j=0;j<MAX_PULSES;j++)
119          {
120             int done = 0;
121             int pulses = j;
122             /* For bands where there's no pitch, id 1 corresponds to intra prediction 
123                with no pulse. id 2 means intra prediction with one pulse, and so on.*/
124             if (eBands[i] >= m->pitchEnd)
125                pulses -= 1;
126             if (pulses < 0)
127                alloc->bits[i][j] = 0;
128             else {
129                alloc->bits[i][j] = log2_frac64(ncwrs64(N, pulses),BITRES);
130                /* FIXME: Could there be a better test for the max number of pulses that fit in 64 bits? */
131                if (alloc->bits[i][j] > (60<<BITRES))
132                   done = 1;
133                /* Add the intra-frame prediction bits */
134                if (eBands[i] >= m->pitchEnd)
135                {
136                   int max_pos = 2*eBands[i]-eBands[i+1];
137                   if (max_pos > 32)
138                      max_pos = 32;
139                   alloc->bits[i][j] += (1<<BITRES) + log2_frac(max_pos,BITRES);
140                }
141             }
142             if (done)
143                break;
144          }
145          for (;j<MAX_PULSES;j++)
146             alloc->bits[i][j] = BITOVERFLOW;
147          prevN = N;
148       }
149    }
150 }
151
152 void alloc_clear(struct alloc_data *alloc)
153 {
154    int i;
155    int *prevPtr = NULL;
156    for (i=0;i<alloc->len;i++)
157    {
158       if (alloc->bits[i] != prevPtr)
159       {
160          prevPtr = alloc->bits[i];
161          celt_free(alloc->bits[i]);
162       }
163    }
164    celt_free(alloc->bits);
165 }
166
167 int bits2pulses(const struct alloc_data *alloc, int band, int bits)
168 {
169    int lo, hi;
170    lo = 0;
171    hi = MAX_PULSES-1;
172    
173    while (hi-lo != 1)
174    {
175       int mid = (lo+hi)>>1;
176       if (alloc->bits[band][mid] >= bits)
177          hi = mid;
178       else
179          lo = mid;
180    }
181    if (bits-alloc->bits[band][lo] <= alloc->bits[band][hi]-bits)
182       return lo;
183    else
184       return hi;
185 }
186
187 int vec_bits2pulses(const struct alloc_data *alloc, const int *bands, int *bits, int *pulses, int len)
188 {
189    int i, BC;
190    int sum=0;
191    BC = alloc->mode->nbMdctBlocks*alloc->mode->nbChannels;
192
193    for (i=0;i<len;i++)
194    {
195       pulses[i] = bits2pulses(alloc, i, bits[i]);
196       sum += alloc->bits[i][pulses[i]];
197    }
198    //printf ("sum = %d\n", sum);
199    return sum;
200 }
201
202 int interp_bits2pulses(const struct alloc_data *alloc, int *bits1, int *bits2, int total, int *pulses, int len)
203 {
204    int lo, hi, out;
205    int j;
206    int bits[len];
207    const int *bands = alloc->bands;
208    lo = 0;
209    hi = 1<<BITRES;
210    while (hi-lo != 1)
211    {
212       int mid = (lo+hi)>>1;
213       for (j=0;j<len;j++)
214          bits[j] = ((1<<BITRES)-mid)*bits1[j] + mid*bits2[j];
215       if (vec_bits2pulses(alloc, bands, bits, pulses, len) > total<<BITRES)
216          hi = mid;
217       else
218          lo = mid;
219    }
220    //printf ("interp bisection gave %d\n", lo);
221    for (j=0;j<len;j++)
222       bits[j] = ((1<<BITRES)-lo)*bits1[j] + lo*bits2[j];
223    out = vec_bits2pulses(alloc, bands, bits, pulses, len);
224    /* Do some refinement to use up all bits. In the first pass, we can only add pulses to 
225       bands that are under their allocated budget. In the second pass, anything goes */
226    int firstpass = 1;
227    while(1)
228    {
229       int incremented = 0;
230       for (j=0;j<len;j++)
231       {
232          if ((!firstpass || alloc->bits[j][pulses[j]] < bits[j]) && pulses[j]<MAX_PULSES-1)
233          {
234             if (out+alloc->bits[j][pulses[j]+1]-alloc->bits[j][pulses[j]] <= total<<BITRES)
235             {
236                out = out+alloc->bits[j][pulses[j]+1]-alloc->bits[j][pulses[j]];
237                pulses[j] += 1;
238                incremented = 1;
239                //printf ("INCREMENT %d\n", j);
240             }
241          }
242       }
243       if (!incremented)
244       {
245          if (firstpass)
246             firstpass = 0;
247          else
248             break;
249       }
250    }
251    return (out+BITROUND) >> BITRES;
252 }
253
254 int compute_allocation(const struct alloc_data *alloc, int *offsets, int total, int *pulses)
255 {
256    int lo, hi, len;
257    const CELTMode *m;
258
259    m = alloc->mode;
260    len = m->nbEBands;
261    lo = 0;
262    hi = m->nbAllocVectors - 1;
263    while (hi-lo != 1)
264    {
265       int j;
266       int bits[len];
267       int pulses[len];
268       int mid = (lo+hi) >> 1;
269       for (j=0;j<len;j++)
270       {
271          bits[j] = (m->allocVectors[mid*len+j] + offsets[j])<<BITRES;
272          if (bits[j] < 0)
273             bits[j] = 0;
274          //printf ("%d ", bits[j]);
275       }
276       //printf ("\n");
277       if (vec_bits2pulses(alloc, alloc->bands, bits, pulses, len) > total<<BITRES)
278          hi = mid;
279       else
280          lo = mid;
281       //printf ("lo = %d, hi = %d\n", lo, hi);
282    }
283    {
284       int bits1[len];
285       int bits2[len];
286       int j;
287       for (j=0;j<len;j++)
288       {
289          bits1[j] = m->allocVectors[lo*len+j] + offsets[j];
290          bits2[j] = m->allocVectors[hi*len+j] + offsets[j];
291          if (bits1[j] < 0)
292             bits1[j] = 0;
293          if (bits2[j] < 0)
294             bits2[j] = 0;
295       }
296       return interp_bits2pulses(alloc, bits1, bits2, total, pulses, len);
297    }
298 }
299
300 #if 0
301 int main()
302 {
303    int i;
304    printf ("log(128) = %d\n", EC_ILOG(128));
305    for(i=1;i<2000000000;i+=1738)
306    {
307       printf ("%d %d\n", i, log2_frac(i, 10));
308    }
309    return 0;
310 }
311 #endif
312 #if 0
313 int main()
314 {
315    int i;
316    int offsets[18] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
317    int bits[18] = {10, 9, 9, 8, 8, 8, 8, 8, 8, 8, 9, 10, 8, 9, 10, 11, 6, 7};
318    int bits1[18] = {8, 7, 7, 6, 6, 6, 5, 4, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5};
319    int bits2[18] = {15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15};
320    int bank[20] = {0,  4,  8, 12, 16, 20, 24, 28, 32, 38, 44, 52, 62, 74, 90,112,142,182, 232,256};
321    int pulses[18];
322    struct alloc_data alloc;
323    
324    alloc_init(&alloc, celt_mode0);
325    int b;
326    //b = vec_bits2pulses(&alloc, bank, bits, pulses, 18);
327    //printf ("total: %d bits\n", b);
328    //for (i=0;i<18;i++)
329    //   printf ("%d ", pulses[i]);
330    //printf ("\n");
331    //b = interp_bits2pulses(&alloc, bits1, bits2, 162, pulses, 18);
332    b = compute_allocation(&alloc, offsets, 190, pulses);
333    printf ("total: %d bits\n", b);
334    for (i=0;i<18;i++)
335       printf ("%d ", pulses[i]);
336    printf ("\n");
337
338    alloc_clear(&alloc);
339    return 0;
340 }
341 #endif