optimisation: changed some for() loops to do-while() to give the compiler
[opus.git] / libcelt / rate.c
1 /* (C) 2007-2008 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 #ifndef STATIC_MODES
50 static int log2_frac(ec_uint32 val, int frac)
51 {
52    int i;
53    /* EC_ILOG() actually returns log2()+1, go figure */
54    int L = EC_ILOG(val)-1;
55    /*printf ("in: %d %d ", val, L);*/
56    if (L>14)
57       val >>= L-14;
58    else if (L<14)
59       val <<= 14-L;
60    L <<= frac;
61    /*printf ("%d\n", val);*/
62    for (i=0;i<frac;i++)
63    {
64       val = (val*val) >> 15;
65       /*printf ("%d\n", val);*/
66       if (val > 16384)
67          L |= (1<<(frac-i-1));
68       else   
69          val <<= 1;
70    }
71    return L;
72 }
73
74 static int log2_frac64(ec_uint64 val, int frac)
75 {
76    int i;
77    /* EC_ILOG64() actually returns log2()+1, go figure */
78    int L = EC_ILOG64(val)-1;
79    /*printf ("in: %d %d ", val, L);*/
80    if (L>14)
81       val >>= L-14;
82    else if (L<14)
83       val <<= 14-L;
84    L <<= frac;
85    /*printf ("%d\n", val);*/
86    for (i=0;i<frac;i++)
87    {
88       val = (val*val) >> 15;
89       /*printf ("%d\n", val);*/
90       if (val > 16384)
91          L |= (1<<(frac-i-1));
92       else   
93          val <<= 1;
94    }
95    return L;
96 }
97
98 void compute_alloc_cache(CELTMode *m)
99 {
100    int i, prevN, BC;
101    celt_int16_t **bits;
102    const celt_int16_t *eBands = m->eBands;
103
104    bits = celt_alloc(m->nbEBands*sizeof(celt_int16_t*));
105    
106    BC = m->nbMdctBlocks*m->nbChannels;
107    prevN = -1;
108    for (i=0;i<m->nbEBands;i++)
109    {
110       int N = BC*(eBands[i+1]-eBands[i]);
111       if (N == prevN && eBands[i] < m->pitchEnd)
112       {
113          bits[i] = bits[i-1];
114       } else {
115          int j;
116          /* FIXME: We could save memory here */
117          bits[i] = celt_alloc(MAX_PULSES*sizeof(celt_int16_t));
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                bits[i][j] = 0;
128             else {
129                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 (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                   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             bits[i][j] = BITOVERFLOW;
147          prevN = N;
148       }
149    }
150    m->bits = (const celt_int16_t * const *)bits;
151 }
152
153 #endif /* !STATIC_MODES */
154
155 static int bits2pulses(const CELTMode *m, int band, int bits)
156 {
157    int i;
158    int lo, hi;
159    lo = 0;
160    hi = MAX_PULSES-1;
161    
162    /* Instead of using the "bisection confition" we use a fixed number of 
163       iterations because it should be faster */
164    /*while (hi-lo != 1)*/
165    for (i=0;i<LOG_MAX_PULSES;i++)
166    {
167       int mid = (lo+hi)>>1;
168       /* OPT: Make sure this is implemented with a conditional move */
169       if (m->bits[band][mid] >= bits)
170          hi = mid;
171       else
172          lo = mid;
173    }
174    if (bits-m->bits[band][lo] <= m->bits[band][hi]-bits)
175       return lo;
176    else
177       return hi;
178 }
179
180 static int vec_bits2pulses(const CELTMode *m, int *bits, int *pulses, int len)
181 {
182    int i;
183    int sum=0;
184
185    for (i=0;i<len;i++)
186    {
187       pulses[i] = bits2pulses(m, i, bits[i]);
188       sum += m->bits[i][pulses[i]];
189    }
190    /*printf ("sum = %d\n", sum);*/
191    return sum;
192 }
193
194 static int interp_bits2pulses(const CELTMode *m, int *bits1, int *bits2, int total, int *pulses, int len)
195 {
196    int lo, hi, out;
197    int j;
198    int firstpass;
199    VARDECL(int, bits);
200    SAVE_STACK;
201    ALLOC(bits, len, int);
202    lo = 0;
203    hi = 1<<BITRES;
204    while (hi-lo != 1)
205    {
206       int mid = (lo+hi)>>1;
207       for (j=0;j<len;j++)
208          bits[j] = ((1<<BITRES)-mid)*bits1[j] + mid*bits2[j];
209       if (vec_bits2pulses(m, bits, pulses, len) > total<<BITRES)
210          hi = mid;
211       else
212          lo = mid;
213    }
214    /*printf ("interp bisection gave %d\n", lo);*/
215    for (j=0;j<len;j++)
216       bits[j] = ((1<<BITRES)-lo)*bits1[j] + lo*bits2[j];
217    out = vec_bits2pulses(m, bits, pulses, len);
218    /* Do some refinement to use up all bits. In the first pass, we can only add pulses to 
219       bands that are under their allocated budget. In the second pass, anything goes */
220    firstpass = 1;
221    while(1)
222    {
223       int incremented = 0;
224       for (j=0;j<len;j++)
225       {
226          if ((!firstpass || m->bits[j][pulses[j]] < bits[j]) && pulses[j]<MAX_PULSES-1)
227          {
228             if (out+m->bits[j][pulses[j]+1]-m->bits[j][pulses[j]] <= total<<BITRES)
229             {
230                out = out+m->bits[j][pulses[j]+1]-m->bits[j][pulses[j]];
231                pulses[j] += 1;
232                incremented = 1;
233             }
234          }
235       }
236       if (!incremented)
237       {
238          if (firstpass)
239             firstpass = 0;
240          else
241             break;
242       }
243    }
244    RESTORE_STACK;
245    return (out+BITROUND) >> BITRES;
246 }
247
248 int compute_allocation(const CELTMode *m, int *offsets, int total, int *pulses)
249 {
250    int lo, hi, len, ret;
251    VARDECL(int, bits1);
252    VARDECL(int, bits2);
253    SAVE_STACK;
254    
255    len = m->nbEBands;
256    ALLOC(bits1, len, int);
257    ALLOC(bits2, len, int);
258    lo = 0;
259    hi = m->nbAllocVectors - 1;
260    while (hi-lo != 1)
261    {
262       int j;
263       int mid = (lo+hi) >> 1;
264       for (j=0;j<len;j++)
265       {
266          bits1[j] = (m->allocVectors[mid*len+j] + offsets[j])<<BITRES;
267          if (bits1[j] < 0)
268             bits1[j] = 0;
269          /*printf ("%d ", bits[j]);*/
270       }
271       /*printf ("\n");*/
272       if (vec_bits2pulses(m, bits1, pulses, len) > total<<BITRES)
273          hi = mid;
274       else
275          lo = mid;
276       /*printf ("lo = %d, hi = %d\n", lo, hi);*/
277    }
278    {
279       int j;
280       for (j=0;j<len;j++)
281       {
282          bits1[j] = m->allocVectors[lo*len+j] + offsets[j];
283          bits2[j] = m->allocVectors[hi*len+j] + offsets[j];
284          if (bits1[j] < 0)
285             bits1[j] = 0;
286          if (bits2[j] < 0)
287             bits2[j] = 0;
288       }
289       ret = interp_bits2pulses(m, bits1, bits2, total, pulses, len);
290       RESTORE_STACK;
291       return ret;
292    }
293 }
294