Optimisation: got rid of about 10% of the 32-bit divisions by using ec_enc_uint
[opus.git] / libcelt / quant_bands.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 "quant_bands.h"
37 #include "laplace.h"
38 #include <math.h>
39 #include "os_support.h"
40 #include "arch.h"
41 #include "mathops.h"
42 #include "stack_alloc.h"
43
44 #ifdef FIXED_POINT
45 const celt_word16_t eMeans[24] = {11520, -2048, -3072, -640, 256, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
46 #else
47 const celt_word16_t eMeans[24] = {45.f, -8.f, -12.f, -2.5f, 1.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f};
48 #endif
49
50 /*const int frac[24] = {4, 3, 3, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2};*/
51 const int frac[24] = {8, 6, 5, 4, 3, 3, 3, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2};
52
53 #ifdef FIXED_POINT
54 static inline celt_ener_t dB2Amp(celt_ener_t dB)
55 {
56    celt_ener_t amp;
57    amp = PSHR32(celt_exp2(MULT16_16_Q14(21771,dB)),2)-QCONST16(.3f, 14);
58    if (amp < 0)
59       amp = 0;
60    return amp;
61 }
62
63 #define DBofTWO 24661
64 static inline celt_word16_t amp2dB(celt_ener_t amp)
65 {
66    /* equivalent to return 6.0207*log2(.3+amp) */
67    return ROUND16(MULT16_16(24661,celt_log2(ADD32(QCONST32(.3f,14),amp))),12);
68    /* return DB_SCALING*20*log10(.3+ENER_SCALING_1*amp); */
69 }
70 #else
71 static inline celt_ener_t dB2Amp(celt_ener_t dB)
72 {
73    celt_ener_t amp;
74    amp = pow(10, .05*dB)-.3;
75    if (amp < 0)
76       amp = 0;
77    return amp;
78 }
79 static inline celt_word16_t amp2dB(celt_ener_t amp)
80 {
81    return 20*log10(.3+amp);
82 }
83 #endif
84
85 /* We could use ec_enc_uint directly, but this saves some divisions */
86 static inline void enc_frac(ec_enc *enc, int val, int ft)
87 {
88    switch(ft)
89    {
90       case 1:
91          break;
92       case 2:
93          ec_enc_bits(enc, val, 1);
94          break;
95       case 4:
96          ec_enc_bits(enc, val, 2);
97          break;
98       case 8:
99          ec_enc_bits(enc, val, 3);
100          break;
101       default:
102          ec_enc_uint(enc, val, ft);
103          break;
104    }
105 }
106
107 /* We could use ec_enc_uint directly, but this saves some divisions */
108 static inline int dec_frac(ec_dec *dec, int ft)
109 {
110    switch(ft)
111    {
112       case 1:
113          return 0;
114          break;
115       case 2:
116          return ec_dec_bits(dec, 1);
117          break;
118       case 4:
119          return ec_dec_bits(dec, 2);
120          break;
121       case 8:
122          return ec_dec_bits(dec, 3);
123          break;
124       default:
125          return ec_dec_uint(dec, ft);
126          break;
127    }
128 }
129
130 static const celt_word16_t base_resolution = QCONST16(6.f,8);
131
132 static void quant_energy_mono(const CELTMode *m, celt_ener_t *eBands, celt_word16_t *oldEBands, unsigned budget, ec_enc *enc)
133 {
134    int i;
135    unsigned bits;
136    celt_word16_t prev = 0;
137    celt_word16_t coef = m->ePredCoef;
138    celt_word16_t beta;
139    VARDECL(celt_word16_t, error);
140    SAVE_STACK;
141    /* The .7 is a heuristic */
142    beta = MULT16_16_Q15(QCONST16(.7f,15),coef);
143    
144    ALLOC(error, m->nbEBands, celt_word16_t);
145    bits = ec_enc_tell(enc, 0);
146    for (i=0;i<m->nbEBands;i++)
147    {
148       int qi;
149       celt_word16_t q;   /* dB */
150       celt_word16_t x;   /* dB */
151       celt_word16_t f;   /* Q8 */
152       celt_word16_t mean = MULT16_16_Q15(Q15ONE-coef,eMeans[i]);
153       x = amp2dB(eBands[i]);
154       f = EXTRACT16(celt_div(SHL32(EXTEND32(x-mean-MULT16_16_Q15(coef,oldEBands[i])-prev),8),base_resolution));
155 #ifdef FIXED_POINT
156       /* Rounding to nearest integer here is really important! */
157       qi = (f+128)>>8;
158 #else
159       qi = (int)floor(.5+f);
160 #endif
161       /*ec_laplace_encode(enc, qi, i==0?11192:6192);*/
162       /*ec_laplace_encode(enc, qi, 8500-i*200);*/
163       /* If we don't have enough bits to encode all the energy, just assume something safe. */
164       if (ec_enc_tell(enc, 0) - bits > budget)
165          qi = -1;
166       else
167          ec_laplace_encode(enc, qi, 6000-i*200);
168       q = qi*base_resolution;
169       error[i] = f - SHL16(qi,8);
170       
171       /*printf("%d ", qi);*/
172       /*printf("%f %f ", pred+prev+q, x);*/
173       /*printf("%f ", x-pred);*/
174       
175       oldEBands[i] = mean+MULT16_16_Q15(coef,oldEBands[i])+prev+q;
176       
177       prev = mean+prev+MULT16_16_Q15(Q15ONE-beta,q);
178    }
179    /*bits = ec_enc_tell(enc, 0) - bits;*/
180    /*printf ("%d\n", bits);*/
181    for (i=0;i<m->nbEBands;i++)
182    {
183       int q2;
184       celt_word16_t offset = (error[i]+QCONST16(.5f,8))*frac[i];
185       /* FIXME: Instead of giving up without warning, we should degrade everything gracefully */
186       if (ec_enc_tell(enc, 0) - bits + celt_ilog2(frac[i]) >= budget)
187          break;
188 #ifdef FIXED_POINT
189       /* Has to be without rounding */
190       q2 = offset>>8;
191 #else
192       q2 = (int)floor(offset);
193 #endif
194       if (q2 > frac[i]-1)
195          q2 = frac[i]-1;
196       enc_frac(enc, q2, frac[i]);
197       offset = EXTRACT16(celt_div(SHL16(q2,8)+QCONST16(.5,8),frac[i])-QCONST16(.5f,8));
198       oldEBands[i] += PSHR32(MULT16_16(DB_SCALING*6,offset),8);
199       /*printf ("%f ", error[i] - offset);*/
200    }
201    for (i=0;i<m->nbEBands;i++)
202    {
203       eBands[i] = dB2Amp(oldEBands[i]);
204    }
205    /*printf ("%d\n", ec_enc_tell(enc, 0)-9);*/
206
207    /*printf ("\n");*/
208    RESTORE_STACK;
209 }
210
211 static void unquant_energy_mono(const CELTMode *m, celt_ener_t *eBands, celt_word16_t *oldEBands, unsigned budget, ec_dec *dec)
212 {
213    int i;
214    unsigned bits;
215    celt_word16_t prev = 0;
216    celt_word16_t coef = m->ePredCoef;
217    /* The .7 is a heuristic */
218    celt_word16_t beta = MULT16_16_Q15(QCONST16(.7f,15),coef);
219    bits = ec_dec_tell(dec, 0);
220    for (i=0;i<m->nbEBands;i++)
221    {
222       int qi;
223       celt_word16_t q;
224       celt_word16_t mean = MULT16_16_Q15(Q15ONE-coef,eMeans[i]);
225       /* If we didn't have enough bits to encode all the energy, just assume something safe. */
226       if (ec_dec_tell(dec, 0) - bits > budget)
227          qi = -1;
228       else
229          qi = ec_laplace_decode(dec, 6000-i*200);
230       q = qi*base_resolution;
231       
232       /*printf("%d ", qi);*/
233       /*printf("%f %f ", pred+prev+q, x);*/
234       /*printf("%f ", x-pred);*/
235       
236       oldEBands[i] = mean+MULT16_16_Q15(coef,oldEBands[i])+prev+q;
237       
238       prev = mean+prev+MULT16_16_Q15(Q15ONE-beta,q);
239    }
240    for (i=0;i<m->nbEBands;i++)
241    {
242       int q2;
243       celt_word16_t offset;
244       if (ec_dec_tell(dec, 0) - bits + celt_ilog2(frac[i]) >= budget)
245          break;
246       q2 = dec_frac(dec, frac[i]);
247       offset = EXTRACT16(celt_div(SHL16(q2,8)+QCONST16(.5,8),frac[i])-QCONST16(.5f,8));
248       oldEBands[i] += PSHR32(MULT16_16(DB_SCALING*6,offset),8);
249    }
250    for (i=0;i<m->nbEBands;i++)
251    {
252       eBands[i] = dB2Amp(oldEBands[i]);
253    }
254    /*printf ("\n");*/
255 }
256
257
258
259 void quant_energy(const CELTMode *m, celt_ener_t *eBands, celt_word16_t *oldEBands, int budget, ec_enc *enc)
260 {
261    int C;
262    SAVE_STACK;
263    
264    C = m->nbChannels;
265
266    if (C==1)
267       quant_energy_mono(m, eBands, oldEBands, budget, enc);
268    else 
269 #if 1
270    {
271       int c;
272       VARDECL(celt_ener_t, E);
273       ALLOC(E, m->nbEBands, celt_ener_t);
274       for (c=0;c<C;c++)
275       {
276          int i;
277          for (i=0;i<m->nbEBands;i++)
278             E[i] = eBands[C*i+c];
279          quant_energy_mono(m, E, oldEBands+c*m->nbEBands, budget/C, enc);
280          for (i=0;i<m->nbEBands;i++)
281             eBands[C*i+c] = E[i];
282       }
283    }
284 #else
285       if (C==2)
286    {
287       int i;
288       int NB = m->nbEBands;
289       celt_ener_t mid[NB];
290       celt_ener_t side[NB];
291       for (i=0;i<NB;i++)
292       {
293          //left = eBands[C*i];
294          //right = eBands[C*i+1];
295          mid[i] = ENER_SCALING_1*sqrt(eBands[C*i]*eBands[C*i] + eBands[C*i+1]*eBands[C*i+1]);
296          side[i] = 20*log10((ENER_SCALING_1*eBands[2*i]+.3)/(ENER_SCALING_1*eBands[2*i+1]+.3));
297          //printf ("%f %f ", mid[i], side[i]);
298       }
299       //printf ("\n");
300       quant_energy_mono(m, mid, oldEBands, enc);
301       for (i=0;i<NB;i++)
302          side[i] = pow(10.f,floor(.5f+side[i])/10.f);
303          
304       //quant_energy_side(m, side, oldEBands+NB, enc);
305       for (i=0;i<NB;i++)
306       {
307          eBands[C*i] = ENER_SCALING*mid[i]*sqrt(side[i]/(1.f+side[i]));
308          eBands[C*i+1] = ENER_SCALING*mid[i]*sqrt(1.f/(1.f+side[i]));
309          //printf ("%f %f ", mid[i], side[i]);
310       }
311
312    } else {
313       celt_fatal("more than 2 channels not supported");
314    }
315 #endif
316    RESTORE_STACK;
317 }
318
319
320
321 void unquant_energy(const CELTMode *m, celt_ener_t *eBands, celt_word16_t *oldEBands, int budget, ec_dec *dec)
322 {
323    int C;   
324    SAVE_STACK;
325    C = m->nbChannels;
326
327    if (C==1)
328       unquant_energy_mono(m, eBands, oldEBands, budget, dec);
329    else {
330       int c;
331       VARDECL(celt_ener_t, E);
332       ALLOC(E, m->nbEBands, celt_ener_t);
333       for (c=0;c<C;c++)
334       {
335          int i;
336          unquant_energy_mono(m, E, oldEBands+c*m->nbEBands, budget/C, dec);
337          for (i=0;i<m->nbEBands;i++)
338             eBands[C*i+c] = E[i];
339       }
340    }
341    RESTORE_STACK;
342 }