Just commenting -- nothing to see.
[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    /* Encode at a fixed coarse resolution */
147    for (i=0;i<m->nbEBands;i++)
148    {
149       int qi;
150       celt_word16_t q;   /* dB */
151       celt_word16_t x;   /* dB */
152       celt_word16_t f;   /* Q8 */
153       celt_word16_t mean = MULT16_16_Q15(Q15ONE-coef,eMeans[i]);
154       x = amp2dB(eBands[i]);
155       f = EXTRACT16(celt_div(SHL32(EXTEND32(x-mean-MULT16_16_Q15(coef,oldEBands[i])-prev),8),base_resolution));
156 #ifdef FIXED_POINT
157       /* Rounding to nearest integer here is really important! */
158       qi = (f+128)>>8;
159 #else
160       qi = (int)floor(.5+f);
161 #endif
162       /* If we don't have enough bits to encode all the energy, just assume something safe. */
163       if (ec_enc_tell(enc, 0) - bits > budget)
164          qi = -1;
165       else
166          ec_laplace_encode(enc, qi, 6000-i*200);
167       q = qi*base_resolution;
168       error[i] = f - SHL16(qi,8);
169       
170       oldEBands[i] = mean+MULT16_16_Q15(coef,oldEBands[i])+prev+q;
171       
172       prev = mean+prev+MULT16_16_Q15(Q15ONE-beta,q);
173    }
174    /* Encode finer resolution */
175    for (i=0;i<m->nbEBands;i++)
176    {
177       int q2;
178       celt_word16_t offset = (error[i]+QCONST16(.5f,8))*frac[i];
179       /* FIXME: Instead of giving up without warning, we should degrade everything gracefully */
180       if (ec_enc_tell(enc, 0) - bits + celt_ilog2(frac[i]) >= budget)
181          break;
182 #ifdef FIXED_POINT
183       /* Has to be without rounding */
184       q2 = offset>>8;
185 #else
186       q2 = (int)floor(offset);
187 #endif
188       if (q2 > frac[i]-1)
189          q2 = frac[i]-1;
190       enc_frac(enc, q2, frac[i]);
191       offset = EXTRACT16(celt_div(SHL16(q2,8)+QCONST16(.5,8),frac[i])-QCONST16(.5f,8));
192       oldEBands[i] += PSHR32(MULT16_16(DB_SCALING*6,offset),8);
193       /*printf ("%f ", error[i] - offset);*/
194    }
195    for (i=0;i<m->nbEBands;i++)
196    {
197       eBands[i] = dB2Amp(oldEBands[i]);
198    }
199    /*printf ("%d\n", ec_enc_tell(enc, 0)-9);*/
200
201    /*printf ("\n");*/
202    RESTORE_STACK;
203 }
204
205 static void unquant_energy_mono(const CELTMode *m, celt_ener_t *eBands, celt_word16_t *oldEBands, unsigned budget, ec_dec *dec)
206 {
207    int i;
208    unsigned bits;
209    celt_word16_t prev = 0;
210    celt_word16_t coef = m->ePredCoef;
211    /* The .7 is a heuristic */
212    celt_word16_t beta = MULT16_16_Q15(QCONST16(.7f,15),coef);
213    bits = ec_dec_tell(dec, 0);
214    
215    /* Decode at a fixed coarse resolution */
216    for (i=0;i<m->nbEBands;i++)
217    {
218       int qi;
219       celt_word16_t q;
220       celt_word16_t mean = MULT16_16_Q15(Q15ONE-coef,eMeans[i]);
221       /* If we didn't have enough bits to encode all the energy, just assume something safe. */
222       if (ec_dec_tell(dec, 0) - bits > budget)
223          qi = -1;
224       else
225          qi = ec_laplace_decode(dec, 6000-i*200);
226       q = qi*base_resolution;
227       
228       oldEBands[i] = mean+MULT16_16_Q15(coef,oldEBands[i])+prev+q;
229       
230       prev = mean+prev+MULT16_16_Q15(Q15ONE-beta,q);
231    }
232    /* Decode finer resolution */
233    for (i=0;i<m->nbEBands;i++)
234    {
235       int q2;
236       celt_word16_t offset;
237       if (ec_dec_tell(dec, 0) - bits + celt_ilog2(frac[i]) >= budget)
238          break;
239       q2 = dec_frac(dec, frac[i]);
240       offset = EXTRACT16(celt_div(SHL16(q2,8)+QCONST16(.5,8),frac[i])-QCONST16(.5f,8));
241       oldEBands[i] += PSHR32(MULT16_16(DB_SCALING*6,offset),8);
242    }
243    for (i=0;i<m->nbEBands;i++)
244    {
245       eBands[i] = dB2Amp(oldEBands[i]);
246    }
247    /*printf ("\n");*/
248 }
249
250
251
252 void quant_energy(const CELTMode *m, celt_ener_t *eBands, celt_word16_t *oldEBands, int budget, ec_enc *enc)
253 {
254    int C;
255    SAVE_STACK;
256    
257    C = m->nbChannels;
258
259    if (C==1)
260       quant_energy_mono(m, eBands, oldEBands, budget, enc);
261    else 
262 #if 1
263    {
264       int c;
265       VARDECL(celt_ener_t, E);
266       ALLOC(E, m->nbEBands, celt_ener_t);
267       for (c=0;c<C;c++)
268       {
269          int i;
270          for (i=0;i<m->nbEBands;i++)
271             E[i] = eBands[C*i+c];
272          quant_energy_mono(m, E, oldEBands+c*m->nbEBands, budget/C, enc);
273          for (i=0;i<m->nbEBands;i++)
274             eBands[C*i+c] = E[i];
275       }
276    }
277 #else
278       if (C==2)
279    {
280       int i;
281       int NB = m->nbEBands;
282       celt_ener_t mid[NB];
283       celt_ener_t side[NB];
284       for (i=0;i<NB;i++)
285       {
286          //left = eBands[C*i];
287          //right = eBands[C*i+1];
288          mid[i] = ENER_SCALING_1*sqrt(eBands[C*i]*eBands[C*i] + eBands[C*i+1]*eBands[C*i+1]);
289          side[i] = 20*log10((ENER_SCALING_1*eBands[2*i]+.3)/(ENER_SCALING_1*eBands[2*i+1]+.3));
290          //printf ("%f %f ", mid[i], side[i]);
291       }
292       //printf ("\n");
293       quant_energy_mono(m, mid, oldEBands, enc);
294       for (i=0;i<NB;i++)
295          side[i] = pow(10.f,floor(.5f+side[i])/10.f);
296          
297       //quant_energy_side(m, side, oldEBands+NB, enc);
298       for (i=0;i<NB;i++)
299       {
300          eBands[C*i] = ENER_SCALING*mid[i]*sqrt(side[i]/(1.f+side[i]));
301          eBands[C*i+1] = ENER_SCALING*mid[i]*sqrt(1.f/(1.f+side[i]));
302          //printf ("%f %f ", mid[i], side[i]);
303       }
304
305    } else {
306       celt_fatal("more than 2 channels not supported");
307    }
308 #endif
309    RESTORE_STACK;
310 }
311
312
313
314 void unquant_energy(const CELTMode *m, celt_ener_t *eBands, celt_word16_t *oldEBands, int budget, ec_dec *dec)
315 {
316    int C;   
317    SAVE_STACK;
318    C = m->nbChannels;
319
320    if (C==1)
321       unquant_energy_mono(m, eBands, oldEBands, budget, dec);
322    else {
323       int c;
324       VARDECL(celt_ener_t, E);
325       ALLOC(E, m->nbEBands, celt_ener_t);
326       for (c=0;c<C;c++)
327       {
328          int i;
329          unquant_energy_mono(m, E, oldEBands+c*m->nbEBands, budget/C, dec);
330          for (i=0;i<m->nbEBands;i++)
331             eBands[C*i+c] = E[i];
332       }
333    }
334    RESTORE_STACK;
335 }