libavcodec/celp_math.c

   1 /*
   2  * Various fixed-point math operations
   3  *
   4  * Copyright (c) 2008 Vladimir Voroshilov
   5  *
   6  * This file is part of FFmpeg.
   7  *
   8  * FFmpeg is free software; you can redistribute it and/or
   9  * modify it under the terms of the GNU Lesser General Public
  10  * License as published by the Free Software Foundation; either
  11  * version 2.1 of the License, or (at your option) any later version.
  12  *
  13  * FFmpeg is distributed in the hope that it will be useful,
  14  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  15  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  16  * Lesser General Public License for more details.
  17  *
  18  * You should have received a copy of the GNU Lesser General Public
  19  * License along with FFmpeg; if not, write to the Free Software
  20  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
  21  */
  22
  23 #include <inttypes.h>
  24 #include <limits.h>
  25 #include <assert.h>
  26
  27 #include "avcodec.h"
  28 #include "celp_math.h"
  29
  30 #ifdef G729_BITEXACT
  31 /**
  32  * Cosine table: base_cos[i] = (1<<15) * cos(i*PI/64)
  33  */
  34 static const int16_t base_cos[64] =
  35 {
  36   32767,  32729,  32610,  32413,  32138,  31786,  31357,  30853,
  37   30274,  29622,  28899,  28106,  27246,  26320,  25330,  24279,
  38   23170,  22006,  20788,  19520,  18205,  16846,  15447,  14010,
  39   12540,  11039,   9512,   7962,   6393,   4808,   3212,   1608,
  40       0,  -1608,  -3212,  -4808,  -6393,  -7962,  -9512, -11039,
  41  -12540, -14010, -15447, -16846, -18205, -19520, -20788, -22006,
  42  -23170, -24279, -25330, -26320, -27246, -28106, -28899, -29622,
  43  -30274, -30853, -31357, -31786, -32138, -32413, -32610, -32729
  44 };
  45
  46 /**
  47  * Slope used to compute cos(x)
  48  *
  49  * cos(ind*64+offset) = base_cos[ind]+offset*slope_cos[ind]
  50  * values multiplied by 1<<19
  51  */
  52 static const int16_t slope_cos[64] =
  53 {
  54    -632,  -1893,  -3150,  -4399,  -5638,  -6863,  -8072,  -9261,
  55  -10428, -11570, -12684, -13767, -14817, -15832, -16808, -17744,
  56  -18637, -19486, -20287, -21039, -21741, -22390, -22986, -23526,
  57  -24009, -24435, -24801, -25108, -25354, -25540, -25664, -25726,
  58  -25726, -25664, -25540, -25354, -25108, -24801, -24435, -24009,
  59  -23526, -22986, -22390, -21741, -21039, -20287, -19486, -18637,
  60  -17744, -16808, -15832, -14817, -13767, -12684, -11570, -10428,
  61   -9261,  -8072,  -6863,  -5638,  -4399,  -3150,  -1893,   -632
  62 };
  63
  64 /**
  65  * Table used to compute exp2(x)
  66  *
  67  * tab_exp2[i] = (1<<14) * exp2(i/32) = 2^(i/32) i=0..32
  68  */
  69 static const uint16_t tab_exp2[33] =
  70 {
  71   16384, 16743, 17109, 17484, 17867, 18258, 18658, 19066, 19484, 19911,
  72   20347, 20792, 21247, 21713, 22188, 22674, 23170, 23678, 24196, 24726,
  73   25268, 25821, 26386, 26964, 27554, 28158, 28774, 29405, 30048, 30706,
  74   31379, 32066, 32767
  75 };
  76
  77 int16_t ff_cos(uint16_t arg)
  78 {
  79     uint8_t offset= arg;
  80     uint8_t ind = arg >> 8;
  81
  82     assert(arg < 0x4000);
  83
  84     return FFMAX(base_cos[ind] + ((slope_cos[ind] * offset) >> 12), -0x8000);
  85 }
  86
  87 int ff_exp2(uint16_t power)
  88 {
  89     uint16_t frac_x0;
  90     uint16_t frac_dx;
  91     int result;
  92
  93     assert(power <= 0x7fff);
  94
  95     frac_x0 = power >> 10;
  96     frac_dx = (power & 0x03ff) << 5;
  97
  98     result = tab_exp2[frac_x0] << 15;
  99     result += frac_dx * (tab_exp2[frac_x0+1] - tab_exp2[frac_x0]);
 100
 101     return result >> 10;
 102 }
 103
 104 #else // G729_BITEXACT
 105
 106 /**
 107  * Cosine table: base_cos[i] = (1<<15) * cos(i*PI/64)
 108  */
 109 static const int16_t tab_cos[65] =
 110 {
 111   32767,  32738,  32617,  32421,  32145,  31793,  31364,  30860,
 112   30280,  29629,  28905,  28113,  27252,  26326,  25336,  24285,
 113   23176,  22011,  20793,  19525,  18210,  16851,  15451,  14014,
 114   12543,  11043,   9515,   7965,   6395,   4810,   3214,   1609,
 115       1,  -1607,  -3211,  -4808,  -6393,  -7962,  -9513, -11040,
 116  -12541, -14012, -15449, -16848, -18207, -19523, -20791, -22009,
 117  -23174, -24283, -25334, -26324, -27250, -28111, -28904, -29627,
 118  -30279, -30858, -31363, -31792, -32144, -32419, -32616, -32736, -32768,
 119 };
 120
 121 static const uint16_t exp2a[]=
 122 {
 123      0,  1435,  2901,  4400,  5931,  7496,  9096, 10730,
 124  12400, 14106, 15850, 17632, 19454, 21315, 23216, 25160,
 125  27146, 29175, 31249, 33368, 35534, 37747, 40009, 42320,
 126  44682, 47095, 49562, 52082, 54657, 57289, 59979, 62727,
 127 };
 128
 129 static const uint16_t exp2b[]=
 130 {
 131      3,   712,  1424,  2134,  2845,  3557,  4270,  4982,
 132   5696,  6409,  7124,  7839,  8554,  9270,  9986, 10704,
 133  11421, 12138, 12857, 13576, 14295, 15014, 15734, 16455,
 134  17176, 17898, 18620, 19343, 20066, 20790, 21514, 22238,
 135 };
 136
 137 int16_t ff_cos(uint16_t arg)
 138 {
 139     uint8_t offset= arg;
 140     uint8_t ind = arg >> 8;
 141
 142     assert(arg <= 0x3fff);
 143
 144     return tab_cos[ind] + (offset * (tab_cos[ind+1] - tab_cos[ind]) >> 8);
 145 }
 146
 147 int ff_exp2(uint16_t power)
 148 {
 149     unsigned int result= exp2a[power>>10] + 0x10000;
 150
 151     assert(power <= 0x7fff);
 152
 153     result= (result<<3) + ((result*exp2b[(power>>5)&31])>>17);
 154     return result + ((result*(power&31)*89)>>22);
 155 }
 156
 157 #endif // else G729_BITEXACT
 158
 159 /**
 160  * Table used to compute log2(x)
 161  *
 162  * tab_log2[i] = (1<<15) * log2(1 + i/32), i=0..32
 163  */
 164 static const uint16_t tab_log2[33] =
 165 {
 166 #ifdef G729_BITEXACT
 167       0,   1455,   2866,   4236,   5568,   6863,   8124,   9352,
 168   10549,  11716,  12855,  13967,  15054,  16117,  17156,  18172,
 169   19167,  20142,  21097,  22033,  22951,  23852,  24735,  25603,
 170   26455,  27291,  28113,  28922,  29716,  30497,  31266,  32023,  32767,
 171 #else
 172       4,   1459,   2870,   4240,   5572,   6867,   8127,   9355,
 173   10552,  11719,  12858,  13971,  15057,  16120,  17158,  18175,
 174   19170,  20145,  21100,  22036,  22954,  23854,  24738,  25605,
 175   26457,  27294,  28116,  28924,  29719,  30500,  31269,  32025,  32769,
 176 #endif
 177 };
 178
 179 int ff_log2(uint32_t value)
 180 {
 181     uint8_t  power_int;
 182     uint8_t  frac_x0;
 183     uint16_t frac_dx;
 184
 185     // Stripping zeros from beginning
 186     power_int = av_log2(value);
 187     value <<= (31 - power_int);
 188
 189     // b31 is always non-zero now
 190     frac_x0 = (value & 0x7c000000) >> 26; // b26-b31 and [32..63] -> [0..31]
 191     frac_dx = (value & 0x03fff800) >> 11;
 192
 193     value = tab_log2[frac_x0];
 194     value += (frac_dx * (tab_log2[frac_x0+1] - tab_log2[frac_x0])) >> 15;
 195
 196     return (power_int << 15) + value;
 197 }
 198
 199 float ff_dot_productf(const float* a, const float* b, int length)
 200 {
 201     float sum = 0;
 202     int i;
 203
 204     for(i=0; i<length; i++)
 205         sum += a[i] * b[i];
 206
 207     return sum;
 208 }