Using the real spectral means instead of the ones
[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 const celt_word16 eMeans[25] = {
47       7941, 7777, 7344, 6791, 6397,
48       6076, 5825, 5773, 6305, 6151,
49       6030, 5922, 6290, 5842, 5525,
50       5733, 5604, 5659, 5732, 5445,
51       4082, 4082, 4082, 4082, 4082};
52 #else
53 const celt_word16 eMeans[25] = {
54       7.755326, 7.594506, 7.172360, 6.632112, 6.247387,
55       5.933998, 5.688906, 5.637953, 6.157458, 6.006739,
56       5.889151, 5.783105, 6.142725, 5.704652, 5.395896,
57       5.598698, 5.472708, 5.526389, 5.597547, 5.317134,
58       3.986353, 3.986353, 3.986353, 3.986353, 3.986353};
59 #endif
60 /* prediction coefficients: 0.9, 0.8, 0.65, 0.5 */
61 #ifdef FIXED_POINT
62 static const celt_word16 pred_coef[4] = {29440, 26112, 21248, 16384};
63 #else
64 static const celt_word16 pred_coef[4] = {29440/32768., 26112/32768., 21248/32768., 16384/32768.};
65 #endif
66
67 int intra_decision(celt_word16 *eBands, celt_word16 *oldEBands, int start, int end, int len, int C)
68 {
69    int c, i;
70    celt_word32 dist = 0;
71    for (c=0;c<C;c++)
72    {
73       for (i=start;i<end;i++)
74       {
75          celt_word16 d = SUB16(eBands[i+c*len], oldEBands[i+c*len]);
76          dist = MAC16_16(dist, d,d);
77       }
78    }
79    return SHR32(dist,2*DB_SHIFT) > 2*C*(end-start);
80 }
81
82 int *quant_prob_alloc(const CELTMode *m)
83 {
84    int i;
85    int *prob;
86    prob = celt_alloc(4*m->nbEBands*sizeof(int));
87    if (prob==NULL)
88      return NULL;
89    for (i=0;i<m->nbEBands;i++)
90    {
91       prob[2*i] = 7000-i*200;
92       prob[2*i+1] = ec_laplace_get_start_freq(prob[2*i]);
93    }
94    for (i=0;i<m->nbEBands;i++)
95    {
96       prob[2*m->nbEBands+2*i] = 9000-i*220;
97       prob[2*m->nbEBands+2*i+1] = ec_laplace_get_start_freq(prob[2*m->nbEBands+2*i]);
98    }
99    return prob;
100 }
101
102 void quant_prob_free(int *freq)
103 {
104    celt_free(freq);
105 }
106
107 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)
108 {
109    int i, c;
110    celt_word32 prev[2] = {0,0};
111    celt_word16 coef;
112    celt_word16 beta;
113    const int C = CHANNELS(_C);
114
115    coef = pred_coef[LM];
116
117    if (intra)
118    {
119       coef = 0;
120       prob += 2*m->nbEBands;
121    }
122    /* No theoretical justification for this, it just works */
123    beta = MULT16_16_P15(coef,coef);
124    /* Encode at a fixed coarse resolution */
125    for (i=start;i<end;i++)
126    {
127       c=0;
128       do {
129          int bits_left;
130          int qi;
131          celt_word16 q;
132          celt_word16 x;
133          celt_word32 f;
134          x = eBands[i+c*m->nbEBands];
135 #ifdef FIXED_POINT
136          f = SHL32(EXTEND32(x),15) -MULT16_16(coef,oldEBands[i+c*m->nbEBands])-prev[c];
137          /* Rounding to nearest integer here is really important! */
138          qi = (f+QCONST32(.5,DB_SHIFT+15))>>(DB_SHIFT+15);
139 #else
140          f = x-coef*oldEBands[i+c*m->nbEBands]-prev[c];
141          /* Rounding to nearest integer here is really important! */
142          qi = (int)floor(.5f+f);
143 #endif
144          /* Prevent the energy from going down too quickly (e.g. for bands
145             that have just one bin) */
146          if (qi < 0 && x < oldEBands[i+c*m->nbEBands]-max_decay)
147          {
148             qi += SHR16(oldEBands[i+c*m->nbEBands]-max_decay-x, DB_SHIFT);
149             if (qi > 0)
150                qi = 0;
151          }
152          /* If we don't have enough bits to encode all the energy, just assume something safe.
153             We allow slightly busting the budget here */
154          bits_left = budget-(int)ec_enc_tell(enc, 0)-2*C*(end-i);
155          if (bits_left < 24)
156          {
157             if (qi > 1)
158                qi = 1;
159             if (qi < -1)
160                qi = -1;
161             if (bits_left<8)
162                qi = 0;
163          }
164          ec_laplace_encode_start(enc, &qi, prob[2*i], prob[2*i+1]);
165          error[i+c*m->nbEBands] = PSHR32(f,15) - SHL16(qi,DB_SHIFT);
166          q = SHL16(qi,DB_SHIFT);
167          
168          oldEBands[i+c*m->nbEBands] = PSHR32(MULT16_16(coef,oldEBands[i+c*m->nbEBands]) + prev[c] + SHL32(EXTEND32(q),15), 15);
169          prev[c] = prev[c] + SHL32(EXTEND32(q),15) - MULT16_16(beta,q);
170       } while (++c < C);
171    }
172 }
173
174 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)
175 {
176    int i, c;
177    const int C = CHANNELS(_C);
178
179    /* Encode finer resolution */
180    for (i=start;i<end;i++)
181    {
182       celt_int16 frac = 1<<fine_quant[i];
183       if (fine_quant[i] <= 0)
184          continue;
185       c=0;
186       do {
187          int q2;
188          celt_word16 offset;
189 #ifdef FIXED_POINT
190          /* Has to be without rounding */
191          q2 = (error[i+c*m->nbEBands]+QCONST16(.5f,DB_SHIFT))>>(DB_SHIFT-fine_quant[i]);
192 #else
193          q2 = (int)floor((error[i+c*m->nbEBands]+.5f)*frac);
194 #endif
195          if (q2 > frac-1)
196             q2 = frac-1;
197          if (q2<0)
198             q2 = 0;
199          ec_enc_bits(enc, q2, fine_quant[i]);
200 #ifdef FIXED_POINT
201          offset = SUB16(SHR16(SHL16(q2,DB_SHIFT)+QCONST16(.5,DB_SHIFT),fine_quant[i]),QCONST16(.5f,DB_SHIFT));
202 #else
203          offset = (q2+.5f)*(1<<(14-fine_quant[i]))*(1.f/16384) - .5f;
204 #endif
205          oldEBands[i+c*m->nbEBands] += offset;
206          error[i+c*m->nbEBands] -= offset;
207          /*printf ("%f ", error[i] - offset);*/
208       } while (++c < C);
209    }
210 }
211
212 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)
213 {
214    int i, prio, c;
215    const int C = CHANNELS(_C);
216
217    /* Use up the remaining bits */
218    for (prio=0;prio<2;prio++)
219    {
220       for (i=start;i<end && bits_left>=C ;i++)
221       {
222          if (fine_quant[i] >= 7 || fine_priority[i]!=prio)
223             continue;
224          c=0;
225          do {
226             int q2;
227             celt_word16 offset;
228             q2 = error[i+c*m->nbEBands]<0 ? 0 : 1;
229             ec_enc_bits(enc, q2, 1);
230 #ifdef FIXED_POINT
231             offset = SHR16(SHL16(q2,DB_SHIFT)-QCONST16(.5,DB_SHIFT),fine_quant[i]+1);
232 #else
233             offset = (q2-.5f)*(1<<(14-fine_quant[i]-1))*(1.f/16384);
234 #endif
235             oldEBands[i+c*m->nbEBands] += offset;
236             bits_left--;
237          } while (++c < C);
238       }
239    }
240 }
241
242 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)
243 {
244    int i, c;
245    celt_word32 prev[2] = {0, 0};
246    celt_word16 coef;
247    celt_word16 beta;
248    const int C = CHANNELS(_C);
249
250    coef = pred_coef[LM];
251
252    if (intra)
253    {
254       coef = 0;
255       prob += 2*m->nbEBands;
256    }
257    /* No theoretical justification for this, it just works */
258    beta = MULT16_16_P15(coef,coef);
259
260    /* Decode at a fixed coarse resolution */
261    for (i=start;i<end;i++)
262    {
263       c=0;
264       do {
265          int qi;
266          celt_word16 q;
267          qi = ec_laplace_decode_start(dec, prob[2*i], prob[2*i+1]);
268          q = SHL16(qi,DB_SHIFT);
269
270          oldEBands[i+c*m->nbEBands] = PSHR32(MULT16_16(coef,oldEBands[i+c*m->nbEBands]) + prev[c] + SHL32(EXTEND32(q),15), 15);
271          prev[c] = prev[c] + SHL32(EXTEND32(q),15) - MULT16_16(beta,q);
272       } while (++c < C);
273    }
274 }
275
276 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)
277 {
278    int i, c;
279    const int C = CHANNELS(_C);
280    /* Decode finer resolution */
281    for (i=start;i<end;i++)
282    {
283       if (fine_quant[i] <= 0)
284          continue;
285       c=0; 
286       do {
287          int q2;
288          celt_word16 offset;
289          q2 = ec_dec_bits(dec, fine_quant[i]);
290 #ifdef FIXED_POINT
291          offset = SUB16(SHR16(SHL16(q2,DB_SHIFT)+QCONST16(.5,DB_SHIFT),fine_quant[i]),QCONST16(.5f,DB_SHIFT));
292 #else
293          offset = (q2+.5f)*(1<<(14-fine_quant[i]))*(1.f/16384) - .5f;
294 #endif
295          oldEBands[i+c*m->nbEBands] += offset;
296       } while (++c < C);
297    }
298 }
299
300 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)
301 {
302    int i, prio, c;
303    const int C = CHANNELS(_C);
304
305    /* Use up the remaining bits */
306    for (prio=0;prio<2;prio++)
307    {
308       for (i=start;i<end && bits_left>=C ;i++)
309       {
310          if (fine_quant[i] >= 7 || fine_priority[i]!=prio)
311             continue;
312          c=0;
313          do {
314             int q2;
315             celt_word16 offset;
316             q2 = ec_dec_bits(dec, 1);
317 #ifdef FIXED_POINT
318             offset = SHR16(SHL16(q2,DB_SHIFT)-QCONST16(.5,DB_SHIFT),fine_quant[i]+1);
319 #else
320             offset = (q2-.5f)*(1<<(14-fine_quant[i]-1))*(1.f/16384);
321 #endif
322             oldEBands[i+c*m->nbEBands] += offset;
323             bits_left--;
324          } while (++c < C);
325       }
326    }
327 }
328
329 void log2Amp(const CELTMode *m, int start, int end,
330       celt_ener *eBands, celt_word16 *oldEBands, int _C)
331 {
332    int c, i;
333    const int C = CHANNELS(_C);
334    c=0;
335    do {
336       for (i=start;i<m->nbEBands;i++)
337       {
338          celt_word16 lg = oldEBands[i+c*m->nbEBands]+eMeans[i];
339          eBands[i+c*m->nbEBands] = PSHR32(celt_exp2(SHL16(lg,11-DB_SHIFT)),4);
340          if (oldEBands[i+c*m->nbEBands] < -QCONST16(14.f,DB_SHIFT))
341             oldEBands[i+c*m->nbEBands] = -QCONST16(14.f,DB_SHIFT);
342       }
343    } while (++c < C);
344 }
345
346 void amp2Log2(const CELTMode *m, int effEnd, int end,
347       celt_ener *bandE, celt_word16 *bandLogE, int _C)
348 {
349    int c, i;
350    const int C = CHANNELS(_C);
351    c=0;
352    do {
353       for (i=0;i<effEnd;i++)
354          bandLogE[i+c*m->nbEBands] =
355                celt_log2(MAX32(QCONST32(.001f,14),SHL32(bandE[i+c*m->nbEBands],2)))
356                - eMeans[i];
357       for (i=effEnd;i<end;i++)
358          bandLogE[c*m->nbEBands+i] = -QCONST16(14.f,DB_SHIFT);
359    } while (++c < C);
360 }