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