summaryrefslogtreecommitdiff
path: root/apps/codecs/libspeex/math_approx.h
diff options
context:
space:
mode:
Diffstat (limited to 'apps/codecs/libspeex/math_approx.h')
-rw-r--r--apps/codecs/libspeex/math_approx.h275
1 files changed, 262 insertions, 13 deletions
diff --git a/apps/codecs/libspeex/math_approx.h b/apps/codecs/libspeex/math_approx.h
index 49cfda6ebe..8421d634bb 100644
--- a/apps/codecs/libspeex/math_approx.h
+++ b/apps/codecs/libspeex/math_approx.h
@@ -37,19 +37,7 @@
37 37
38#include "misc.h" 38#include "misc.h"
39 39
40spx_word16_t spx_cos(spx_word16_t x); 40#ifndef FIXED_POINT
41spx_int16_t spx_ilog2(spx_uint32_t x);
42spx_int16_t spx_ilog4(spx_uint32_t x);
43#ifdef FIXED_POINT
44spx_word16_t spx_sqrt(spx_word32_t x);
45spx_word16_t spx_acos(spx_word16_t x);
46spx_word32_t spx_exp(spx_word16_t x);
47spx_word16_t spx_cos_norm(spx_word32_t x);
48
49/* Input in Q15, output in Q14 */
50spx_word16_t spx_atan(spx_word32_t x);
51
52#else
53 41
54#define spx_sqrt sqrt 42#define spx_sqrt sqrt
55#define spx_acos acos 43#define spx_acos acos
@@ -59,4 +47,265 @@ spx_word16_t spx_atan(spx_word32_t x);
59 47
60#endif 48#endif
61 49
50
51
52static inline spx_int16_t spx_ilog2(spx_uint32_t x)
53{
54 int r=0;
55 if (x>=(spx_int32_t)65536)
56 {
57 x >>= 16;
58 r += 16;
59 }
60 if (x>=256)
61 {
62 x >>= 8;
63 r += 8;
64 }
65 if (x>=16)
66 {
67 x >>= 4;
68 r += 4;
69 }
70 if (x>=4)
71 {
72 x >>= 2;
73 r += 2;
74 }
75 if (x>=2)
76 {
77 r += 1;
78 }
79 return r;
80}
81
82static inline spx_int16_t spx_ilog4(spx_uint32_t x)
83{
84 int r=0;
85 if (x>=(spx_int32_t)65536)
86 {
87 x >>= 16;
88 r += 8;
89 }
90 if (x>=256)
91 {
92 x >>= 8;
93 r += 4;
94 }
95 if (x>=16)
96 {
97 x >>= 4;
98 r += 2;
99 }
100 if (x>=4)
101 {
102 r += 1;
103 }
104 return r;
105}
106
107#ifdef FIXED_POINT
108
109/* sqrt(x) ~= 0.22178 + 1.29227*x - 0.77070*x^2 + 0.25723*x^3 (for .25 < x < 1) */
110/*#define C0 3634
111#define C1 21173
112#define C2 -12627
113#define C3 4215*/
114
115/* sqrt(x) ~= 0.22178 + 1.29227*x - 0.77070*x^2 + 0.25659*x^3 (for .25 < x < 1) */
116#define C0 3634
117#define C1 21173
118#define C2 -12627
119#define C3 4204
120
121static inline spx_word16_t spx_sqrt(spx_word32_t x)
122{
123 int k;
124 spx_word32_t rt;
125 k = spx_ilog4(x)-6;
126 x = VSHR32(x, (k<<1));
127 rt = ADD16(C0, MULT16_16_Q14(x, ADD16(C1, MULT16_16_Q14(x, ADD16(C2, MULT16_16_Q14(x, (C3)))))));
128 rt = VSHR32(rt,7-k);
129 return rt;
130}
131
132/* log(x) ~= -2.18151 + 4.20592*x - 2.88938*x^2 + 0.86535*x^3 (for .5 < x < 1) */
133
134
135#define A1 16469
136#define A2 2242
137#define A3 1486
138
139static inline spx_word16_t spx_acos(spx_word16_t x)
140{
141 int s=0;
142 spx_word16_t ret;
143 spx_word16_t sq;
144 if (x<0)
145 {
146 s=1;
147 x = NEG16(x);
148 }
149 x = SUB16(16384,x);
150
151 x = x >> 1;
152 sq = MULT16_16_Q13(x, ADD16(A1, MULT16_16_Q13(x, ADD16(A2, MULT16_16_Q13(x, (A3))))));
153 ret = spx_sqrt(SHL32(EXTEND32(sq),13));
154
155 /*ret = spx_sqrt(67108864*(-1.6129e-04 + 2.0104e+00*f + 2.7373e-01*f*f + 1.8136e-01*f*f*f));*/
156 if (s)
157 ret = SUB16(25736,ret);
158 return ret;
159}
160
161
162#define K1 8192
163#define K2 -4096
164#define K3 340
165#define K4 -10
166
167static inline spx_word16_t spx_cos(spx_word16_t x)
168{
169 spx_word16_t x2;
170
171 if (x<12868)
172 {
173 x2 = MULT16_16_P13(x,x);
174 return ADD32(K1, MULT16_16_P13(x2, ADD32(K2, MULT16_16_P13(x2, ADD32(K3, MULT16_16_P13(K4, x2))))));
175 } else {
176 x = SUB16(25736,x);
177 x2 = MULT16_16_P13(x,x);
178 return SUB32(-K1, MULT16_16_P13(x2, ADD32(K2, MULT16_16_P13(x2, ADD32(K3, MULT16_16_P13(K4, x2))))));
179 }
180}
181
182#define L1 32767
183#define L2 -7651
184#define L3 8277
185#define L4 -626
186
187static inline spx_word16_t _spx_cos_pi_2(spx_word16_t x)
188{
189 spx_word16_t x2;
190
191 x2 = MULT16_16_P15(x,x);
192 return ADD16(1,MIN16(32766,ADD32(SUB16(L1,x2), MULT16_16_P15(x2, ADD32(L2, MULT16_16_P15(x2, ADD32(L3, MULT16_16_P15(L4, x2))))))));
193}
194
195static inline spx_word16_t spx_cos_norm(spx_word32_t x)
196{
197 x = x&0x0001ffff;
198 if (x>SHL32(EXTEND32(1), 16))
199 x = SUB32(SHL32(EXTEND32(1), 17),x);
200 if (x&0x00007fff)
201 {
202 if (x<SHL32(EXTEND32(1), 15))
203 {
204 return _spx_cos_pi_2(EXTRACT16(x));
205 } else {
206 return NEG32(_spx_cos_pi_2(EXTRACT16(65536-x)));
207 }
208 } else {
209 if (x&0x0000ffff)
210 return 0;
211 else if (x&0x0001ffff)
212 return -32767;
213 else
214 return 32767;
215 }
216}
217
218/*
219 K0 = 1
220 K1 = log(2)
221 K2 = 3-4*log(2)
222 K3 = 3*log(2) - 2
223*/
224#define D0 16384
225#define D1 11356
226#define D2 3726
227#define D3 1301
228/* Input in Q11 format, output in Q16 */
229static inline spx_word32_t spx_exp2(spx_word16_t x)
230{
231 int integer;
232 spx_word16_t frac;
233 integer = SHR16(x,11);
234 if (integer>14)
235 return 0x7fffffff;
236 else if (integer < -15)
237 return 0;
238 frac = SHL16(x-SHL16(integer,11),3);
239 frac = ADD16(D0, MULT16_16_Q14(frac, ADD16(D1, MULT16_16_Q14(frac, ADD16(D2 , MULT16_16_Q14(D3,frac))))));
240 return VSHR32(EXTEND32(frac), -integer-2);
241}
242
243/* Input in Q11 format, output in Q16 */
244static inline spx_word32_t spx_exp(spx_word16_t x)
245{
246 if (x>21290)
247 return 0x7fffffff;
248 else if (x<-21290)
249 return 0;
250 else
251 return spx_exp2(MULT16_16_P14(23637,x));
252}
253#define M1 32767
254#define M2 -21
255#define M3 -11943
256#define M4 4936
257
258static inline spx_word16_t spx_atan01(spx_word16_t x)
259{
260 return MULT16_16_P15(x, ADD32(M1, MULT16_16_P15(x, ADD32(M2, MULT16_16_P15(x, ADD32(M3, MULT16_16_P15(M4, x)))))));
261}
262
263#undef M1
264#undef M2
265#undef M3
266#undef M4
267
268/* Input in Q15, output in Q14 */
269static inline spx_word16_t spx_atan(spx_word32_t x)
270{
271 if (x <= 32767)
272 {
273 return SHR16(spx_atan01(x),1);
274 } else {
275 int e = spx_ilog2(x);
276 if (e>=29)
277 return 25736;
278 x = DIV32_16(SHL32(EXTEND32(32767),29-e), EXTRACT16(SHR32(x, e-14)));
279 return SUB16(25736, SHR16(spx_atan01(x),1));
280 }
281}
282#else
283
284#ifndef M_PI
285#define M_PI 3.14159265358979323846 /* pi */
286#endif
287
288#define C1 0.9999932946f
289#define C2 -0.4999124376f
290#define C3 0.0414877472f
291#define C4 -0.0012712095f
292
293
294#define SPX_PI_2 1.5707963268
295static inline spx_word16_t spx_cos(spx_word16_t x)
296{
297 if (x<SPX_PI_2)
298 {
299 x *= x;
300 return C1 + x*(C2+x*(C3+C4*x));
301 } else {
302 x = M_PI-x;
303 x *= x;
304 return NEG16(C1 + x*(C2+x*(C3+C4*x)));
305 }
306}
307
308#endif
309
310
62#endif 311#endif