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