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