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 */ 34static 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 */ 52static 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 */ 69static 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 77int16_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 87int 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 */ 109static 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 121static 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 129static 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 137int16_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 147int 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 */ 164static 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 179int 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 199float 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} 209