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