Max delta: +/- 16384
[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 static const signed char eMeans[25] = {
48       103,100, 92, 85, 81,
49        77, 72, 70, 78, 75,
50        73, 71, 78, 74, 69,
51        72, 70, 74, 76, 71,
52        60, 60, 60, 60, 60
53 };
54 #else
55 /* Mean energy in each band quantized in Q6 and converted back to float */
56 static const celt_word16 eMeans[25] = {
57       6.437500f, 6.250000f, 5.750000f, 5.312500f, 5.062500f,
58       4.812500f, 4.500000f, 4.375000f, 4.875000f, 4.687500f,
59       4.562500f, 4.437500f, 4.875000f, 4.625000f, 4.312500f,
60       4.500000f, 4.375000f, 4.625000f, 4.750000f, 4.437500f,
61       3.750000f, 3.750000f, 3.750000f, 3.750000f, 3.750000f
62 };
63 #endif
64 /* prediction coefficients: 0.9, 0.8, 0.65, 0.5 */
65 #ifdef FIXED_POINT
66 static const celt_word16 pred_coef[4] = {29440, 26112, 21248, 16384};
67 static const celt_word16 beta_coef[4] = {30147, 22282, 12124, 6554};
68 #else
69 static const celt_word16 pred_coef[4] = {29440/32768., 26112/32768., 21248/32768., 16384/32768.};
70 static const celt_word16 beta_coef[4] = {30147/32768., 22282/32768., 12124/32768., 6554/32768.};
71 #endif
72
73 /*Parameters of the Laplace-like probability models used for the coarse energy.
74   There is one pair of parameters for each frame size, prediction type
75    (inter/intra), and band number.
76   The first number of each pair is the probability of 0, and the second is the
77    decay rate, both in Q8 precision.*/
78 static const unsigned char e_prob_model[4][2][42] = {
79    /*120 sample frames.*/
80    {
81       /*Inter*/
82       {
83           72, 127,  65, 129,  66, 128,  65, 128,  64, 128,  62, 128,  64, 128,
84           64, 128,  92,  78,  92,  79,  92,  78,  90,  79, 116,  41, 115,  40,
85          114,  40, 132,  26, 132,  26, 145,  17, 161,  12, 176,  10, 177,  11
86       },
87       /*Intra*/
88       {
89           24, 179,  48, 138,  54, 135,  54, 132,  53, 134,  56, 133,  55, 132,
90           55, 132,  61, 114,  70,  96,  74,  88,  75,  88,  87,  74,  89,  66,
91           91,  67, 100,  59, 108,  50, 120,  40, 122,  37,  97,  43,  78,  50
92       }
93    },
94    /*240 sample frames.*/
95    {
96       /*Inter*/
97       {
98           83,  78,  84,  81,  88,  75,  86,  74,  87,  71,  90,  73,  93,  74,
99           93,  74, 109,  40, 114,  36, 117,  34, 117,  34, 143,  17, 145,  18,
100          146,  19, 162,  12, 165,  10, 178,   7, 189,   6, 190,   8, 177,   9
101       },
102       /*Intra*/
103       {
104           23, 178,  54, 115,  63, 102,  66,  98,  69,  99,  74,  89,  71,  91,
105           73,  91,  78,  89,  86,  80,  92,  66,  93,  64, 102,  59, 103,  60,
106          104,  60, 117,  52, 123,  44, 138,  35, 133,  31,  97,  38,  77,  45
107       }
108    },
109    /*480 sample frames.*/
110    {
111       /*Inter*/
112       {
113           61,  90,  93,  60, 105,  42, 107,  41, 110,  45, 116,  38, 113,  38,
114          112,  38, 124,  26, 132,  27, 136,  19, 140,  20, 155,  14, 159,  16,
115          158,  18, 170,  13, 177,  10, 187,   8, 192,   6, 175,   9, 159,  10
116       },
117       /*Intra*/
118       {
119           21, 178,  59, 110,  71,  86,  75,  85,  84,  83,  91,  66,  88,  73,
120           87,  72,  92,  75,  98,  72, 105,  58, 107,  54, 115,  52, 114,  55,
121          112,  56, 129,  51, 132,  40, 150,  33, 140,  29,  98,  35,  77,  42
122       }
123    },
124    /*960 sample frames.*/
125    {
126       /*Inter*/
127       {
128           42, 121,  96,  66, 108,  43, 111,  40, 117,  44, 123,  32, 120,  36,
129          119,  33, 127,  33, 134,  34, 139,  21, 147,  23, 152,  20, 158,  25,
130          154,  26, 166,  21, 173,  16, 184,  13, 184,  10, 150,  13, 139,  15
131       },
132       /*Intra*/
133       {
134           22, 178,  63, 114,  74,  82,  84,  83,  92,  82, 103,  62,  96,  72,
135           96,  67, 101,  73, 107,  72, 113,  55, 118,  52, 125,  52, 118,  52,
136          117,  55, 135,  49, 137,  39, 157,  32, 145,  29,  97,  33,  77,  40
137       }
138    }
139 };
140
141 static const unsigned char small_energy_icdf[3]={2,1,0};
142
143 static int intra_decision(const celt_word16 *eBands, celt_word16 *oldEBands, int start, int end, int len, int C)
144 {
145    int c, i;
146    celt_word32 dist = 0;
147    c=0; do {
148       for (i=start;i<end;i++)
149       {
150          celt_word16 d = SHR16(SUB16(eBands[i+c*len], oldEBands[i+c*len]),2);
151          dist = MAC16_16(dist, d,d);
152       }
153    } while (++c<C);
154    return SHR32(dist,2*DB_SHIFT-4) > 2*C*(end-start);
155 }
156
157 static void quant_coarse_energy_impl(const CELTMode *m, int start, int end,
158       const celt_word16 *eBands, celt_word16 *oldEBands,
159       ec_int32 budget, ec_int32 tell,
160       const unsigned char *prob_model, celt_word16 *error, ec_enc *enc,
161       int _C, int LM, int intra, celt_word16 max_decay)
162 {
163    const int C = CHANNELS(_C);
164    int i, c;
165    celt_word32 prev[2] = {0,0};
166    celt_word16 coef;
167    celt_word16 beta;
168
169    if (tell+3 <= budget)
170       ec_enc_bit_logp(enc, intra, 3);
171    if (intra)
172    {
173       coef = 0;
174       beta = QCONST16(.15f,15);
175    } else {
176       beta = beta_coef[LM];
177       coef = pred_coef[LM];
178    }
179
180    /* Encode at a fixed coarse resolution */
181    for (i=start;i<end;i++)
182    {
183       c=0;
184       do {
185          int bits_left;
186          int qi;
187          celt_word16 q;
188          celt_word16 x;
189          celt_word32 f;
190          x = eBands[i+c*m->nbEBands];
191 #ifdef FIXED_POINT
192          f = SHL32(EXTEND32(x),15) -MULT16_16(coef,oldEBands[i+c*m->nbEBands])-prev[c];
193          /* Rounding to nearest integer here is really important! */
194          qi = (f+QCONST32(.5,DB_SHIFT+15))>>(DB_SHIFT+15);
195 #else
196          f = x-coef*oldEBands[i+c*m->nbEBands]-prev[c];
197          /* Rounding to nearest integer here is really important! */
198          qi = (int)floor(.5f+f);
199 #endif
200          /* Prevent the energy from going down too quickly (e.g. for bands
201             that have just one bin) */
202          if (qi < 0 && x < oldEBands[i+c*m->nbEBands]-max_decay)
203          {
204             qi += (int)SHR16(oldEBands[i+c*m->nbEBands]-max_decay-x, DB_SHIFT);
205             if (qi > 0)
206                qi = 0;
207          }
208          /* If we don't have enough bits to encode all the energy, just assume
209              something safe. */
210          tell = ec_enc_tell(enc, 0);
211          bits_left = budget-tell-3*C*(end-i);
212          if (i!=start && bits_left < 30)
213          {
214             if (bits_left < 24)
215                qi = IMIN(1, qi);
216             if (bits_left < 16)
217                qi = IMAX(-1, qi);
218          }
219          if (budget-tell >= 15)
220          {
221             int pi;
222             pi = 2*IMIN(i,20);
223             ec_laplace_encode(enc, &qi,
224                   prob_model[pi]<<7, prob_model[pi+1]<<6);
225          }
226          else if(budget-tell >= 2)
227          {
228             qi = IMAX(-1, IMIN(qi, 1));
229             ec_enc_icdf(enc, 2*qi^-(qi<0), small_energy_icdf, 2);
230          }
231          else if(budget-tell >= 1)
232          {
233             qi = IMIN(0, qi);
234             ec_enc_bit_logp(enc, -qi, 1);
235          }
236          else
237             qi = -1;
238          error[i+c*m->nbEBands] = PSHR32(f,15) - SHL16(qi,DB_SHIFT);
239          q = SHL16(qi,DB_SHIFT);
240          
241          oldEBands[i+c*m->nbEBands] = PSHR32(MULT16_16(coef,oldEBands[i+c*m->nbEBands]) + prev[c] + SHL32(EXTEND32(q),15), 15);
242          prev[c] = prev[c] + SHL32(EXTEND32(q),15) - MULT16_16(beta,q);
243       } while (++c < C);
244    }
245 }
246
247 void quant_coarse_energy(const CELTMode *m, int start, int end, int effEnd,
248       const celt_word16 *eBands, celt_word16 *oldEBands, ec_uint32 budget,
249       celt_word16 *error, ec_enc *enc, int _C, int LM, int nbAvailableBytes,
250       int force_intra, int *delayedIntra, int two_pass)
251 {
252    const int C = CHANNELS(_C);
253    int intra;
254    celt_word16 max_decay;
255    VARDECL(celt_word16, oldEBands_intra);
256    VARDECL(celt_word16, error_intra);
257    ec_enc enc_start_state;
258    ec_byte_buffer buf_start_state;
259    ec_uint32 tell;
260    SAVE_STACK;
261
262    intra = force_intra || (*delayedIntra && nbAvailableBytes > end*C);
263    if (/*shortBlocks || */intra_decision(eBands, oldEBands, start, effEnd, m->nbEBands, C))
264       *delayedIntra = 1;
265    else
266       *delayedIntra = 0;
267
268    tell = ec_enc_tell(enc, 0);
269    if (tell+3 > budget)
270       two_pass = intra = 0;
271
272    /* Encode the global flags using a simple probability model
273       (first symbols in the stream) */
274
275 #ifdef FIXED_POINT
276       max_decay = MIN32(QCONST16(16,DB_SHIFT), SHL32(EXTEND32(nbAvailableBytes),DB_SHIFT-3));
277 #else
278    max_decay = MIN32(16.f, .125f*nbAvailableBytes);
279 #endif
280
281    enc_start_state = *enc;
282    buf_start_state = *(enc->buf);
283
284    ALLOC(oldEBands_intra, C*m->nbEBands, celt_word16);
285    ALLOC(error_intra, C*m->nbEBands, celt_word16);
286    CELT_COPY(oldEBands_intra, oldEBands, C*end);
287
288    if (two_pass || intra)
289    {
290       quant_coarse_energy_impl(m, start, end, eBands, oldEBands_intra, budget,
291             tell, e_prob_model[LM][1], error_intra, enc, C, LM, 1, max_decay);
292    }
293
294    if (!intra)
295    {
296       ec_enc enc_intra_state;
297       ec_byte_buffer buf_intra_state;
298       int tell_intra;
299       ec_uint32 nstart_bytes;
300       ec_uint32 nintra_bytes;
301       VARDECL(unsigned char, intra_bits);
302
303       tell_intra = ec_enc_tell(enc, 3);
304
305       enc_intra_state = *enc;
306       buf_intra_state = *(enc->buf);
307
308       nstart_bytes = ec_byte_bytes(&buf_start_state);
309       nintra_bytes = ec_byte_bytes(&buf_intra_state);
310       ALLOC(intra_bits, nintra_bytes-nstart_bytes, unsigned char);
311       /* Copy bits from intra bit-stream */
312       CELT_COPY(intra_bits,
313             ec_byte_get_buffer(&buf_intra_state) + nstart_bytes,
314             nintra_bytes - nstart_bytes);
315
316       *enc = enc_start_state;
317       *(enc->buf) = buf_start_state;
318
319       quant_coarse_energy_impl(m, start, end, eBands, oldEBands, budget,
320             tell, e_prob_model[LM][intra], error, enc, C, LM, 0, max_decay);
321
322       if (two_pass && ec_enc_tell(enc, 3) > tell_intra)
323       {
324          *enc = enc_intra_state;
325          *(enc->buf) = buf_intra_state;
326          /* Copy intra bits to bit-stream */
327          CELT_COPY(ec_byte_get_buffer(&buf_intra_state) + nstart_bytes,
328                intra_bits, nintra_bytes - nstart_bytes);
329          CELT_COPY(oldEBands, oldEBands_intra, C*end);
330          CELT_COPY(error, error_intra, C*end);
331       }
332    } else {
333       CELT_COPY(oldEBands, oldEBands_intra, C*end);
334       CELT_COPY(error, error_intra, C*end);
335    }
336    RESTORE_STACK;
337 }
338
339 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)
340 {
341    int i, c;
342    const int C = CHANNELS(_C);
343
344    /* Encode finer resolution */
345    for (i=start;i<end;i++)
346    {
347       celt_int16 frac = 1<<fine_quant[i];
348       if (fine_quant[i] <= 0)
349          continue;
350       c=0;
351       do {
352          int q2;
353          celt_word16 offset;
354 #ifdef FIXED_POINT
355          /* Has to be without rounding */
356          q2 = (error[i+c*m->nbEBands]+QCONST16(.5f,DB_SHIFT))>>(DB_SHIFT-fine_quant[i]);
357 #else
358          q2 = (int)floor((error[i+c*m->nbEBands]+.5f)*frac);
359 #endif
360          if (q2 > frac-1)
361             q2 = frac-1;
362          if (q2<0)
363             q2 = 0;
364          ec_enc_bits(enc, q2, fine_quant[i]);
365 #ifdef FIXED_POINT
366          offset = SUB16(SHR32(SHL32(EXTEND32(q2),DB_SHIFT)+QCONST16(.5,DB_SHIFT),fine_quant[i]),QCONST16(.5f,DB_SHIFT));
367 #else
368          offset = (q2+.5f)*(1<<(14-fine_quant[i]))*(1.f/16384) - .5f;
369 #endif
370          oldEBands[i+c*m->nbEBands] += offset;
371          error[i+c*m->nbEBands] -= offset;
372          /*printf ("%f ", error[i] - offset);*/
373       } while (++c < C);
374    }
375 }
376
377 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)
378 {
379    int i, prio, c;
380    const int C = CHANNELS(_C);
381
382    /* Use up the remaining bits */
383    for (prio=0;prio<2;prio++)
384    {
385       for (i=start;i<end && bits_left>=C ;i++)
386       {
387          if (fine_quant[i] >= 7 || fine_priority[i]!=prio)
388             continue;
389          c=0;
390          do {
391             int q2;
392             celt_word16 offset;
393             q2 = error[i+c*m->nbEBands]<0 ? 0 : 1;
394             ec_enc_bits(enc, q2, 1);
395 #ifdef FIXED_POINT
396             offset = SHR16(SHL16(q2,DB_SHIFT)-QCONST16(.5,DB_SHIFT),fine_quant[i]+1);
397 #else
398             offset = (q2-.5f)*(1<<(14-fine_quant[i]-1))*(1.f/16384);
399 #endif
400             oldEBands[i+c*m->nbEBands] += offset;
401             bits_left--;
402          } while (++c < C);
403       }
404    }
405 }
406
407 void unquant_coarse_energy(const CELTMode *m, int start, int end, celt_ener *eBands, celt_word16 *oldEBands, int intra, ec_dec *dec, int _C, int LM)
408 {
409    const unsigned char *prob_model = e_prob_model[LM][intra];
410    int i, c;
411    celt_word32 prev[2] = {0, 0};
412    celt_word16 coef;
413    celt_word16 beta;
414    const int C = CHANNELS(_C);
415    ec_int32 budget;
416    ec_int32 tell;
417
418
419    if (intra)
420    {
421       coef = 0;
422       beta = QCONST16(.15f,15);
423    } else {
424       beta = beta_coef[LM];
425       coef = pred_coef[LM];
426    }
427
428    budget = dec->buf->storage*8;
429
430    /* Decode at a fixed coarse resolution */
431    for (i=start;i<end;i++)
432    {
433       c=0;
434       do {
435          int qi;
436          celt_word16 q;
437          tell = ec_dec_tell(dec, 0);
438          if(budget-tell>=15)
439          {
440             int pi;
441             pi = 2*IMIN(i,20);
442             qi = ec_laplace_decode(dec,
443                   prob_model[pi]<<7, prob_model[pi+1]<<6);
444          }
445          else if(budget-tell>=2)
446          {
447             qi = ec_dec_icdf(dec, small_energy_icdf, 2);
448             qi = (qi>>1)^-(qi&1);
449          }
450          else if(budget-tell>=1)
451          {
452             qi = -ec_dec_bit_logp(dec, 1);
453          }
454          else
455             qi = -1;
456          q = SHL16(qi,DB_SHIFT);
457
458          oldEBands[i+c*m->nbEBands] = PSHR32(MULT16_16(coef,oldEBands[i+c*m->nbEBands]) + prev[c] + SHL32(EXTEND32(q),15), 15);
459          prev[c] = prev[c] + SHL32(EXTEND32(q),15) - MULT16_16(beta,q);
460       } while (++c < C);
461    }
462 }
463
464 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)
465 {
466    int i, c;
467    const int C = CHANNELS(_C);
468    /* Decode finer resolution */
469    for (i=start;i<end;i++)
470    {
471       if (fine_quant[i] <= 0)
472          continue;
473       c=0; 
474       do {
475          int q2;
476          celt_word16 offset;
477          q2 = ec_dec_bits(dec, fine_quant[i]);
478 #ifdef FIXED_POINT
479          offset = SUB16(SHR32(SHL32(EXTEND32(q2),DB_SHIFT)+QCONST16(.5,DB_SHIFT),fine_quant[i]),QCONST16(.5f,DB_SHIFT));
480 #else
481          offset = (q2+.5f)*(1<<(14-fine_quant[i]))*(1.f/16384) - .5f;
482 #endif
483          oldEBands[i+c*m->nbEBands] += offset;
484       } while (++c < C);
485    }
486 }
487
488 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)
489 {
490    int i, prio, c;
491    const int C = CHANNELS(_C);
492
493    /* Use up the remaining bits */
494    for (prio=0;prio<2;prio++)
495    {
496       for (i=start;i<end && bits_left>=C ;i++)
497       {
498          if (fine_quant[i] >= 8 || fine_priority[i]!=prio)
499             continue;
500          c=0;
501          do {
502             int q2;
503             celt_word16 offset;
504             q2 = ec_dec_bits(dec, 1);
505 #ifdef FIXED_POINT
506             offset = SHR16(SHL16(q2,DB_SHIFT)-QCONST16(.5,DB_SHIFT),fine_quant[i]+1);
507 #else
508             offset = (q2-.5f)*(1<<(14-fine_quant[i]-1))*(1.f/16384);
509 #endif
510             oldEBands[i+c*m->nbEBands] += offset;
511             bits_left--;
512          } while (++c < C);
513       }
514    }
515 }
516
517 void log2Amp(const CELTMode *m, int start, int end,
518       celt_ener *eBands, celt_word16 *oldEBands, int _C)
519 {
520    int c, i;
521    const int C = CHANNELS(_C);
522    c=0;
523    do {
524       for (i=start;i<m->nbEBands;i++)
525       {
526          celt_word16 lg = oldEBands[i+c*m->nbEBands]
527                         + SHL16((celt_word16)eMeans[i],6);
528          eBands[i+c*m->nbEBands] = PSHR32(celt_exp2(SHL16(lg,11-DB_SHIFT)),4);
529          if (oldEBands[i+c*m->nbEBands] < -QCONST16(14.f,DB_SHIFT))
530             oldEBands[i+c*m->nbEBands] = -QCONST16(14.f,DB_SHIFT);
531       }
532    } while (++c < C);
533 }
534
535 void amp2Log2(const CELTMode *m, int effEnd, int end,
536       celt_ener *bandE, celt_word16 *bandLogE, int _C)
537 {
538    int c, i;
539    const int C = CHANNELS(_C);
540    c=0;
541    do {
542       for (i=0;i<effEnd;i++)
543          bandLogE[i+c*m->nbEBands] =
544                celt_log2(MAX32(QCONST32(.001f,14),SHL32(bandE[i+c*m->nbEBands],2)))
545                - SHL16((celt_word16)eMeans[i],6);
546       for (i=effEnd;i<end;i++)
547          bandLogE[c*m->nbEBands+i] = -QCONST16(14.f,DB_SHIFT);
548    } while (++c < C);
549 }