Trying to clean up celt_ilog2() vs. EC_ILOG a bit.
[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 static const celt_word16_t base_resolution = QCONST16(6.f,8);
86
87 static void quant_energy_mono(const CELTMode *m, celt_ener_t *eBands, celt_word16_t *oldEBands, unsigned budget, ec_enc *enc)
88 {
89    int i;
90    unsigned bits;
91    celt_word16_t prev = 0;
92    celt_word16_t coef = m->ePredCoef;
93    celt_word16_t beta;
94    VARDECL(celt_word16_t, error);
95    SAVE_STACK;
96    /* The .7 is a heuristic */
97    beta = MULT16_16_Q15(QCONST16(.7f,15),coef);
98    
99    ALLOC(error, m->nbEBands, celt_word16_t);
100    bits = ec_enc_tell(enc, 0);
101    for (i=0;i<m->nbEBands;i++)
102    {
103       int qi;
104       celt_word16_t q;   /* dB */
105       celt_word16_t x;   /* dB */
106       celt_word16_t f;   /* Q8 */
107       celt_word16_t mean = MULT16_16_Q15(Q15ONE-coef,eMeans[i]);
108       x = amp2dB(eBands[i]);
109       f = EXTRACT16(celt_div(SHL32(EXTEND32(x-mean-MULT16_16_Q15(coef,oldEBands[i])-prev),8),base_resolution));
110 #ifdef FIXED_POINT
111       /* Rounding to nearest integer here is really important! */
112       qi = (f+128)>>8;
113 #else
114       qi = (int)floor(.5+f);
115 #endif
116       /*ec_laplace_encode(enc, qi, i==0?11192:6192);*/
117       /*ec_laplace_encode(enc, qi, 8500-i*200);*/
118       /* If we don't have enough bits to encode all the energy, just assume something safe. */
119       if (ec_enc_tell(enc, 0) - bits > budget)
120          qi = -1;
121       else
122          ec_laplace_encode(enc, qi, 6000-i*200);
123       q = qi*base_resolution;
124       error[i] = f - SHL16(qi,8);
125       
126       /*printf("%d ", qi);*/
127       /*printf("%f %f ", pred+prev+q, x);*/
128       /*printf("%f ", x-pred);*/
129       
130       oldEBands[i] = mean+MULT16_16_Q15(coef,oldEBands[i])+prev+q;
131       
132       prev = mean+prev+MULT16_16_Q15(Q15ONE-beta,q);
133    }
134    /*bits = ec_enc_tell(enc, 0) - bits;*/
135    /*printf ("%d\n", bits);*/
136    for (i=0;i<m->nbEBands;i++)
137    {
138       int q2;
139       celt_word16_t offset = (error[i]+QCONST16(.5f,8))*frac[i];
140       /* FIXME: Instead of giving up without warning, we should degrade everything gracefully */
141       if (ec_enc_tell(enc, 0) - bits + celt_ilog2(frac[i]) >= budget)
142          break;
143 #ifdef FIXED_POINT
144       /* Has to be without rounding */
145       q2 = offset>>8;
146 #else
147       q2 = (int)floor(offset);
148 #endif
149       if (q2 > frac[i]-1)
150          q2 = frac[i]-1;
151       ec_enc_uint(enc, q2, frac[i]);
152       offset = EXTRACT16(celt_div(SHL16(q2,8)+QCONST16(.5,8),frac[i])-QCONST16(.5f,8));
153       oldEBands[i] += PSHR32(MULT16_16(DB_SCALING*6,offset),8);
154       /*printf ("%f ", error[i] - offset);*/
155    }
156    for (i=0;i<m->nbEBands;i++)
157    {
158       eBands[i] = dB2Amp(oldEBands[i]);
159    }
160    /*printf ("%d\n", ec_enc_tell(enc, 0)-9);*/
161
162    /*printf ("\n");*/
163    RESTORE_STACK;
164 }
165
166 static void unquant_energy_mono(const CELTMode *m, celt_ener_t *eBands, celt_word16_t *oldEBands, unsigned budget, ec_dec *dec)
167 {
168    int i;
169    unsigned bits;
170    celt_word16_t prev = 0;
171    celt_word16_t coef = m->ePredCoef;
172    /* The .7 is a heuristic */
173    celt_word16_t beta = MULT16_16_Q15(QCONST16(.7f,15),coef);
174    bits = ec_dec_tell(dec, 0);
175    for (i=0;i<m->nbEBands;i++)
176    {
177       int qi;
178       celt_word16_t q;
179       celt_word16_t mean = MULT16_16_Q15(Q15ONE-coef,eMeans[i]);
180       /* If we didn't have enough bits to encode all the energy, just assume something safe. */
181       if (ec_dec_tell(dec, 0) - bits > budget)
182          qi = -1;
183       else
184          qi = ec_laplace_decode(dec, 6000-i*200);
185       q = qi*base_resolution;
186       
187       /*printf("%d ", qi);*/
188       /*printf("%f %f ", pred+prev+q, x);*/
189       /*printf("%f ", x-pred);*/
190       
191       oldEBands[i] = mean+MULT16_16_Q15(coef,oldEBands[i])+prev+q;
192       
193       prev = mean+prev+MULT16_16_Q15(Q15ONE-beta,q);
194    }
195    for (i=0;i<m->nbEBands;i++)
196    {
197       int q2;
198       celt_word16_t offset;
199       if (ec_dec_tell(dec, 0) - bits + celt_ilog2(frac[i]) >= budget)
200          break;
201       q2 = ec_dec_uint(dec, frac[i]);
202       offset = EXTRACT16(celt_div(SHL16(q2,8)+QCONST16(.5,8),frac[i])-QCONST16(.5f,8));
203       oldEBands[i] += PSHR32(MULT16_16(DB_SCALING*6,offset),8);
204    }
205    for (i=0;i<m->nbEBands;i++)
206    {
207       eBands[i] = dB2Amp(oldEBands[i]);
208    }
209    /*printf ("\n");*/
210 }
211
212
213
214 void quant_energy(const CELTMode *m, celt_ener_t *eBands, celt_word16_t *oldEBands, int budget, ec_enc *enc)
215 {
216    int C;
217    SAVE_STACK;
218    
219    C = m->nbChannels;
220
221    if (C==1)
222       quant_energy_mono(m, eBands, oldEBands, budget, enc);
223    else 
224 #if 1
225    {
226       int c;
227       VARDECL(celt_ener_t, E);
228       ALLOC(E, m->nbEBands, celt_ener_t);
229       for (c=0;c<C;c++)
230       {
231          int i;
232          for (i=0;i<m->nbEBands;i++)
233             E[i] = eBands[C*i+c];
234          quant_energy_mono(m, E, oldEBands+c*m->nbEBands, budget/C, enc);
235          for (i=0;i<m->nbEBands;i++)
236             eBands[C*i+c] = E[i];
237       }
238    }
239 #else
240       if (C==2)
241    {
242       int i;
243       int NB = m->nbEBands;
244       celt_ener_t mid[NB];
245       celt_ener_t side[NB];
246       for (i=0;i<NB;i++)
247       {
248          //left = eBands[C*i];
249          //right = eBands[C*i+1];
250          mid[i] = ENER_SCALING_1*sqrt(eBands[C*i]*eBands[C*i] + eBands[C*i+1]*eBands[C*i+1]);
251          side[i] = 20*log10((ENER_SCALING_1*eBands[2*i]+.3)/(ENER_SCALING_1*eBands[2*i+1]+.3));
252          //printf ("%f %f ", mid[i], side[i]);
253       }
254       //printf ("\n");
255       quant_energy_mono(m, mid, oldEBands, enc);
256       for (i=0;i<NB;i++)
257          side[i] = pow(10.f,floor(.5f+side[i])/10.f);
258          
259       //quant_energy_side(m, side, oldEBands+NB, enc);
260       for (i=0;i<NB;i++)
261       {
262          eBands[C*i] = ENER_SCALING*mid[i]*sqrt(side[i]/(1.f+side[i]));
263          eBands[C*i+1] = ENER_SCALING*mid[i]*sqrt(1.f/(1.f+side[i]));
264          //printf ("%f %f ", mid[i], side[i]);
265       }
266
267    } else {
268       celt_fatal("more than 2 channels not supported");
269    }
270 #endif
271    RESTORE_STACK;
272 }
273
274
275
276 void unquant_energy(const CELTMode *m, celt_ener_t *eBands, celt_word16_t *oldEBands, int budget, ec_dec *dec)
277 {
278    int C;   
279    SAVE_STACK;
280    C = m->nbChannels;
281
282    if (C==1)
283       unquant_energy_mono(m, eBands, oldEBands, budget, dec);
284    else {
285       int c;
286       VARDECL(celt_ener_t, E);
287       ALLOC(E, m->nbEBands, celt_ener_t);
288       for (c=0;c<C;c++)
289       {
290          int i;
291          unquant_energy_mono(m, E, oldEBands+c*m->nbEBands, budget/C, dec);
292          for (i=0;i<m->nbEBands;i++)
293             eBands[C*i+c] = E[i];
294       }
295    }
296    RESTORE_STACK;
297 }