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