03da8524ee618638fa5cdfc1f7c1feeffce69c19
[opus.git] / libcelt / quant_bands.c
1 /* Copyright (c) 2007-2008 CSIRO
2    Copyright (c) 2007-2009 Xiph.Org Foundation
3    Written by Jean-Marc Valin */
4 /*
5    Redistribution and use in source and binary forms, with or without
6    modification, are permitted provided that the following conditions
7    are met:
8    
9    - Redistributions of source code must retain the above copyright
10    notice, this list of conditions and the following disclaimer.
11    
12    - Redistributions in binary form must reproduce the above copyright
13    notice, this list of conditions and the following disclaimer in the
14    documentation and/or other materials provided with the distribution.
15    
16    - Neither the name of the Xiph.org Foundation nor the names of its
17    contributors may be used to endorse or promote products derived from
18    this software without specific prior written permission.
19    
20    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
21    ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
22    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
23    A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE FOUNDATION OR
24    CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
25    EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
26    PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
27    PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
28    LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
29    NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
30    SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
31 */
32
33 #ifdef HAVE_CONFIG_H
34 #include "config.h"
35 #endif
36
37 #include "quant_bands.h"
38 #include "laplace.h"
39 #include <math.h>
40 #include "os_support.h"
41 #include "arch.h"
42 #include "mathops.h"
43 #include "stack_alloc.h"
44
45 #ifdef FIXED_POINT
46 /* Mean energy in each band quantized in Q6 */
47 const signed char eMeans[25] = {
48      124,122,115,106,100,
49       95, 91, 90, 99, 96,
50       94, 93, 98, 91, 86,
51       90, 88, 88, 90, 85,
52       64, 64, 64, 64, 64};
53 #else
54 /* Mean energy in each band quantized in Q6 and converted back to float */
55 const celt_word16 eMeans[25] = {
56       7.750000f, 7.625000f, 7.187500f, 6.625000f, 6.250000f,
57       5.937500f, 5.687500f, 5.625000f, 6.187500f, 6.000000f,
58       5.875000f, 5.812500f, 6.125000f, 5.687500f, 5.375000f,
59       5.625000f, 5.500000f, 5.500000f, 5.625000f, 5.312500f,
60       4.000000f, 4.000000f, 4.000000f, 4.000000f, 4.000000f};
61 #endif
62 /* prediction coefficients: 0.9, 0.8, 0.65, 0.5 */
63 #ifdef FIXED_POINT
64 static const celt_word16 pred_coef[4] = {29440, 26112, 21248, 16384};
65 #else
66 static const celt_word16 pred_coef[4] = {29440/32768., 26112/32768., 21248/32768., 16384/32768.};
67 #endif
68
69 int intra_decision(celt_word16 *eBands, celt_word16 *oldEBands, int start, int end, int len, int C)
70 {
71    int c, i;
72    celt_word32 dist = 0;
73    for (c=0;c<C;c++)
74    {
75       for (i=start;i<end;i++)
76       {
77          celt_word16 d = SUB16(eBands[i+c*len], oldEBands[i+c*len]);
78          dist = MAC16_16(dist, d,d);
79       }
80    }
81    return SHR32(dist,2*DB_SHIFT) > 2*C*(end-start);
82 }
83
84 int *quant_prob_alloc(const CELTMode *m)
85 {
86    int i;
87    int *prob;
88    prob = celt_alloc(4*m->nbEBands*sizeof(int));
89    if (prob==NULL)
90      return NULL;
91    for (i=0;i<m->nbEBands;i++)
92    {
93       prob[2*i] = 7000-i*200;
94       prob[2*i+1] = ec_laplace_get_start_freq(prob[2*i]);
95    }
96    for (i=0;i<m->nbEBands;i++)
97    {
98       prob[2*m->nbEBands+2*i] = 9000-i*220;
99       prob[2*m->nbEBands+2*i+1] = ec_laplace_get_start_freq(prob[2*m->nbEBands+2*i]);
100    }
101    return prob;
102 }
103
104 void quant_prob_free(int *freq)
105 {
106    celt_free(freq);
107 }
108
109 void quant_coarse_energy(const CELTMode *m, int start, int end, const celt_word16 *eBands, celt_word16 *oldEBands, int budget, int intra, int *prob, celt_word16 *error, ec_enc *enc, int _C, int LM, celt_word16 max_decay)
110 {
111    int i, c;
112    celt_word32 prev[2] = {0,0};
113    celt_word16 coef;
114    celt_word16 beta;
115    const int C = CHANNELS(_C);
116
117    coef = pred_coef[LM];
118
119    if (intra)
120    {
121       coef = 0;
122       prob += 2*m->nbEBands;
123    }
124    /* No theoretical justification for this, it just works */
125    beta = MULT16_16_P15(coef,coef);
126    /* Encode at a fixed coarse resolution */
127    for (i=start;i<end;i++)
128    {
129       c=0;
130       do {
131          int bits_left;
132          int qi;
133          celt_word16 q;
134          celt_word16 x;
135          celt_word32 f;
136          x = eBands[i+c*m->nbEBands];
137 #ifdef FIXED_POINT
138          f = SHL32(EXTEND32(x),15) -MULT16_16(coef,oldEBands[i+c*m->nbEBands])-prev[c];
139          /* Rounding to nearest integer here is really important! */
140          qi = (f+QCONST32(.5,DB_SHIFT+15))>>(DB_SHIFT+15);
141 #else
142          f = x-coef*oldEBands[i+c*m->nbEBands]-prev[c];
143          /* Rounding to nearest integer here is really important! */
144          qi = (int)floor(.5f+f);
145 #endif
146          /* Prevent the energy from going down too quickly (e.g. for bands
147             that have just one bin) */
148          if (qi < 0 && x < oldEBands[i+c*m->nbEBands]-max_decay)
149          {
150             qi += SHR16(oldEBands[i+c*m->nbEBands]-max_decay-x, DB_SHIFT);
151             if (qi > 0)
152                qi = 0;
153          }
154          /* If we don't have enough bits to encode all the energy, just assume something safe.
155             We allow slightly busting the budget here */
156          bits_left = budget-(int)ec_enc_tell(enc, 0)-2*C*(end-i);
157          if (bits_left < 24)
158          {
159             if (qi > 1)
160                qi = 1;
161             if (qi < -1)
162                qi = -1;
163             if (bits_left<8)
164                qi = 0;
165          }
166          ec_laplace_encode_start(enc, &qi, prob[2*i], prob[2*i+1]);
167          error[i+c*m->nbEBands] = PSHR32(f,15) - SHL16(qi,DB_SHIFT);
168          q = SHL16(qi,DB_SHIFT);
169          
170          oldEBands[i+c*m->nbEBands] = PSHR32(MULT16_16(coef,oldEBands[i+c*m->nbEBands]) + prev[c] + SHL32(EXTEND32(q),15), 15);
171          prev[c] = prev[c] + SHL32(EXTEND32(q),15) - MULT16_16(beta,q);
172       } while (++c < C);
173    }
174 }
175
176 void quant_fine_energy(const CELTMode *m, int start, int end, celt_ener *eBands, celt_word16 *oldEBands, celt_word16 *error, int *fine_quant, ec_enc *enc, int _C)
177 {
178    int i, c;
179    const int C = CHANNELS(_C);
180
181    /* Encode finer resolution */
182    for (i=start;i<end;i++)
183    {
184       celt_int16 frac = 1<<fine_quant[i];
185       if (fine_quant[i] <= 0)
186          continue;
187       c=0;
188       do {
189          int q2;
190          celt_word16 offset;
191 #ifdef FIXED_POINT
192          /* Has to be without rounding */
193          q2 = (error[i+c*m->nbEBands]+QCONST16(.5f,DB_SHIFT))>>(DB_SHIFT-fine_quant[i]);
194 #else
195          q2 = (int)floor((error[i+c*m->nbEBands]+.5f)*frac);
196 #endif
197          if (q2 > frac-1)
198             q2 = frac-1;
199          if (q2<0)
200             q2 = 0;
201          ec_enc_bits(enc, q2, fine_quant[i]);
202 #ifdef FIXED_POINT
203          offset = SUB16(SHR16(SHL16(q2,DB_SHIFT)+QCONST16(.5,DB_SHIFT),fine_quant[i]),QCONST16(.5f,DB_SHIFT));
204 #else
205          offset = (q2+.5f)*(1<<(14-fine_quant[i]))*(1.f/16384) - .5f;
206 #endif
207          oldEBands[i+c*m->nbEBands] += offset;
208          error[i+c*m->nbEBands] -= offset;
209          /*printf ("%f ", error[i] - offset);*/
210       } while (++c < C);
211    }
212 }
213
214 void quant_energy_finalise(const CELTMode *m, int start, int end, celt_ener *eBands, celt_word16 *oldEBands, celt_word16 *error, int *fine_quant, int *fine_priority, int bits_left, ec_enc *enc, int _C)
215 {
216    int i, prio, c;
217    const int C = CHANNELS(_C);
218
219    /* Use up the remaining bits */
220    for (prio=0;prio<2;prio++)
221    {
222       for (i=start;i<end && bits_left>=C ;i++)
223       {
224          if (fine_quant[i] >= 7 || fine_priority[i]!=prio)
225             continue;
226          c=0;
227          do {
228             int q2;
229             celt_word16 offset;
230             q2 = error[i+c*m->nbEBands]<0 ? 0 : 1;
231             ec_enc_bits(enc, q2, 1);
232 #ifdef FIXED_POINT
233             offset = SHR16(SHL16(q2,DB_SHIFT)-QCONST16(.5,DB_SHIFT),fine_quant[i]+1);
234 #else
235             offset = (q2-.5f)*(1<<(14-fine_quant[i]-1))*(1.f/16384);
236 #endif
237             oldEBands[i+c*m->nbEBands] += offset;
238             bits_left--;
239          } while (++c < C);
240       }
241    }
242 }
243
244 void unquant_coarse_energy(const CELTMode *m, int start, int end, celt_ener *eBands, celt_word16 *oldEBands, int intra, int *prob, ec_dec *dec, int _C, int LM)
245 {
246    int i, c;
247    celt_word32 prev[2] = {0, 0};
248    celt_word16 coef;
249    celt_word16 beta;
250    const int C = CHANNELS(_C);
251
252    coef = pred_coef[LM];
253
254    if (intra)
255    {
256       coef = 0;
257       prob += 2*m->nbEBands;
258    }
259    /* No theoretical justification for this, it just works */
260    beta = MULT16_16_P15(coef,coef);
261
262    /* Decode at a fixed coarse resolution */
263    for (i=start;i<end;i++)
264    {
265       c=0;
266       do {
267          int qi;
268          celt_word16 q;
269          qi = ec_laplace_decode_start(dec, prob[2*i], prob[2*i+1]);
270          q = SHL16(qi,DB_SHIFT);
271
272          oldEBands[i+c*m->nbEBands] = PSHR32(MULT16_16(coef,oldEBands[i+c*m->nbEBands]) + prev[c] + SHL32(EXTEND32(q),15), 15);
273          prev[c] = prev[c] + SHL32(EXTEND32(q),15) - MULT16_16(beta,q);
274       } while (++c < C);
275    }
276 }
277
278 void unquant_fine_energy(const CELTMode *m, int start, int end, celt_ener *eBands, celt_word16 *oldEBands, int *fine_quant, ec_dec *dec, int _C)
279 {
280    int i, c;
281    const int C = CHANNELS(_C);
282    /* Decode finer resolution */
283    for (i=start;i<end;i++)
284    {
285       if (fine_quant[i] <= 0)
286          continue;
287       c=0; 
288       do {
289          int q2;
290          celt_word16 offset;
291          q2 = ec_dec_bits(dec, fine_quant[i]);
292 #ifdef FIXED_POINT
293          offset = SUB16(SHR16(SHL16(q2,DB_SHIFT)+QCONST16(.5,DB_SHIFT),fine_quant[i]),QCONST16(.5f,DB_SHIFT));
294 #else
295          offset = (q2+.5f)*(1<<(14-fine_quant[i]))*(1.f/16384) - .5f;
296 #endif
297          oldEBands[i+c*m->nbEBands] += offset;
298       } while (++c < C);
299    }
300 }
301
302 void unquant_energy_finalise(const CELTMode *m, int start, int end, celt_ener *eBands, celt_word16 *oldEBands, int *fine_quant,  int *fine_priority, int bits_left, ec_dec *dec, int _C)
303 {
304    int i, prio, c;
305    const int C = CHANNELS(_C);
306
307    /* Use up the remaining bits */
308    for (prio=0;prio<2;prio++)
309    {
310       for (i=start;i<end && bits_left>=C ;i++)
311       {
312          if (fine_quant[i] >= 7 || fine_priority[i]!=prio)
313             continue;
314          c=0;
315          do {
316             int q2;
317             celt_word16 offset;
318             q2 = ec_dec_bits(dec, 1);
319 #ifdef FIXED_POINT
320             offset = SHR16(SHL16(q2,DB_SHIFT)-QCONST16(.5,DB_SHIFT),fine_quant[i]+1);
321 #else
322             offset = (q2-.5f)*(1<<(14-fine_quant[i]-1))*(1.f/16384);
323 #endif
324             oldEBands[i+c*m->nbEBands] += offset;
325             bits_left--;
326          } while (++c < C);
327       }
328    }
329 }
330
331 void log2Amp(const CELTMode *m, int start, int end,
332       celt_ener *eBands, celt_word16 *oldEBands, int _C)
333 {
334    int c, i;
335    const int C = CHANNELS(_C);
336    c=0;
337    do {
338       for (i=start;i<m->nbEBands;i++)
339       {
340          celt_word16 lg = oldEBands[i+c*m->nbEBands]
341                         + SHL16((celt_word16)eMeans[i],6);
342          eBands[i+c*m->nbEBands] = PSHR32(celt_exp2(SHL16(lg,11-DB_SHIFT)),4);
343          if (oldEBands[i+c*m->nbEBands] < -QCONST16(14.f,DB_SHIFT))
344             oldEBands[i+c*m->nbEBands] = -QCONST16(14.f,DB_SHIFT);
345       }
346    } while (++c < C);
347 }
348
349 void amp2Log2(const CELTMode *m, int effEnd, int end,
350       celt_ener *bandE, celt_word16 *bandLogE, int _C)
351 {
352    int c, i;
353    const int C = CHANNELS(_C);
354    c=0;
355    do {
356       for (i=0;i<effEnd;i++)
357          bandLogE[i+c*m->nbEBands] =
358                celt_log2(MAX32(QCONST32(.001f,14),SHL32(bandE[i+c*m->nbEBands],2)))
359                - SHL16((celt_word16)eMeans[i],6);
360       for (i=effEnd;i<end;i++)
361          bandLogE[c*m->nbEBands+i] = -QCONST16(14.f,DB_SHIFT);
362    } while (++c < C);
363 }