fixed-point: copied the exp2 implementation from Speex, using it for dB2Amp()
[opus.git] / libcelt / quant_bands.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 "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
43 #ifdef FIXED_POINT
44 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};
45 #else
46 const float 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};
47 #endif
48
49 /*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};*/
50 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};
51
52 #ifdef FIXED_POINT
53 static inline celt_ener_t dB2Amp(celt_ener_t dB)
54 {
55    celt_ener_t amp;
56    amp = PSHR32(celt_exp2(dB/6.0207*8),2)-QCONST16(.3f, 14);
57    if (amp < 0)
58       amp = 0;
59    return amp;
60 }
61
62 #define DBofTWO 24661
63 static inline celt_word16_t amp2dB(celt_ener_t amp)
64 {
65    /* equivalent to return 6.0207*log2(.3+amp) */
66    return ROUND(MULT16_16(24661,celt_log2(ADD32(QCONST32(.3f,14),amp))),12);
67    /* return DB_SCALING*20*log10(.3+ENER_SCALING_1*amp); */
68 }
69 #else
70 static inline celt_ener_t dB2Amp(celt_ener_t dB)
71 {
72    celt_ener_t amp;
73    amp = pow(10, .05*dB)-.3;
74    if (amp < 0)
75       amp = 0;
76    return amp;
77 }
78 static inline celt_word16_t amp2dB(celt_ener_t amp)
79 {
80    return 20*log10(.3+amp);
81 }
82 #endif
83
84
85 #ifdef FIXED_POINT
86 #define Q8 256.f
87 #define Q8_1 (1.f/256.f)
88 #else
89 #define Q8 1.f
90 #define Q8_1 1.f
91 #endif
92 static void quant_energy_mono(const CELTMode *m, celt_ener_t *eBands, celt_word16_t *oldEBands, int budget, ec_enc *enc)
93 {
94    int i;
95    int bits;
96    celt_word16_t prev = 0;
97    float coef = m->ePredCoef;
98    VARDECL(celt_word16_t *error);
99    /* The .7 is a heuristic */
100    float beta = .7*coef;
101    
102    ALLOC(error, m->nbEBands, celt_word16_t);
103    bits = ec_enc_tell(enc, 0);
104    for (i=0;i<m->nbEBands;i++)
105    {
106       int qi;
107       celt_word16_t q;   /* dB */
108       celt_word16_t res; /* dB */
109       celt_word16_t x;   /* dB */
110       celt_word16_t f;   /* Q8 */
111       celt_word16_t mean = (1-coef)*eMeans[i];
112       x = amp2dB(eBands[i]);
113       res = DB_SCALING*6;
114       f = QCONST16(1.f,8)*(x-mean-coef*oldEBands[i]-prev*1.f)/res;
115 #ifdef FIXED_POINT
116       /* Rounding to nearest integer here is really important! */
117       qi = (f+128)>>8;
118 #else
119       qi = (int)floor(.5+f);
120 #endif
121       /*ec_laplace_encode(enc, qi, i==0?11192:6192);*/
122       /*ec_laplace_encode(enc, qi, 8500-i*200);*/
123       /* If we don't have enough bits to encode all the energy, just assume something safe. */
124       if (ec_enc_tell(enc, 0) - bits > budget)
125          qi = -1;
126       else
127          ec_laplace_encode(enc, qi, 6000-i*200);
128       q = qi*res;
129       error[i] = f - SHL16(qi,8);
130       
131       /*printf("%d ", qi);*/
132       /*printf("%f %f ", pred+prev+q, x);*/
133       /*printf("%f ", x-pred);*/
134       
135       oldEBands[i] = mean+coef*oldEBands[i]+prev+q;
136       
137       prev = mean+prev+(1-beta)*q;
138    }
139    /*bits = ec_enc_tell(enc, 0) - bits;*/
140    /*printf ("%d\n", bits);*/
141    for (i=0;i<m->nbEBands;i++)
142    {
143       int q2;
144       celt_word16_t offset = (error[i]+QCONST16(.5f,8))*frac[i];
145       /* FIXME: Instead of giving up without warning, we should degrade everything gracefully */
146       if (ec_enc_tell(enc, 0) - bits +EC_ILOG(frac[i])> budget)
147          break;
148 #ifdef FIXED_POINT
149       /* Has to be without rounding */
150       q2 = offset>>8;
151 #else
152       q2 = (int)floor(Q8_1*offset);
153 #endif
154       if (q2 > frac[i]-1)
155          q2 = frac[i]-1;
156       ec_enc_uint(enc, q2, frac[i]);
157       offset = (Q8*(q2+.5)/frac[i])-QCONST16(.5f,8);
158       oldEBands[i] += PSHR32(MULT16_16(DB_SCALING*6,offset),8);
159       /*printf ("%f ", error[i] - offset);*/
160    }
161    for (i=0;i<m->nbEBands;i++)
162    {
163       eBands[i] = dB2Amp(oldEBands[i]);
164    }
165    /*printf ("%d\n", ec_enc_tell(enc, 0)-9);*/
166
167    /*printf ("\n");*/
168 }
169
170 static void unquant_energy_mono(const CELTMode *m, celt_ener_t *eBands, celt_word16_t *oldEBands, int budget, ec_dec *dec)
171 {
172    int i;
173    int bits;
174    celt_word16_t prev = 0;
175    float coef = m->ePredCoef;
176    /* The .7 is a heuristic */
177    float beta = .7*coef;
178    bits = ec_dec_tell(dec, 0);
179    for (i=0;i<m->nbEBands;i++)
180    {
181       int qi;
182       celt_word16_t q;
183       celt_word16_t res;
184       celt_word16_t mean = (1-coef)*eMeans[i];
185       res = DB_SCALING*6.;
186       /* If we didn't have enough bits to encode all the energy, just assume something safe. */
187       if (ec_dec_tell(dec, 0) - bits > budget)
188          qi = -1;
189       else
190          qi = ec_laplace_decode(dec, 6000-i*200);
191       q = qi*res;
192       
193       /*printf("%d ", qi);*/
194       /*printf("%f %f ", pred+prev+q, x);*/
195       /*printf("%f ", x-pred);*/
196       
197       oldEBands[i] = mean+coef*oldEBands[i]+prev+q;
198       
199       prev = mean+prev+(1-beta)*q;
200    }
201    for (i=0;i<m->nbEBands;i++)
202    {
203       int q2;
204       celt_word16_t offset;
205       if (ec_dec_tell(dec, 0) - bits +EC_ILOG(frac[i])> budget)
206          break;
207       q2 = ec_dec_uint(dec, frac[i]);
208       offset = (Q8*(q2+.5)/frac[i])-QCONST16(.5f,8);
209       oldEBands[i] += PSHR32(MULT16_16(DB_SCALING*6,offset),8);
210    }
211    for (i=0;i<m->nbEBands;i++)
212    {
213       eBands[i] = dB2Amp(oldEBands[i]);
214    }
215    /*printf ("\n");*/
216 }
217
218
219
220 void quant_energy(const CELTMode *m, celt_ener_t *eBands, celt_word16_t *oldEBands, int budget, ec_enc *enc)
221 {
222    int C;
223    
224    C = m->nbChannels;
225
226    if (C==1)
227       quant_energy_mono(m, eBands, oldEBands, budget, enc);
228    else 
229 #if 1
230    {
231       int c;
232       VARDECL(celt_ener_t *E);
233       ALLOC(E, m->nbEBands, celt_ener_t);
234       for (c=0;c<C;c++)
235       {
236          int i;
237          for (i=0;i<m->nbEBands;i++)
238             E[i] = eBands[C*i+c];
239          quant_energy_mono(m, E, oldEBands+c*m->nbEBands, budget/C, enc);
240          for (i=0;i<m->nbEBands;i++)
241             eBands[C*i+c] = E[i];
242       }
243    }
244 #else
245       if (C==2)
246    {
247       int i;
248       int NB = m->nbEBands;
249       float mid[NB];
250       float side[NB];
251       float left;
252       float right;
253       for (i=0;i<NB;i++)
254       {
255          //left = eBands[C*i];
256          //right = eBands[C*i+1];
257          mid[i] = ENER_SCALING_1*sqrt(eBands[C*i]*eBands[C*i] + eBands[C*i+1]*eBands[C*i+1]);
258          side[i] = 20*log10((ENER_SCALING_1*eBands[2*i]+.3)/(ENER_SCALING_1*eBands[2*i+1]+.3));
259          //printf ("%f %f ", mid[i], side[i]);
260       }
261       //printf ("\n");
262       quant_energy_mono(m, mid, oldEBands, enc);
263       for (i=0;i<NB;i++)
264          side[i] = pow(10.f,floor(.5f+side[i])/10.f);
265          
266       //quant_energy_side(m, side, oldEBands+NB, enc);
267       for (i=0;i<NB;i++)
268       {
269          eBands[C*i] = ENER_SCALING*mid[i]*sqrt(side[i]/(1.f+side[i]));
270          eBands[C*i+1] = ENER_SCALING*mid[i]*sqrt(1.f/(1.f+side[i]));
271          //printf ("%f %f ", mid[i], side[i]);
272       }
273
274    } else {
275       celt_fatal("more than 2 channels not supported");
276    }
277 #endif
278 }
279
280
281
282 void unquant_energy(const CELTMode *m, celt_ener_t *eBands, celt_word16_t *oldEBands, int budget, ec_dec *dec)
283 {
284    int C;   
285    C = m->nbChannels;
286
287    if (C==1)
288       unquant_energy_mono(m, eBands, oldEBands, budget, dec);
289    else {
290       int c;
291       VARDECL(celt_ener_t *E);
292       ALLOC(E, m->nbEBands, celt_ener_t);
293       for (c=0;c<C;c++)
294       {
295          int i;
296          unquant_energy_mono(m, E, oldEBands+c*m->nbEBands, budget/C, dec);
297          for (i=0;i<m->nbEBands;i++)
298             eBands[C*i+c] = E[i];
299       }
300    }
301 }