1 /***************************************************************************
3 * Open \______ \ ____ ____ | | _\_ |__ _______ ___
4 * Source | _// _ \_/ ___\| |/ /| __ \ / _ \ \/ /
5 * Jukebox | | ( <_> ) \___| < | \_\ ( <_> > < <
6 * Firmware |____|_ /\____/ \___ >__|_ \|___ /\____/__/\_ \
8 * $Id: fixedpoint.c -1 $
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"
29 #define BIT_N(n) (1U << (n))
32 /** TAKEN FROM ORIGINAL fixedpoint.h */
33 /* Inverse gain of circular cordic rotation in s0.31 format. */
34 static const long cordic_circular_gain
= 0xb2458939; /* 0.607252929 */
36 /* Table of values of atan(2^-i) in 0.32 format fractions of pi where pi = 0xffffffff / 2 */
37 static const unsigned long atan_table
[] = {
38 0x1fffffff, /* +0.785398163 (or pi/4) */
39 0x12e4051d, /* +0.463647609 */
40 0x09fb385b, /* +0.244978663 */
41 0x051111d4, /* +0.124354995 */
42 0x028b0d43, /* +0.062418810 */
43 0x0145d7e1, /* +0.031239833 */
44 0x00a2f61e, /* +0.015623729 */
45 0x00517c55, /* +0.007812341 */
46 0x0028be53, /* +0.003906230 */
47 0x00145f2e, /* +0.001953123 */
48 0x000a2f98, /* +0.000976562 */
49 0x000517cc, /* +0.000488281 */
50 0x00028be6, /* +0.000244141 */
51 0x000145f3, /* +0.000122070 */
52 0x0000a2f9, /* +0.000061035 */
53 0x0000517c, /* +0.000030518 */
54 0x000028be, /* +0.000015259 */
55 0x0000145f, /* +0.000007629 */
56 0x00000a2f, /* +0.000003815 */
57 0x00000517, /* +0.000001907 */
58 0x0000028b, /* +0.000000954 */
59 0x00000145, /* +0.000000477 */
60 0x000000a2, /* +0.000000238 */
61 0x00000051, /* +0.000000119 */
62 0x00000028, /* +0.000000060 */
63 0x00000014, /* +0.000000030 */
64 0x0000000a, /* +0.000000015 */
65 0x00000005, /* +0.000000007 */
66 0x00000002, /* +0.000000004 */
67 0x00000001, /* +0.000000002 */
68 0x00000000, /* +0.000000001 */
69 0x00000000, /* +0.000000000 */
72 /* Precalculated sine and cosine * 16384 (2^14) (fixed point 18.14) */
73 static const short sin_table
[91] =
75 0, 285, 571, 857, 1142, 1427, 1712, 1996, 2280, 2563,
76 2845, 3126, 3406, 3685, 3963, 4240, 4516, 4790, 5062, 5334,
77 5603, 5871, 6137, 6401, 6663, 6924, 7182, 7438, 7691, 7943,
78 8191, 8438, 8682, 8923, 9161, 9397, 9630, 9860, 10086, 10310,
79 10531, 10748, 10963, 11173, 11381, 11585, 11785, 11982, 12175, 12365,
80 12550, 12732, 12910, 13084, 13254, 13420, 13582, 13740, 13894, 14043,
81 14188, 14329, 14466, 14598, 14725, 14848, 14967, 15081, 15190, 15295,
82 15395, 15491, 15582, 15668, 15749, 15825, 15897, 15964, 16025, 16082,
83 16135, 16182, 16224, 16261, 16294, 16321, 16344, 16361, 16374, 16381,
88 * Implements sin and cos using CORDIC rotation.
90 * @param phase has range from 0 to 0xffffffff, representing 0 and
92 * @param cos return address for cos
93 * @return sin of phase, value is a signed value from LONG_MIN to LONG_MAX,
94 * representing -1 and 1 respectively.
96 long fsincos(unsigned long phase
, long *cos
)
102 /* Setup initial vector */
103 x
= cordic_circular_gain
;
107 /* The phase has to be somewhere between 0..pi for this to work right */
108 if (z
< 0xffffffff / 4) {
109 /* z in first quadrant, z += pi/2 to correct */
112 } else if (z
< 3 * (0xffffffff / 4)) {
113 /* z in third quadrant, z -= pi/2 to correct */
116 /* z in fourth quadrant, z -= 3pi/2 to correct */
118 z
-= 3 * (0xffffffff / 4);
121 /* Each iteration adds roughly 1-bit of extra precision */
122 for (i
= 0; i
< 31; i
++) {
127 /* Decided which direction to rotate vector. Pivot point is pi/2 */
128 if (z
>= 0xffffffff / 4) {
146 * Fixed point square root via Newton-Raphson.
147 * @param x square root argument.
148 * @param fracbits specifies number of fractional bits in argument.
149 * @return Square root of argument in same fixed point format as input.
151 * This routine has been modified to run longer for greater precision,
152 * but cuts calculation short if the answer is reached sooner. In
153 * general, the closer x is to 1, the quicker the calculation.
155 long fsqrt(long x
, unsigned int fracbits
)
157 long b
= x
/2 + BIT_N(fracbits
); /* initial approximation */
160 const unsigned iterations
= 8;
162 for (n
= 0; n
< iterations
; ++n
)
164 c
= DIV64(x
, b
, fracbits
);
173 * Fixed point sinus using a lookup table
174 * don't forget to divide the result by 16384 to get the actual sinus value
175 * @param val sinus argument in degree
176 * @return sin(val)*16384
178 long sin_int(int val
)
183 if (val
< 91)/* phase 0-90 degree */
184 return (long)sin_table
[val
];
185 else/* phase 91-180 degree */
186 return (long)sin_table
[180-val
];
190 if (val
< 271)/* phase 181-270 degree */
191 return -(long)sin_table
[val
-180];
192 else/* phase 270-359 degree */
193 return -(long)sin_table
[360-val
];
199 * Fixed point cosinus using a lookup table
200 * don't forget to divide the result by 16384 to get the actual cosinus value
201 * @param val sinus argument in degree
202 * @return cos(val)*16384
204 long cos_int(int val
)
209 if (val
< 91)/* phase 0-90 degree */
210 return (long)sin_table
[90-val
];
211 else/* phase 91-180 degree */
212 return -(long)sin_table
[val
-90];
216 if (val
< 271)/* phase 181-270 degree */
217 return -(long)sin_table
[270-val
];
218 else/* phase 270-359 degree */
219 return (long)sin_table
[val
-270];
225 * Fixed-point natural log
226 * taken from http://www.quinapalus.com/efunc.html
227 * "The code assumes integers are at least 32 bits long. The (positive)
228 * argument and the result of the function are both expressed as fixed-point
229 * values with 16 fractional bits, although intermediates are kept with 28
230 * bits of precision to avoid loss of accuracy during shifts."
237 if(x
<0x00008000) x
<<=16, y
-=0xb1721;
238 if(x
<0x00800000) x
<<= 8, y
-=0x58b91;
239 if(x
<0x08000000) x
<<= 4, y
-=0x2c5c8;
240 if(x
<0x20000000) x
<<= 2, y
-=0x162e4;
241 if(x
<0x40000000) x
<<= 1, y
-=0x0b172;
242 t
=x
+(x
>>1); if((t
&0x80000000)==0) x
=t
,y
-=0x067cd;
243 t
=x
+(x
>>2); if((t
&0x80000000)==0) x
=t
,y
-=0x03920;
244 t
=x
+(x
>>3); if((t
&0x80000000)==0) x
=t
,y
-=0x01e27;
245 t
=x
+(x
>>4); if((t
&0x80000000)==0) x
=t
,y
-=0x00f85;
246 t
=x
+(x
>>5); if((t
&0x80000000)==0) x
=t
,y
-=0x007e1;
247 t
=x
+(x
>>6); if((t
&0x80000000)==0) x
=t
,y
-=0x003f8;
248 t
=x
+(x
>>7); if((t
&0x80000000)==0) x
=t
,y
-=0x001fe;
254 /** MODIFIED FROM replaygain.c */
255 /* These math routines have 64-bit internal precision to avoid overflows.
256 * Arguments and return values are 32-bit (long) precision.
259 #define FP_MUL64(x, y) (((x) * (y)) >> (fracbits))
260 #define FP_DIV64(x, y) (((x) << (fracbits)) / (y))
262 static long long fp_exp10(long long x
, unsigned int fracbits
);
263 static long long fp_log10(long long n
, unsigned int fracbits
);
265 /* constants in fixed point format, 28 fractional bits */
266 #define FP28_LN2 (186065279LL) /* ln(2) */
267 #define FP28_LN2_INV (387270501LL) /* 1/ln(2) */
268 #define FP28_EXP_ZERO (44739243LL) /* 1/6 */
269 #define FP28_EXP_ONE (-745654LL) /* -1/360 */
270 #define FP28_EXP_TWO (12428LL) /* 1/21600 */
271 #define FP28_LN10 (618095479LL) /* ln(10) */
272 #define FP28_LOG10OF2 (80807124LL) /* log10(2) */
274 #define TOL_BITS 2 /* log calculation tolerance */
277 /* The fpexp10 fixed point math routine is based
278 * on oMathFP by Dan Carter (http://orbisstudios.com).
281 /** FIXED POINT EXP10
282 * Return 10^x as FP integer. Argument is FP integer.
284 static long long fp_exp10(long long x
, unsigned int fracbits
)
291 /* scale constants */
292 const long long fp_one
= (1 << fracbits
);
293 const long long fp_half
= (1 << (fracbits
- 1));
294 const long long fp_two
= (2 << fracbits
);
295 const long long fp_mask
= (fp_one
- 1);
296 const long long fp_ln2_inv
= (FP28_LN2_INV
>> (28 - fracbits
));
297 const long long fp_ln2
= (FP28_LN2
>> (28 - fracbits
));
298 const long long fp_ln10
= (FP28_LN10
>> (28 - fracbits
));
299 const long long fp_exp_zero
= (FP28_EXP_ZERO
>> (28 - fracbits
));
300 const long long fp_exp_one
= (FP28_EXP_ONE
>> (28 - fracbits
));
301 const long long fp_exp_two
= (FP28_EXP_TWO
>> (28 - fracbits
));
309 /* convert from base 10 to base e */
310 x
= FP_MUL64(x
, fp_ln10
);
312 /* calculate exp(x) */
313 k
= (FP_MUL64(abs(x
), fp_ln2_inv
) + fp_half
) & ~fp_mask
;
320 x
-= FP_MUL64(k
, fp_ln2
);
322 R
= fp_two
+ FP_MUL64(z
, fp_exp_zero
+ FP_MUL64(z
, fp_exp_one
323 + FP_MUL64(z
, fp_exp_two
)));
324 xp
= fp_one
+ FP_DIV64(FP_MUL64(fp_two
, x
), R
- x
);
328 k
= fp_one
>> (-k
>> fracbits
);
332 k
= fp_one
<< (k
>> fracbits
);
335 return FP_MUL64(k
, xp
);
339 /** FIXED POINT LOG10
340 * Return log10(x) as FP integer. Argument is FP integer.
342 static long long fp_log10(long long n
, unsigned int fracbits
)
344 /* Calculate log2 of argument */
346 long long log2
, frac
;
347 const long long fp_one
= (1 << fracbits
);
348 const long long fp_two
= (2 << fracbits
);
349 const long tolerance
= (1 << ((fracbits
/ 2) + 2));
351 if (n
<=0) return FP_NEGINF
;
366 /* fractional part */
368 while (frac
> tolerance
)
379 /* convert log2 to log10 */
380 return FP_MUL64(log2
, (FP28_LOG10OF2
>> (28 - fracbits
)));
384 /** CONVERT FACTOR TO DECIBELS */
385 long fp_decibels(unsigned long factor
, unsigned int fracbits
)
388 long long f
= (long long)factor
;
391 /* keep factor in signed long range */
392 if (f
>= (1LL << 31))
395 /* decibels = 20 * log10(factor) */
396 decibels
= FP_MUL64((20LL << fracbits
), fp_log10(f
, fracbits
));
398 /* keep result in signed long range */
399 if ((neg
= (decibels
< 0)))
400 decibels
= -decibels
;
401 if (decibels
>= (1LL << 31))
402 return neg
? FP_NEGINF
: FP_INF
;
404 return neg
? (long)-decibels
: (long)decibels
;
408 /** CONVERT DECIBELS TO FACTOR */
409 long fp_factor(long decibels
, unsigned int fracbits
)
413 long long db
= (long long)decibels
;
415 /* if decibels is 0, factor is 1 */
417 return (1L << fracbits
);
419 /* calculate for positive decibels only */
420 if ((neg
= (db
< 0)))
423 /* factor = 10 ^ (decibels / 20) */
424 factor
= fp_exp10(FP_DIV64(db
, (20LL << fracbits
)), fracbits
);
426 /* keep result in signed long range, return 0 if very small */
427 if (factor
>= (1LL << 31))
435 /* if negative argument, factor is 1 / result */
437 factor
= FP_DIV64((1LL << fracbits
), factor
);