Add ff_dot_productf() to celp_math.{c,h}
[libav.git] / 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(arg <= 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 }