Completed the separation of coarse and fine energy quantisation
[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 void compute_fine_allocation(const CELTMode *m, celt_int16_t *bits, int budget)
54 {
55    int i,j;
56    int len;
57    len = m->nbEBands;
58    for (i=0;i<m->nbAllocVectors;i++)
59    {
60       if (m->energy_alloc[i*(len+1)+len] > budget)
61          break;
62    }
63    if (i==0)
64    {
65       for (j=0;j<len;j++)
66          bits[j] = 0;
67    } else {
68       for (j=0;j<len;j++)
69          bits[j] = m->energy_alloc[(i-1)*(len+1)+j];
70       budget -= m->energy_alloc[(i-1)*(len+1)+len];
71    }
72    if (i<m->nbAllocVectors)
73    {
74       j=0;
75       while (budget>0)
76       {
77          if (m->energy_alloc[i*(len+1)+j]>bits[j])
78          {
79             bits[j]++;
80             budget--;
81          }
82          j++;
83          if (j>=len)
84             j=0;
85       }
86    }
87    
88    /*for (j=0;j<len;j++)
89       printf ("%d ", bits[j]);
90    printf ("\n");*/
91 }
92
93 #ifdef FIXED_POINT
94 static inline celt_ener_t dB2Amp(celt_ener_t dB)
95 {
96    celt_ener_t amp;
97    amp = PSHR32(celt_exp2(MULT16_16_Q14(21771,dB)),2)-QCONST16(.3f, 14);
98    if (amp < 0)
99       amp = 0;
100    return amp;
101 }
102
103 #define DBofTWO 24661
104 static inline celt_word16_t amp2dB(celt_ener_t amp)
105 {
106    /* equivalent to return 6.0207*log2(.3+amp) */
107    return ROUND16(MULT16_16(24661,celt_log2(ADD32(QCONST32(.3f,14),amp))),12);
108    /* return DB_SCALING*20*log10(.3+ENER_SCALING_1*amp); */
109 }
110 #else
111 static inline celt_ener_t dB2Amp(celt_ener_t dB)
112 {
113    celt_ener_t amp;
114    amp = pow(10, .05*dB)-.3;
115    if (amp < 0)
116       amp = 0;
117    return amp;
118 }
119 static inline celt_word16_t amp2dB(celt_ener_t amp)
120 {
121    return 20*log10(.3+amp);
122 }
123 #endif
124
125 static const celt_word16_t base_resolution = QCONST16(6.f,8);
126
127 int *quant_prob_alloc(const CELTMode *m)
128 {
129    int i;
130    int *prob;
131    prob = celt_alloc(2*m->nbEBands*sizeof(int));
132    for (i=0;i<m->nbEBands;i++)
133    {
134       prob[2*i] = 6000-i*200;
135       prob[2*i+1] = ec_laplace_get_start_freq(prob[2*i]);
136    }
137    return prob;
138 }
139
140 void quant_prob_free(int *freq)
141 {
142    celt_free(freq);
143 }
144
145 static void quant_coarse_energy_mono(const CELTMode *m, celt_ener_t *eBands, celt_word16_t *oldEBands, unsigned budget, int *prob, celt_word16_t *error, ec_enc *enc)
146 {
147    int i;
148    unsigned bits;
149    celt_word16_t prev = 0;
150    celt_word16_t coef = m->ePredCoef;
151    celt_word16_t beta;
152    /* The .7 is a heuristic */
153    beta = MULT16_16_Q15(QCONST16(.8f,15),coef);
154    
155    bits = ec_enc_tell(enc, 0);
156    /* Encode at a fixed coarse resolution */
157    for (i=0;i<m->nbEBands;i++)
158    {
159       int qi;
160       celt_word16_t q;   /* dB */
161       celt_word16_t x;   /* dB */
162       celt_word16_t f;   /* Q8 */
163       celt_word16_t mean = MULT16_16_Q15(Q15ONE-coef,eMeans[i]);
164       x = amp2dB(eBands[i]);
165       f = EXTRACT16(celt_div(SHL32(EXTEND32(x-mean-MULT16_16_Q15(coef,oldEBands[i])-prev),8),base_resolution));
166 #ifdef FIXED_POINT
167       /* Rounding to nearest integer here is really important! */
168       qi = (f+128)>>8;
169 #else
170       qi = (int)floor(.5+f);
171 #endif
172       /* If we don't have enough bits to encode all the energy, just assume something safe.
173          We allow slightly busting the budget here */
174       if (ec_enc_tell(enc, 0) - bits > budget+16)
175          qi = -1;
176       else
177          ec_laplace_encode_start(enc, &qi, prob[2*i], prob[2*i+1]);
178       q = qi*base_resolution;
179       error[i] = f - SHL16(qi,8);
180       
181       oldEBands[i] = mean+MULT16_16_Q15(coef,oldEBands[i])+prev+q;
182       if (oldEBands[i] < -QCONST16(12.f,8))
183          oldEBands[i] = -QCONST16(12.f,8);
184       prev = mean+prev+MULT16_16_Q15(Q15ONE-beta,q);
185    }
186 }
187
188 static void quant_fine_energy_mono(const CELTMode *m, celt_ener_t *eBands, celt_word16_t *oldEBands, celt_word16_t *error, celt_int16_t *fine_quant, ec_enc *enc)
189 {
190    int i;
191    /* Encode finer resolution */
192    for (i=0;i<m->nbEBands;i++)
193    {
194       int q2;
195       celt_int16_t frac = 1<<fine_quant[i];
196       celt_word16_t offset = (error[i]+QCONST16(.5f,8))*frac;
197       if (fine_quant[i] <= 0)
198          continue;
199 #ifdef FIXED_POINT
200       /* Has to be without rounding */
201       q2 = offset>>8;
202 #else
203       q2 = (int)floor(offset);
204 #endif
205       if (q2 > frac-1)
206          q2 = frac-1;
207       ec_enc_bits(enc, q2, fine_quant[i]);
208       offset = EXTRACT16(celt_div(SHL16(q2,8)+QCONST16(.5,8),frac)-QCONST16(.5f,8));
209       oldEBands[i] += PSHR32(MULT16_16(DB_SCALING*6,offset),8);
210       /*printf ("%f ", error[i] - offset);*/
211    }
212    for (i=0;i<m->nbEBands;i++)
213    {
214       eBands[i] = dB2Amp(oldEBands[i]);
215    }
216    /*printf ("%d\n", ec_enc_tell(enc, 0)-9);*/
217
218    /*printf ("\n");*/
219 }
220
221 static void unquant_coarse_energy_mono(const CELTMode *m, celt_ener_t *eBands, celt_word16_t *oldEBands, unsigned budget, int *prob, ec_dec *dec)
222 {
223    int i;
224    unsigned bits;
225    celt_word16_t prev = 0;
226    celt_word16_t coef = m->ePredCoef;
227    /* The .7 is a heuristic */
228    celt_word16_t beta = MULT16_16_Q15(QCONST16(.8f,15),coef);
229    
230    bits = ec_dec_tell(dec, 0);
231    /* Decode at a fixed coarse resolution */
232    for (i=0;i<m->nbEBands;i++)
233    {
234       int qi;
235       celt_word16_t q;
236       celt_word16_t mean = MULT16_16_Q15(Q15ONE-coef,eMeans[i]);
237       /* If we didn't have enough bits to encode all the energy, just assume something safe.
238          We allow slightly busting the budget here */
239       if (ec_dec_tell(dec, 0) - bits > budget+16)
240          qi = -1;
241       else
242          qi = ec_laplace_decode_start(dec, prob[2*i], prob[2*i+1]);
243       q = qi*base_resolution;
244       
245       oldEBands[i] = mean+MULT16_16_Q15(coef,oldEBands[i])+prev+q;
246       if (oldEBands[i] < -QCONST16(12.f,8))
247          oldEBands[i] = -QCONST16(12.f,8);
248       
249       prev = mean+prev+MULT16_16_Q15(Q15ONE-beta,q);
250    }
251 }
252
253 static void unquant_fine_energy_mono(const CELTMode *m, celt_ener_t *eBands, celt_word16_t *oldEBands, celt_int16_t *fine_quant, ec_dec *dec)
254 {
255    int i;
256    /* Decode finer resolution */
257    for (i=0;i<m->nbEBands;i++)
258    {
259       int q2;
260       celt_int16_t frac = 1<<fine_quant[i];
261       celt_word16_t offset;
262       if (fine_quant[i] <= 0)
263          continue;
264       q2 = ec_dec_bits(dec, fine_quant[i]);
265       offset = EXTRACT16(celt_div(SHL16(q2,8)+QCONST16(.5,8),frac)-QCONST16(.5f,8));
266       oldEBands[i] += PSHR32(MULT16_16(DB_SCALING*6,offset),8);
267    }
268    for (i=0;i<m->nbEBands;i++)
269    {
270       eBands[i] = dB2Amp(oldEBands[i]);
271    }
272    /*printf ("\n");*/
273 }
274
275
276
277 void quant_coarse_energy(const CELTMode *m, celt_ener_t *eBands, celt_word16_t *oldEBands, int budget, int *prob, celt_word16_t *error, ec_enc *enc)
278 {
279    int C;
280    C = m->nbChannels;
281
282    if (C==1)
283    {
284       quant_coarse_energy_mono(m, eBands, oldEBands, budget, prob, error, enc);
285
286    } else {
287       int c;
288       for (c=0;c<C;c++)
289       {
290          int i;
291          VARDECL(celt_ener_t, E);
292          SAVE_STACK;
293          ALLOC(E, m->nbEBands, celt_ener_t);
294          for (i=0;i<m->nbEBands;i++)
295             E[i] = eBands[C*i+c];
296          quant_coarse_energy_mono(m, E, oldEBands+c*m->nbEBands, budget/C, prob, error+c*m->nbEBands, enc);
297       }
298       RESTORE_STACK;
299    }
300 }
301
302 void quant_fine_energy(const CELTMode *m, celt_ener_t *eBands, celt_word16_t *oldEBands, celt_word16_t *error, celt_int16_t *fine_quant, ec_enc *enc)
303 {
304    int C;
305    C = m->nbChannels;
306
307    if (C==1)
308    {
309       quant_fine_energy_mono(m, eBands, oldEBands, error, fine_quant, enc);
310
311    } else {
312       int c;
313       VARDECL(celt_ener_t, E);
314       ALLOC(E, m->nbEBands, celt_ener_t);
315       for (c=0;c<C;c++)
316       {
317          int i;
318          SAVE_STACK;
319          quant_fine_energy_mono(m, E, oldEBands+c*m->nbEBands, error+c*m->nbEBands, fine_quant, enc);
320          for (i=0;i<m->nbEBands;i++)
321             eBands[C*i+c] = E[i];
322       }
323       RESTORE_STACK;
324    }
325 }
326
327
328 void unquant_coarse_energy(const CELTMode *m, celt_ener_t *eBands, celt_word16_t *oldEBands, int budget, int *prob, ec_dec *dec)
329 {
330    int C;   
331
332    C = m->nbChannels;
333    if (C==1)
334    {
335       unquant_coarse_energy_mono(m, eBands, oldEBands, budget, prob, dec);
336    }
337    else {
338       int c;
339       VARDECL(celt_ener_t, E);
340       SAVE_STACK;
341       ALLOC(E, m->nbEBands, celt_ener_t);
342       for (c=0;c<C;c++)
343       {
344          unquant_coarse_energy_mono(m, E, oldEBands+c*m->nbEBands, budget/C, prob, dec);
345       }
346       RESTORE_STACK;
347    }
348 }
349
350 void unquant_fine_energy(const CELTMode *m, celt_ener_t *eBands, celt_word16_t *oldEBands, celt_int16_t *fine_quant, ec_dec *dec)
351 {
352    int C;   
353
354    C = m->nbChannels;
355
356    if (C==1)
357    {
358       unquant_fine_energy_mono(m, eBands, oldEBands, fine_quant, dec);
359    }
360    else {
361       int c;
362       VARDECL(celt_ener_t, E);
363       SAVE_STACK;
364       ALLOC(E, m->nbEBands, celt_ener_t);
365       for (c=0;c<C;c++)
366       {
367          int i;
368          unquant_fine_energy_mono(m, E, oldEBands+c*m->nbEBands, fine_quant, dec);
369          for (i=0;i<m->nbEBands;i++)
370             eBands[C*i+c] = E[i];
371       }
372       RESTORE_STACK;
373    }
374 }