Bump version numbers for 3.13
[maemo-rb.git] / apps / fixedpoint.c
bloba8f606a5db7a64a3cc79f86fcb8fd11824725c9b
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.
157 long fp_sqrt(long x, unsigned int fracbits)
159 unsigned long xfp, b;
160 int n = 8; /* iteration limit (should terminate earlier) */
162 if (x <= 0)
163 return 0; /* no sqrt(neg), or just sqrt(0) = 0 */
165 /* Increase working precision by one bit */
166 xfp = x << 1;
167 fracbits++;
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);
175 unsigned long c = b;
176 b = (fp_div(xfp, b, fracbits) + b) >> 1;
177 if (c == b) break;
179 while (n-- > 0);
181 return b >> 1;
184 /* Accurate int sqrt with only elementary operations.
185 * Snagged from:
186 * http://www.devmaster.net/articles/fixed-point-optimizations/ */
187 unsigned long isqrt(unsigned long x)
189 /* Adding CLZ could optimize this further */
190 unsigned long g = 0;
191 int bshift = 15;
192 unsigned long b = 1ul << bshift;
196 unsigned long temp = (g + g + b) << bshift;
198 if (x > temp)
200 g += b;
201 x -= temp;
204 b >>= 1;
206 while (bshift--);
208 return g;
210 #endif /* PLUGIN or CODEC */
213 #if defined(PLUGIN)
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
220 long fp14_sin(int val)
222 val = (val+360)%360;
223 if (val < 181)
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];
230 else
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];
237 return 0;
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
246 long fp14_cos(int val)
248 val = (val+360)%360;
249 if (val < 181)
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];
256 else
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];
263 return 0;
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."
275 long fp16_log(int x) {
276 long t,y;
278 y=0xa65af;
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 return y;
297 * Fixed-point exponential
298 * taken from http://www.quinapalus.com/efunc.html
299 * "The code assumes integers are at least 32 bits long. The (non-negative)
300 * argument and the result of the function are both expressed as fixed-point
301 * values with 16 fractional bits. Notice that after 11 steps of the
302 * algorithm the constants involved become such that the code is simply
303 * doing a multiplication: this is explained in the note below.
304 * The extension to negative arguments is left as an exercise."
306 long fp16_exp(int x)
308 int t,y;
310 y=0x00010000;
311 t=x-0x58b91; if(t>=0) x=t,y<<=8;
312 t=x-0x2c5c8; if(t>=0) x=t,y<<=4;
313 t=x-0x162e4; if(t>=0) x=t,y<<=2;
314 t=x-0x0b172; if(t>=0) x=t,y<<=1;
315 t=x-0x067cd; if(t>=0) x=t,y+=y>>1;
316 t=x-0x03920; if(t>=0) x=t,y+=y>>2;
317 t=x-0x01e27; if(t>=0) x=t,y+=y>>3;
318 t=x-0x00f85; if(t>=0) x=t,y+=y>>4;
319 t=x-0x007e1; if(t>=0) x=t,y+=y>>5;
320 t=x-0x003f8; if(t>=0) x=t,y+=y>>6;
321 t=x-0x001fe; if(t>=0) x=t,y+=y>>7;
322 if(x&0x100) y+=y>>8;
323 if(x&0x080) y+=y>>9;
324 if(x&0x040) y+=y>>10;
325 if(x&0x020) y+=y>>11;
326 if(x&0x010) y+=y>>12;
327 if(x&0x008) y+=y>>13;
328 if(x&0x004) y+=y>>14;
329 if(x&0x002) y+=y>>15;
330 if(x&0x001) y+=y>>16;
331 return y;
333 #endif /* PLUGIN */
336 #if (!defined(PLUGIN) && !defined(CODEC))
337 /** MODIFIED FROM replaygain.c */
339 #define FP_MUL_FRAC(x, y) fp_mul(x, y, fracbits)
340 #define FP_DIV_FRAC(x, y) fp_div(x, y, fracbits)
342 /* constants in fixed point format, 28 fractional bits */
343 #define FP28_LN2 (186065279L) /* ln(2) */
344 #define FP28_LN2_INV (387270501L) /* 1/ln(2) */
345 #define FP28_EXP_ZERO (44739243L) /* 1/6 */
346 #define FP28_EXP_ONE (-745654L) /* -1/360 */
347 #define FP28_EXP_TWO (12428L) /* 1/21600 */
348 #define FP28_LN10 (618095479L) /* ln(10) */
349 #define FP28_LOG10OF2 (80807124L) /* log10(2) */
351 #define TOL_BITS 2 /* log calculation tolerance */
354 /* The fpexp10 fixed point math routine is based
355 * on oMathFP by Dan Carter (http://orbisstudios.com).
358 /** FIXED POINT EXP10
359 * Return 10^x as FP integer. Argument is FP integer.
361 static long fp_exp10(long x, unsigned int fracbits)
363 long k;
364 long z;
365 long R;
366 long xp;
368 /* scale constants */
369 const long fp_one = (1 << fracbits);
370 const long fp_half = (1 << (fracbits - 1));
371 const long fp_two = (2 << fracbits);
372 const long fp_mask = (fp_one - 1);
373 const long fp_ln2_inv = (FP28_LN2_INV >> (28 - fracbits));
374 const long fp_ln2 = (FP28_LN2 >> (28 - fracbits));
375 const long fp_ln10 = (FP28_LN10 >> (28 - fracbits));
376 const long fp_exp_zero = (FP28_EXP_ZERO >> (28 - fracbits));
377 const long fp_exp_one = (FP28_EXP_ONE >> (28 - fracbits));
378 const long fp_exp_two = (FP28_EXP_TWO >> (28 - fracbits));
380 /* exp(0) = 1 */
381 if (x == 0)
383 return fp_one;
386 /* convert from base 10 to base e */
387 x = FP_MUL_FRAC(x, fp_ln10);
389 /* calculate exp(x) */
390 k = (FP_MUL_FRAC(abs(x), fp_ln2_inv) + fp_half) & ~fp_mask;
392 if (x < 0)
394 k = -k;
397 x -= FP_MUL_FRAC(k, fp_ln2);
398 z = FP_MUL_FRAC(x, x);
399 R = fp_two + FP_MUL_FRAC(z, fp_exp_zero + FP_MUL_FRAC(z, fp_exp_one
400 + FP_MUL_FRAC(z, fp_exp_two)));
401 xp = fp_one + FP_DIV_FRAC(FP_MUL_FRAC(fp_two, x), R - x);
403 if (k < 0)
405 k = fp_one >> (-k >> fracbits);
407 else
409 k = fp_one << (k >> fracbits);
412 return FP_MUL_FRAC(k, xp);
416 #if 0 /* useful code, but not currently used */
417 /** FIXED POINT LOG10
418 * Return log10(x) as FP integer. Argument is FP integer.
420 static long fp_log10(long n, unsigned int fracbits)
422 /* Calculate log2 of argument */
424 long log2, frac;
425 const long fp_one = (1 << fracbits);
426 const long fp_two = (2 << fracbits);
427 const long tolerance = (1 << ((fracbits / 2) + 2));
429 if (n <=0) return FP_NEGINF;
430 log2 = 0;
432 /* integer part */
433 while (n < fp_one)
435 log2 -= fp_one;
436 n <<= 1;
438 while (n >= fp_two)
440 log2 += fp_one;
441 n >>= 1;
444 /* fractional part */
445 frac = fp_one;
446 while (frac > tolerance)
448 frac >>= 1;
449 n = FP_MUL_FRAC(n, n);
450 if (n >= fp_two)
452 n >>= 1;
453 log2 += frac;
457 /* convert log2 to log10 */
458 return FP_MUL_FRAC(log2, (FP28_LOG10OF2 >> (28 - fracbits)));
462 /** CONVERT FACTOR TO DECIBELS */
463 long fp_decibels(unsigned long factor, unsigned int fracbits)
465 /* decibels = 20 * log10(factor) */
466 return FP_MUL_FRAC((20L << fracbits), fp_log10(factor, fracbits));
468 #endif /* unused code */
471 /** CONVERT DECIBELS TO FACTOR */
472 long fp_factor(long decibels, unsigned int fracbits)
474 /* factor = 10 ^ (decibels / 20) */
475 return fp_exp10(FP_DIV_FRAC(decibels, (20L << fracbits)), fracbits);
477 #endif /* !PLUGIN and !CODEC */