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