More C89 fixes, making sure to include config.h from all source files.
[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 #ifdef HAVE_CONFIG_H
33 #include "config.h"
34 #endif
35
36 #include <math.h>
37 #include "modes.h"
38 #include "cwrs.h"
39 #include "arch.h"
40 #include "os_support.h"
41
42 #include "entcode.h"
43 #include "rate.h"
44
45 #define BITRES 4
46 #define BITROUND 8
47 #define BITOVERFLOW 10000
48
49 #define MAX_PULSES 64
50
51 static int log2_frac(ec_uint32 val, int frac)
52 {
53    int i;
54    /* EC_ILOG() actually returns log2()+1, go figure */
55    int L = EC_ILOG(val)-1;
56    /*printf ("in: %d %d ", val, L);*/
57    if (L>14)
58       val >>= L-14;
59    else if (L<14)
60       val <<= 14-L;
61    L <<= frac;
62    /*printf ("%d\n", val);*/
63    for (i=0;i<frac;i++)
64    {
65       val = (val*val) >> 15;
66       /*printf ("%d\n", val);*/
67       if (val > 16384)
68          L |= (1<<(frac-i-1));
69       else   
70          val <<= 1;
71    }
72    return L;
73 }
74
75 static int log2_frac64(ec_uint64 val, int frac)
76 {
77    int i;
78    /* EC_ILOG64() actually returns log2()+1, go figure */
79    int L = EC_ILOG64(val)-1;
80    /*printf ("in: %d %d ", val, L);*/
81    if (L>14)
82       val >>= L-14;
83    else if (L<14)
84       val <<= 14-L;
85    L <<= frac;
86    /*printf ("%d\n", val);*/
87    for (i=0;i<frac;i++)
88    {
89       val = (val*val) >> 15;
90       /*printf ("%d\n", val);*/
91       if (val > 16384)
92          L |= (1<<(frac-i-1));
93       else   
94          val <<= 1;
95    }
96    return L;
97 }
98
99 void compute_alloc_cache(CELTMode *m)
100 {
101    int i, prevN, BC;
102    int **bits;
103    const int *eBands = m->eBands;
104
105    bits = celt_alloc(m->nbEBands*sizeof(int*));
106    
107    BC = m->nbMdctBlocks*m->nbChannels;
108    prevN = -1;
109    for (i=0;i<m->nbEBands;i++)
110    {
111       int N = BC*(eBands[i+1]-eBands[i]);
112       if (N == prevN && eBands[i] < m->pitchEnd)
113       {
114          bits[i] = bits[i-1];
115       } else {
116          int j;
117          /* FIXME: We could save memory here */
118          bits[i] = celt_alloc(MAX_PULSES*sizeof(int));
119          for (j=0;j<MAX_PULSES;j++)
120          {
121             int done = 0;
122             int pulses = j;
123             /* For bands where there's no pitch, id 1 corresponds to intra prediction 
124             with no pulse. id 2 means intra prediction with one pulse, and so on.*/
125             if (eBands[i] >= m->pitchEnd)
126                pulses -= 1;
127             if (pulses < 0)
128                bits[i][j] = 0;
129             else {
130                bits[i][j] = log2_frac64(ncwrs64(N, pulses),BITRES);
131                /* FIXME: Could there be a better test for the max number of pulses that fit in 64 bits? */
132                if (bits[i][j] > (60<<BITRES))
133                   done = 1;
134                /* Add the intra-frame prediction bits */
135                if (eBands[i] >= m->pitchEnd)
136                {
137                   int max_pos = 2*eBands[i]-eBands[i+1];
138                   if (max_pos > 32)
139                      max_pos = 32;
140                   bits[i][j] += (1<<BITRES) + log2_frac(max_pos,BITRES);
141                }
142             }
143             if (done)
144                break;
145          }
146          for (;j<MAX_PULSES;j++)
147             bits[i][j] = BITOVERFLOW;
148          prevN = N;
149       }
150    }
151    m->bits = (const int * const *)bits;
152 }
153
154
155 int bits2pulses(const CELTMode *m, int band, int bits)
156 {
157    int lo, hi;
158    lo = 0;
159    hi = MAX_PULSES-1;
160    
161    while (hi-lo != 1)
162    {
163       int mid = (lo+hi)>>1;
164       if (m->bits[band][mid] >= bits)
165          hi = mid;
166       else
167          lo = mid;
168    }
169    if (bits-m->bits[band][lo] <= m->bits[band][hi]-bits)
170       return lo;
171    else
172       return hi;
173 }
174
175 int vec_bits2pulses(const CELTMode *m, const int *bands, int *bits, int *pulses, int len)
176 {
177    int i, BC;
178    int sum=0;
179    BC = m->nbMdctBlocks*m->nbChannels;
180
181    for (i=0;i<len;i++)
182    {
183       pulses[i] = bits2pulses(m, i, bits[i]);
184       sum += m->bits[i][pulses[i]];
185    }
186    /*printf ("sum = %d\n", sum);*/
187    return sum;
188 }
189
190 int interp_bits2pulses(const CELTMode *m, int *bits1, int *bits2, int total, int *pulses, int len)
191 {
192    int lo, hi, out;
193    int j;
194    int bits[len];
195    const int *bands = m->eBands;
196    lo = 0;
197    hi = 1<<BITRES;
198    while (hi-lo != 1)
199    {
200       int mid = (lo+hi)>>1;
201       for (j=0;j<len;j++)
202          bits[j] = ((1<<BITRES)-mid)*bits1[j] + mid*bits2[j];
203       if (vec_bits2pulses(m, bands, bits, pulses, len) > total<<BITRES)
204          hi = mid;
205       else
206          lo = mid;
207    }
208    /*printf ("interp bisection gave %d\n", lo);*/
209    for (j=0;j<len;j++)
210       bits[j] = ((1<<BITRES)-lo)*bits1[j] + lo*bits2[j];
211    out = vec_bits2pulses(m, bands, bits, pulses, len);
212    /* Do some refinement to use up all bits. In the first pass, we can only add pulses to 
213       bands that are under their allocated budget. In the second pass, anything goes */
214    int firstpass = 1;
215    while(1)
216    {
217       int incremented = 0;
218       for (j=0;j<len;j++)
219       {
220          if ((!firstpass || m->bits[j][pulses[j]] < bits[j]) && pulses[j]<MAX_PULSES-1)
221          {
222             if (out+m->bits[j][pulses[j]+1]-m->bits[j][pulses[j]] <= total<<BITRES)
223             {
224                out = out+m->bits[j][pulses[j]+1]-m->bits[j][pulses[j]];
225                pulses[j] += 1;
226                incremented = 1;
227             }
228          }
229       }
230       if (!incremented)
231       {
232          if (firstpass)
233             firstpass = 0;
234          else
235             break;
236       }
237    }
238    return (out+BITROUND) >> BITRES;
239 }
240
241 int compute_allocation(const CELTMode *m, int *offsets, int total, int *pulses)
242 {
243    int lo, hi, len;
244
245    len = m->nbEBands;
246    lo = 0;
247    hi = m->nbAllocVectors - 1;
248    while (hi-lo != 1)
249    {
250       int j;
251       int bits[len];
252       int pulses[len];
253       int mid = (lo+hi) >> 1;
254       for (j=0;j<len;j++)
255       {
256          bits[j] = (m->allocVectors[mid*len+j] + offsets[j])<<BITRES;
257          if (bits[j] < 0)
258             bits[j] = 0;
259          /*printf ("%d ", bits[j]);*/
260       }
261       /*printf ("\n");*/
262       if (vec_bits2pulses(m, m->eBands, bits, pulses, len) > total<<BITRES)
263          hi = mid;
264       else
265          lo = mid;
266       /*printf ("lo = %d, hi = %d\n", lo, hi);*/
267    }
268    {
269       int bits1[len];
270       int bits2[len];
271       int j;
272       for (j=0;j<len;j++)
273       {
274          bits1[j] = m->allocVectors[lo*len+j] + offsets[j];
275          bits2[j] = m->allocVectors[hi*len+j] + offsets[j];
276          if (bits1[j] < 0)
277             bits1[j] = 0;
278          if (bits2[j] < 0)
279             bits2[j] = 0;
280       }
281       return interp_bits2pulses(m, bits1, bits2, total, pulses, len);
282    }
283 }
284
285 #if 0
286 int main()
287 {
288    int i;
289    printf ("log(128) = %d\n", EC_ILOG(128));
290    for(i=1;i<2000000000;i+=1738)
291    {
292       printf ("%d %d\n", i, log2_frac(i, 10));
293    }
294    return 0;
295 }
296 #endif
297 #if 0
298 int main()
299 {
300    int i;
301    int offsets[18] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
302    int bits[18] = {10, 9, 9, 8, 8, 8, 8, 8, 8, 8, 9, 10, 8, 9, 10, 11, 6, 7};
303    int bits1[18] = {8, 7, 7, 6, 6, 6, 5, 4, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5};
304    int bits2[18] = {15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15};
305    int bank[20] = {0,  4,  8, 12, 16, 20, 24, 28, 32, 38, 44, 52, 62, 74, 90,112,142,182, 232,256};
306    int pulses[18];
307    struct alloc_data alloc;
308    
309    alloc_init(&alloc, celt_mode0);
310    int b;
311    //b = vec_bits2pulses(&alloc, bank, bits, pulses, 18);
312    //printf ("total: %d bits\n", b);
313    //for (i=0;i<18;i++)
314    //   printf ("%d ", pulses[i]);
315    //printf ("\n");
316    //b = interp_bits2pulses(&alloc, bits1, bits2, 162, pulses, 18);
317    b = compute_allocation(&alloc, offsets, 190, pulses);
318    printf ("total: %d bits\n", b);
319    for (i=0;i<18;i++)
320       printf ("%d ", pulses[i]);
321    printf ("\n");
322
323    alloc_clear(&alloc);
324    return 0;
325 }
326 #endif