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