summaryrefslogtreecommitdiff
path: root/apps/fixedpoint.c
diff options
context:
space:
mode:
Diffstat (limited to 'apps/fixedpoint.c')
-rw-r--r--apps/fixedpoint.c472
1 files changed, 0 insertions, 472 deletions
diff --git a/apps/fixedpoint.c b/apps/fixedpoint.c
deleted file mode 100644
index b212929245..0000000000
--- a/apps/fixedpoint.c
+++ /dev/null
@@ -1,472 +0,0 @@
1/***************************************************************************
2 * __________ __ ___.
3 * Open \______ \ ____ ____ | | _\_ |__ _______ ___
4 * Source | _// _ \_/ ___\| |/ /| __ \ / _ \ \/ /
5 * Jukebox | | ( <_> ) \___| < | \_\ ( <_> > < <
6 * Firmware |____|_ /\____/ \___ >__|_ \|___ /\____/__/\_ \
7 * \/ \/ \/ \/ \/
8 * $Id$
9 *
10 * Copyright (C) 2006 Jens Arnold
11 *
12 * Fixed point library for plugins
13 *
14 * This program is free software; you can redistribute it and/or
15 * modify it under the terms of the GNU General Public License
16 * as published by the Free Software Foundation; either version 2
17 * of the License, or (at your option) any later version.
18 *
19 * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY OF ANY
20 * KIND, either express or implied.
21 *
22 ****************************************************************************/
23
24#include "fixedpoint.h"
25#include <stdlib.h>
26#include <stdbool.h>
27#include <inttypes.h>
28
29#ifndef BIT_N
30#define BIT_N(n) (1U << (n))
31#endif
32
33/** TAKEN FROM ORIGINAL fixedpoint.h */
34/* Inverse gain of circular cordic rotation in s0.31 format. */
35static const long cordic_circular_gain = 0xb2458939; /* 0.607252929 */
36
37/* Table of values of atan(2^-i) in 0.32 format fractions of pi where pi = 0xffffffff / 2 */
38static const unsigned long atan_table[] = {
39 0x1fffffff, /* +0.785398163 (or pi/4) */
40 0x12e4051d, /* +0.463647609 */
41 0x09fb385b, /* +0.244978663 */
42 0x051111d4, /* +0.124354995 */
43 0x028b0d43, /* +0.062418810 */
44 0x0145d7e1, /* +0.031239833 */
45 0x00a2f61e, /* +0.015623729 */
46 0x00517c55, /* +0.007812341 */
47 0x0028be53, /* +0.003906230 */
48 0x00145f2e, /* +0.001953123 */
49 0x000a2f98, /* +0.000976562 */
50 0x000517cc, /* +0.000488281 */
51 0x00028be6, /* +0.000244141 */
52 0x000145f3, /* +0.000122070 */
53 0x0000a2f9, /* +0.000061035 */
54 0x0000517c, /* +0.000030518 */
55 0x000028be, /* +0.000015259 */
56 0x0000145f, /* +0.000007629 */
57 0x00000a2f, /* +0.000003815 */
58 0x00000517, /* +0.000001907 */
59 0x0000028b, /* +0.000000954 */
60 0x00000145, /* +0.000000477 */
61 0x000000a2, /* +0.000000238 */
62 0x00000051, /* +0.000000119 */
63 0x00000028, /* +0.000000060 */
64 0x00000014, /* +0.000000030 */
65 0x0000000a, /* +0.000000015 */
66 0x00000005, /* +0.000000007 */
67 0x00000002, /* +0.000000004 */
68 0x00000001, /* +0.000000002 */
69 0x00000000, /* +0.000000001 */
70 0x00000000, /* +0.000000000 */
71};
72
73/* Precalculated sine and cosine * 16384 (2^14) (fixed point 18.14) */
74static const short sin_table[91] =
75{
76 0, 285, 571, 857, 1142, 1427, 1712, 1996, 2280, 2563,
77 2845, 3126, 3406, 3685, 3963, 4240, 4516, 4790, 5062, 5334,
78 5603, 5871, 6137, 6401, 6663, 6924, 7182, 7438, 7691, 7943,
79 8191, 8438, 8682, 8923, 9161, 9397, 9630, 9860, 10086, 10310,
80 10531, 10748, 10963, 11173, 11381, 11585, 11785, 11982, 12175, 12365,
81 12550, 12732, 12910, 13084, 13254, 13420, 13582, 13740, 13894, 14043,
82 14188, 14329, 14466, 14598, 14725, 14848, 14967, 15081, 15190, 15295,
83 15395, 15491, 15582, 15668, 15749, 15825, 15897, 15964, 16025, 16082,
84 16135, 16182, 16224, 16261, 16294, 16321, 16344, 16361, 16374, 16381,
85 16384
86};
87
88/**
89 * Implements sin and cos using CORDIC rotation.
90 *
91 * @param phase has range from 0 to 0xffffffff, representing 0 and
92 * 2*pi respectively.
93 * @param cos return address for cos
94 * @return sin of phase, value is a signed value from LONG_MIN to LONG_MAX,
95 * representing -1 and 1 respectively.
96 */
97long fp_sincos(unsigned long phase, long *cos)
98{
99 int32_t x, x1, y, y1;
100 unsigned long z, z1;
101 int i;
102
103 /* Setup initial vector */
104 x = cordic_circular_gain;
105 y = 0;
106 z = phase;
107
108 /* The phase has to be somewhere between 0..pi for this to work right */
109 if (z < 0xffffffff / 4) {
110 /* z in first quadrant, z += pi/2 to correct */
111 x = -x;
112 z += 0xffffffff / 4;
113 } else if (z < 3 * (0xffffffff / 4)) {
114 /* z in third quadrant, z -= pi/2 to correct */
115 z -= 0xffffffff / 4;
116 } else {
117 /* z in fourth quadrant, z -= 3pi/2 to correct */
118 x = -x;
119 z -= 3 * (0xffffffff / 4);
120 }
121
122 /* Each iteration adds roughly 1-bit of extra precision */
123 for (i = 0; i < 31; i++) {
124 x1 = x >> i;
125 y1 = y >> i;
126 z1 = atan_table[i];
127
128 /* Decided which direction to rotate vector. Pivot point is pi/2 */
129 if (z >= 0xffffffff / 4) {
130 x -= y1;
131 y += x1;
132 z -= z1;
133 } else {
134 x += y1;
135 y -= x1;
136 z += z1;
137 }
138 }
139
140 if (cos)
141 *cos = x;
142
143 return y;
144}
145
146
147#if defined(PLUGIN) || defined(CODEC)
148/**
149 * Fixed point square root via Newton-Raphson.
150 * @param x square root argument.
151 * @param fracbits specifies number of fractional bits in argument.
152 * @return Square root of argument in same fixed point format as input.
153 *
154 * This routine has been modified to run longer for greater precision,
155 * but cuts calculation short if the answer is reached sooner.
156 */
157long fp_sqrt(long x, unsigned int fracbits)
158{
159 unsigned long xfp, b;
160 int n = 8; /* iteration limit (should terminate earlier) */
161
162 if (x <= 0)
163 return 0; /* no sqrt(neg), or just sqrt(0) = 0 */
164
165 /* Increase working precision by one bit */
166 xfp = x << 1;
167 fracbits++;
168
169 /* Get the midpoint between fracbits index and the highest bit index */
170 b = ((sizeof(xfp)*8-1) - __builtin_clzl(xfp) + fracbits) >> 1;
171 b = BIT_N(b);
172
173 do
174 {
175 unsigned long c = b;
176 b = (fp_div(xfp, b, fracbits) + b) >> 1;
177 if (c == b) break;
178 }
179 while (n-- > 0);
180
181 return b >> 1;
182}
183
184/* Accurate int sqrt with only elementary operations.
185 * Snagged from:
186 * http://www.devmaster.net/articles/fixed-point-optimizations/ */
187unsigned long isqrt(unsigned long x)
188{
189 /* Adding CLZ could optimize this further */
190 unsigned long g = 0;
191 int bshift = 15;
192 unsigned long b = 1ul << bshift;
193
194 do
195 {
196 unsigned long temp = (g + g + b) << bshift;
197
198 if (x > temp)
199 {
200 g += b;
201 x -= temp;
202 }
203
204 b >>= 1;
205 }
206 while (bshift--);
207
208 return g;
209}
210#endif /* PLUGIN or CODEC */
211
212
213#if defined(PLUGIN)
214/**
215 * Fixed point sinus using a lookup table
216 * don't forget to divide the result by 16384 to get the actual sinus value
217 * @param val sinus argument in degree
218 * @return sin(val)*16384
219 */
220long fp14_sin(int val)
221{
222 val = (val+360)%360;
223 if (val < 181)
224 {
225 if (val < 91)/* phase 0-90 degree */
226 return (long)sin_table[val];
227 else/* phase 91-180 degree */
228 return (long)sin_table[180-val];
229 }
230 else
231 {
232 if (val < 271)/* phase 181-270 degree */
233 return -(long)sin_table[val-180];
234 else/* phase 270-359 degree */
235 return -(long)sin_table[360-val];
236 }
237 return 0;
238}
239
240/**
241 * Fixed point cosinus using a lookup table
242 * don't forget to divide the result by 16384 to get the actual cosinus value
243 * @param val sinus argument in degree
244 * @return cos(val)*16384
245 */
246long fp14_cos(int val)
247{
248 val = (val+360)%360;
249 if (val < 181)
250 {
251 if (val < 91)/* phase 0-90 degree */
252 return (long)sin_table[90-val];
253 else/* phase 91-180 degree */
254 return -(long)sin_table[val-90];
255 }
256 else
257 {
258 if (val < 271)/* phase 181-270 degree */
259 return -(long)sin_table[270-val];
260 else/* phase 270-359 degree */
261 return (long)sin_table[val-270];
262 }
263 return 0;
264}
265
266/**
267 * Fixed-point natural log
268 * taken from http://www.quinapalus.com/efunc.html
269 * "The code assumes integers are at least 32 bits long. The (positive)
270 * argument and the result of the function are both expressed as fixed-point
271 * values with 16 fractional bits, although intermediates are kept with 28
272 * bits of precision to avoid loss of accuracy during shifts."
273 */
274long fp16_log(int x)
275{
276 int t;
277 int y = 0xa65af;
278
279 if (x < 0x00008000) x <<=16, y -= 0xb1721;
280 if (x < 0x00800000) x <<= 8, y -= 0x58b91;
281 if (x < 0x08000000) x <<= 4, y -= 0x2c5c8;
282 if (x < 0x20000000) x <<= 2, y -= 0x162e4;
283 if (x < 0x40000000) x <<= 1, y -= 0x0b172;
284 t = x + (x >> 1); if ((t & 0x80000000) == 0) x = t, y -= 0x067cd;
285 t = x + (x >> 2); if ((t & 0x80000000) == 0) x = t, y -= 0x03920;
286 t = x + (x >> 3); if ((t & 0x80000000) == 0) x = t, y -= 0x01e27;
287 t = x + (x >> 4); if ((t & 0x80000000) == 0) x = t, y -= 0x00f85;
288 t = x + (x >> 5); if ((t & 0x80000000) == 0) x = t, y -= 0x007e1;
289 t = x + (x >> 6); if ((t & 0x80000000) == 0) x = t, y -= 0x003f8;
290 t = x + (x >> 7); if ((t & 0x80000000) == 0) x = t, y -= 0x001fe;
291 x = 0x80000000 - x;
292 y -= x >> 15;
293
294 return y;
295}
296
297/**
298 * Fixed-point exponential
299 * taken from http://www.quinapalus.com/efunc.html
300 * "The code assumes integers are at least 32 bits long. The (non-negative)
301 * argument and the result of the function are both expressed as fixed-point
302 * values with 16 fractional bits. Notice that after 11 steps of the
303 * algorithm the constants involved become such that the code is simply
304 * doing a multiplication: this is explained in the note below.
305 * The extension to negative arguments is left as an exercise."
306 */
307long fp16_exp(int x)
308{
309 int t;
310 int y = 0x00010000;
311
312 if (x < 0) x += 0xb1721, y >>= 16;
313 t = x - 0x58b91; if (t >= 0) x = t, y <<= 8;
314 t = x - 0x2c5c8; if (t >= 0) x = t, y <<= 4;
315 t = x - 0x162e4; if (t >= 0) x = t, y <<= 2;
316 t = x - 0x0b172; if (t >= 0) x = t, y <<= 1;
317 t = x - 0x067cd; if (t >= 0) x = t, y += y >> 1;
318 t = x - 0x03920; if (t >= 0) x = t, y += y >> 2;
319 t = x - 0x01e27; if (t >= 0) x = t, y += y >> 3;
320 t = x - 0x00f85; if (t >= 0) x = t, y += y >> 4;
321 t = x - 0x007e1; if (t >= 0) x = t, y += y >> 5;
322 t = x - 0x003f8; if (t >= 0) x = t, y += y >> 6;
323 t = x - 0x001fe; if (t >= 0) x = t, y += y >> 7;
324 y += ((y >> 8) * x) >> 8;
325
326 return y;
327}
328#endif /* PLUGIN */
329
330
331#if (!defined(PLUGIN) && !defined(CODEC))
332/** MODIFIED FROM replaygain.c */
333
334#define FP_MUL_FRAC(x, y) fp_mul(x, y, fracbits)
335#define FP_DIV_FRAC(x, y) fp_div(x, y, fracbits)
336
337/* constants in fixed point format, 28 fractional bits */
338#define FP28_LN2 (186065279L) /* ln(2) */
339#define FP28_LN2_INV (387270501L) /* 1/ln(2) */
340#define FP28_EXP_ZERO (44739243L) /* 1/6 */
341#define FP28_EXP_ONE (-745654L) /* -1/360 */
342#define FP28_EXP_TWO (12428L) /* 1/21600 */
343#define FP28_LN10 (618095479L) /* ln(10) */
344#define FP28_LOG10OF2 (80807124L) /* log10(2) */
345
346#define TOL_BITS 2 /* log calculation tolerance */
347
348
349/* The fpexp10 fixed point math routine is based
350 * on oMathFP by Dan Carter (http://orbisstudios.com).
351 */
352
353/** FIXED POINT EXP10
354 * Return 10^x as FP integer. Argument is FP integer.
355 */
356long fp_exp10(long x, unsigned int fracbits)
357{
358 long k;
359 long z;
360 long R;
361 long xp;
362
363 /* scale constants */
364 const long fp_one = (1 << fracbits);
365 const long fp_half = (1 << (fracbits - 1));
366 const long fp_two = (2 << fracbits);
367 const long fp_mask = (fp_one - 1);
368 const long fp_ln2_inv = (FP28_LN2_INV >> (28 - fracbits));
369 const long fp_ln2 = (FP28_LN2 >> (28 - fracbits));
370 const long fp_ln10 = (FP28_LN10 >> (28 - fracbits));
371 const long fp_exp_zero = (FP28_EXP_ZERO >> (28 - fracbits));
372 const long fp_exp_one = (FP28_EXP_ONE >> (28 - fracbits));
373 const long fp_exp_two = (FP28_EXP_TWO >> (28 - fracbits));
374
375 /* exp(0) = 1 */
376 if (x == 0)
377 {
378 return fp_one;
379 }
380
381 /* convert from base 10 to base e */
382 x = FP_MUL_FRAC(x, fp_ln10);
383
384 /* calculate exp(x) */
385 k = (FP_MUL_FRAC(abs(x), fp_ln2_inv) + fp_half) & ~fp_mask;
386
387 if (x < 0)
388 {
389 k = -k;
390 }
391
392 x -= FP_MUL_FRAC(k, fp_ln2);
393 z = FP_MUL_FRAC(x, x);
394 R = fp_two + FP_MUL_FRAC(z, fp_exp_zero + FP_MUL_FRAC(z, fp_exp_one
395 + FP_MUL_FRAC(z, fp_exp_two)));
396 xp = fp_one + FP_DIV_FRAC(FP_MUL_FRAC(fp_two, x), R - x);
397
398 if (k < 0)
399 {
400 k = fp_one >> (-k >> fracbits);
401 }
402 else
403 {
404 k = fp_one << (k >> fracbits);
405 }
406
407 return FP_MUL_FRAC(k, xp);
408}
409
410
411#if 0 /* useful code, but not currently used */
412/** FIXED POINT LOG10
413 * Return log10(x) as FP integer. Argument is FP integer.
414 */
415static long fp_log10(long n, unsigned int fracbits)
416{
417 /* Calculate log2 of argument */
418
419 long log2, frac;
420 const long fp_one = (1 << fracbits);
421 const long fp_two = (2 << fracbits);
422 const long tolerance = (1 << ((fracbits / 2) + 2));
423
424 if (n <=0) return FP_NEGINF;
425 log2 = 0;
426
427 /* integer part */
428 while (n < fp_one)
429 {
430 log2 -= fp_one;
431 n <<= 1;
432 }
433 while (n >= fp_two)
434 {
435 log2 += fp_one;
436 n >>= 1;
437 }
438
439 /* fractional part */
440 frac = fp_one;
441 while (frac > tolerance)
442 {
443 frac >>= 1;
444 n = FP_MUL_FRAC(n, n);
445 if (n >= fp_two)
446 {
447 n >>= 1;
448 log2 += frac;
449 }
450 }
451
452 /* convert log2 to log10 */
453 return FP_MUL_FRAC(log2, (FP28_LOG10OF2 >> (28 - fracbits)));
454}
455
456
457/** CONVERT FACTOR TO DECIBELS */
458long fp_decibels(unsigned long factor, unsigned int fracbits)
459{
460 /* decibels = 20 * log10(factor) */
461 return FP_MUL_FRAC((20L << fracbits), fp_log10(factor, fracbits));
462}
463#endif /* unused code */
464
465
466/** CONVERT DECIBELS TO FACTOR */
467long fp_factor(long decibels, unsigned int fracbits)
468{
469 /* factor = 10 ^ (decibels / 20) */
470 return fp_exp10(FP_DIV_FRAC(decibels, (20L << fracbits)), fracbits);
471}
472#endif /* !PLUGIN and !CODEC */
generated by cgit v1.2.3 (git 2.39.1) at 2024-10-05 22:31:42 +0000