new plugin: FS#10559 - lrcplayer: a plugin to view .lrc file.
[kugel-rb.git] / apps / fixedpoint.c
blobfb89a8d30f6ed88a0d044c6a1543763066c51490
1 /***************************************************************************
2 * __________ __ ___.
3 * Open \______ \ ____ ____ | | _\_ |__ _______ ___
4 * Source | _// _ \_/ ___\| |/ /| __ \ / _ \ \/ /
5 * Jukebox | | ( <_> ) \___| < | \_\ ( <_> > < <
6 * Firmware |____|_ /\____/ \___ >__|_ \|___ /\____/__/\_ \
7 * \/ \/ \/ \/ \/
8 * $Id$
10 * Copyright (C) 2006 Jens Arnold
12 * Fixed point library for plugins
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.
19 * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY OF ANY
20 * KIND, either express or implied.
22 ****************************************************************************/
24 #include "fixedpoint.h"
25 #include <stdlib.h>
26 #include <stdbool.h>
27 #include <inttypes.h>
29 #ifndef BIT_N
30 #define BIT_N(n) (1U << (n))
31 #endif
33 /** TAKEN FROM ORIGINAL fixedpoint.h */
34 /* Inverse gain of circular cordic rotation in s0.31 format. */
35 static const long cordic_circular_gain = 0xb2458939; /* 0.607252929 */
37 /* Table of values of atan(2^-i) in 0.32 format fractions of pi where pi = 0xffffffff / 2 */
38 static 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 */
73 /* Precalculated sine and cosine * 16384 (2^14) (fixed point 18.14) */
74 static const short sin_table[91] =
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
88 /**
89 * Implements sin and cos using CORDIC rotation.
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.
97 long fp_sincos(unsigned long phase, long *cos)
99 int32_t x, x1, y, y1;
100 unsigned long z, z1;
101 int i;
103 /* Setup initial vector */
104 x = cordic_circular_gain;
105 y = 0;
106 z = phase;
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);
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];
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;
140 if (cos)
141 *cos = x;
143 return y;
147 #if defined(PLUGIN) || defined(CODEC)
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.
154 * This routine has been modified to run longer for greater precision,
155 * but cuts calculation short if the answer is reached sooner. In
156 * general, the closer x is to 1, the quicker the calculation.
158 long fp_sqrt(long x, unsigned int fracbits)
160 long b = x/2 + BIT_N(fracbits); /* initial approximation */
161 long c;
162 unsigned n;
163 const unsigned iterations = 8;
165 for (n = 0; n < iterations; ++n)
167 c = fp_div(x, b, fracbits);
168 if (c == b) break;
169 b = (b + c)/2;
172 return b;
175 /* Accurate int sqrt with only elementary operations. (the above
176 * routine fails badly without enough iterations, more iterations
177 * than this requires -- [give that one a FIXME]).
178 * Snagged from:
179 * http://www.devmaster.net/articles/fixed-point-optimizations/ */
180 unsigned long isqrt(unsigned long x)
182 /* Adding CLZ could optimize this further */
183 unsigned long g = 0;
184 int bshift = 15;
185 unsigned long b = 1ul << bshift;
189 unsigned long temp = (g + g + b) << bshift;
191 if (x > temp)
193 g += b;
194 x -= temp;
197 b >>= 1;
199 while (bshift--);
201 return g;
203 #endif /* PLUGIN or CODEC */
206 #if defined(PLUGIN)
208 * Fixed point sinus using a lookup table
209 * don't forget to divide the result by 16384 to get the actual sinus value
210 * @param val sinus argument in degree
211 * @return sin(val)*16384
213 long fp14_sin(int val)
215 val = (val+360)%360;
216 if (val < 181)
218 if (val < 91)/* phase 0-90 degree */
219 return (long)sin_table[val];
220 else/* phase 91-180 degree */
221 return (long)sin_table[180-val];
223 else
225 if (val < 271)/* phase 181-270 degree */
226 return -(long)sin_table[val-180];
227 else/* phase 270-359 degree */
228 return -(long)sin_table[360-val];
230 return 0;
234 * Fixed point cosinus using a lookup table
235 * don't forget to divide the result by 16384 to get the actual cosinus value
236 * @param val sinus argument in degree
237 * @return cos(val)*16384
239 long fp14_cos(int val)
241 val = (val+360)%360;
242 if (val < 181)
244 if (val < 91)/* phase 0-90 degree */
245 return (long)sin_table[90-val];
246 else/* phase 91-180 degree */
247 return -(long)sin_table[val-90];
249 else
251 if (val < 271)/* phase 181-270 degree */
252 return -(long)sin_table[270-val];
253 else/* phase 270-359 degree */
254 return (long)sin_table[val-270];
256 return 0;
260 * Fixed-point natural log
261 * taken from http://www.quinapalus.com/efunc.html
262 * "The code assumes integers are at least 32 bits long. The (positive)
263 * argument and the result of the function are both expressed as fixed-point
264 * values with 16 fractional bits, although intermediates are kept with 28
265 * bits of precision to avoid loss of accuracy during shifts."
268 long fp16_log(int x) {
269 long t,y;
271 y=0xa65af;
272 if(x<0x00008000) x<<=16, y-=0xb1721;
273 if(x<0x00800000) x<<= 8, y-=0x58b91;
274 if(x<0x08000000) x<<= 4, y-=0x2c5c8;
275 if(x<0x20000000) x<<= 2, y-=0x162e4;
276 if(x<0x40000000) x<<= 1, y-=0x0b172;
277 t=x+(x>>1); if((t&0x80000000)==0) x=t,y-=0x067cd;
278 t=x+(x>>2); if((t&0x80000000)==0) x=t,y-=0x03920;
279 t=x+(x>>3); if((t&0x80000000)==0) x=t,y-=0x01e27;
280 t=x+(x>>4); if((t&0x80000000)==0) x=t,y-=0x00f85;
281 t=x+(x>>5); if((t&0x80000000)==0) x=t,y-=0x007e1;
282 t=x+(x>>6); if((t&0x80000000)==0) x=t,y-=0x003f8;
283 t=x+(x>>7); if((t&0x80000000)==0) x=t,y-=0x001fe;
284 x=0x80000000-x;
285 y-=x>>15;
286 return y;
290 * Fixed-point exponential
291 * taken from http://www.quinapalus.com/efunc.html
292 * "The code assumes integers are at least 32 bits long. The (non-negative)
293 * argument and the result of the function are both expressed as fixed-point
294 * values with 16 fractional bits. Notice that after 11 steps of the
295 * algorithm the constants involved become such that the code is simply
296 * doing a multiplication: this is explained in the note below.
297 * The extension to negative arguments is left as an exercise."
299 long fp16_exp(int x)
301 int t,y;
303 y=0x00010000;
304 t=x-0x58b91; if(t>=0) x=t,y<<=8;
305 t=x-0x2c5c8; if(t>=0) x=t,y<<=4;
306 t=x-0x162e4; if(t>=0) x=t,y<<=2;
307 t=x-0x0b172; if(t>=0) x=t,y<<=1;
308 t=x-0x067cd; if(t>=0) x=t,y+=y>>1;
309 t=x-0x03920; if(t>=0) x=t,y+=y>>2;
310 t=x-0x01e27; if(t>=0) x=t,y+=y>>3;
311 t=x-0x00f85; if(t>=0) x=t,y+=y>>4;
312 t=x-0x007e1; if(t>=0) x=t,y+=y>>5;
313 t=x-0x003f8; if(t>=0) x=t,y+=y>>6;
314 t=x-0x001fe; if(t>=0) x=t,y+=y>>7;
315 if(x&0x100) y+=y>>8;
316 if(x&0x080) y+=y>>9;
317 if(x&0x040) y+=y>>10;
318 if(x&0x020) y+=y>>11;
319 if(x&0x010) y+=y>>12;
320 if(x&0x008) y+=y>>13;
321 if(x&0x004) y+=y>>14;
322 if(x&0x002) y+=y>>15;
323 if(x&0x001) y+=y>>16;
324 return y;
326 #endif /* PLUGIN */
329 #if (!defined(PLUGIN) && !defined(CODEC))
330 /** MODIFIED FROM replaygain.c */
332 #define FP_MUL_FRAC(x, y) fp_mul(x, y, fracbits)
333 #define FP_DIV_FRAC(x, y) fp_div(x, y, fracbits)
335 /* constants in fixed point format, 28 fractional bits */
336 #define FP28_LN2 (186065279L) /* ln(2) */
337 #define FP28_LN2_INV (387270501L) /* 1/ln(2) */
338 #define FP28_EXP_ZERO (44739243L) /* 1/6 */
339 #define FP28_EXP_ONE (-745654L) /* -1/360 */
340 #define FP28_EXP_TWO (12428L) /* 1/21600 */
341 #define FP28_LN10 (618095479L) /* ln(10) */
342 #define FP28_LOG10OF2 (80807124L) /* log10(2) */
344 #define TOL_BITS 2 /* log calculation tolerance */
347 /* The fpexp10 fixed point math routine is based
348 * on oMathFP by Dan Carter (http://orbisstudios.com).
351 /** FIXED POINT EXP10
352 * Return 10^x as FP integer. Argument is FP integer.
354 static long fp_exp10(long x, unsigned int fracbits)
356 long k;
357 long z;
358 long R;
359 long xp;
361 /* scale constants */
362 const long fp_one = (1 << fracbits);
363 const long fp_half = (1 << (fracbits - 1));
364 const long fp_two = (2 << fracbits);
365 const long fp_mask = (fp_one - 1);
366 const long fp_ln2_inv = (FP28_LN2_INV >> (28 - fracbits));
367 const long fp_ln2 = (FP28_LN2 >> (28 - fracbits));
368 const long fp_ln10 = (FP28_LN10 >> (28 - fracbits));
369 const long fp_exp_zero = (FP28_EXP_ZERO >> (28 - fracbits));
370 const long fp_exp_one = (FP28_EXP_ONE >> (28 - fracbits));
371 const long fp_exp_two = (FP28_EXP_TWO >> (28 - fracbits));
373 /* exp(0) = 1 */
374 if (x == 0)
376 return fp_one;
379 /* convert from base 10 to base e */
380 x = FP_MUL_FRAC(x, fp_ln10);
382 /* calculate exp(x) */
383 k = (FP_MUL_FRAC(abs(x), fp_ln2_inv) + fp_half) & ~fp_mask;
385 if (x < 0)
387 k = -k;
390 x -= FP_MUL_FRAC(k, fp_ln2);
391 z = FP_MUL_FRAC(x, x);
392 R = fp_two + FP_MUL_FRAC(z, fp_exp_zero + FP_MUL_FRAC(z, fp_exp_one
393 + FP_MUL_FRAC(z, fp_exp_two)));
394 xp = fp_one + FP_DIV_FRAC(FP_MUL_FRAC(fp_two, x), R - x);
396 if (k < 0)
398 k = fp_one >> (-k >> fracbits);
400 else
402 k = fp_one << (k >> fracbits);
405 return FP_MUL_FRAC(k, xp);
409 #if 0 /* useful code, but not currently used */
410 /** FIXED POINT LOG10
411 * Return log10(x) as FP integer. Argument is FP integer.
413 static long fp_log10(long n, unsigned int fracbits)
415 /* Calculate log2 of argument */
417 long log2, frac;
418 const long fp_one = (1 << fracbits);
419 const long fp_two = (2 << fracbits);
420 const long tolerance = (1 << ((fracbits / 2) + 2));
422 if (n <=0) return FP_NEGINF;
423 log2 = 0;
425 /* integer part */
426 while (n < fp_one)
428 log2 -= fp_one;
429 n <<= 1;
431 while (n >= fp_two)
433 log2 += fp_one;
434 n >>= 1;
437 /* fractional part */
438 frac = fp_one;
439 while (frac > tolerance)
441 frac >>= 1;
442 n = FP_MUL_FRAC(n, n);
443 if (n >= fp_two)
445 n >>= 1;
446 log2 += frac;
450 /* convert log2 to log10 */
451 return FP_MUL_FRAC(log2, (FP28_LOG10OF2 >> (28 - fracbits)));
455 /** CONVERT FACTOR TO DECIBELS */
456 long fp_decibels(unsigned long factor, unsigned int fracbits)
458 /* decibels = 20 * log10(factor) */
459 return FP_MUL_FRAC((20L << fracbits), fp_log10(factor, fracbits));
461 #endif /* unused code */
464 /** CONVERT DECIBELS TO FACTOR */
465 long fp_factor(long decibels, unsigned int fracbits)
467 /* factor = 10 ^ (decibels / 20) */
468 return fp_exp10(FP_DIV_FRAC(decibels, (20L << fracbits)), fracbits);
470 #endif /* !PLUGIN and !CODEC */