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