summaryrefslogtreecommitdiff
path: root/apps
diff options
context:
space:
mode:
Diffstat (limited to 'apps')
-rw-r--r--apps/SOURCES1
-rw-r--r--apps/codecs/adx.c123
-rw-r--r--apps/codecs/lib/SOURCES2
-rw-r--r--apps/codecs/lib/fixedpoint.c1
-rw-r--r--apps/codecs/lib/fixedpoint.h49
-rw-r--r--apps/codecs/spc.c1
-rw-r--r--apps/dsp.c4
-rw-r--r--apps/dsp.h80
-rw-r--r--apps/eq.c120
-rw-r--r--apps/eq.h1
-rw-r--r--apps/fixedpoint.c452
-rw-r--r--apps/fixedpoint.h69
-rw-r--r--apps/fracmul.h93
-rw-r--r--apps/plugins/bounce.c4
-rw-r--r--apps/plugins/bubbles.c20
-rw-r--r--apps/plugins/clock/clock_draw_analog.c6
-rw-r--r--apps/plugins/cube.c12
-rw-r--r--apps/plugins/lib/fixedpoint.c239
-rw-r--r--apps/plugins/lib/fixedpoint.h45
-rw-r--r--apps/plugins/plasma.c2
-rw-r--r--apps/plugins/vu_meter.c4
-rw-r--r--apps/replaygain.c181
22 files changed, 754 insertions, 755 deletions
diff --git a/apps/SOURCES b/apps/SOURCES
index 7475826015..f3acef1739 100644
--- a/apps/SOURCES
+++ b/apps/SOURCES
@@ -125,6 +125,7 @@ recorder/recording.c
125#if INPUT_SRC_CAPS != 0 125#if INPUT_SRC_CAPS != 0
126audio_path.c 126audio_path.c
127#endif /* INPUT_SRC_CAPS != 0 */ 127#endif /* INPUT_SRC_CAPS != 0 */
128fixedpoint.c
128pcmbuf.c 129pcmbuf.c
129playback.c 130playback.c
130codecs.c 131codecs.c
diff --git a/apps/codecs/adx.c b/apps/codecs/adx.c
index cc36f6a908..0e50054753 100644
--- a/apps/codecs/adx.c
+++ b/apps/codecs/adx.c
@@ -21,6 +21,7 @@
21#include "codeclib.h" 21#include "codeclib.h"
22#include "inttypes.h" 22#include "inttypes.h"
23#include "math.h" 23#include "math.h"
24#include "lib/fixedpoint.h"
24 25
25CODEC_HEADER 26CODEC_HEADER
26 27
@@ -41,124 +42,6 @@ const long cutoff = 500;
41 42
42static int16_t samples[WAV_CHUNK_SIZE] IBSS_ATTR; 43static int16_t samples[WAV_CHUNK_SIZE] IBSS_ATTR;
43 44
44/* fixed point stuff from apps/plugins/lib/fixedpoint.c */
45
46/* Inverse gain of circular cordic rotation in s0.31 format. */
47static const long cordic_circular_gain = 0xb2458939; /* 0.607252929 */
48
49/* Table of values of atan(2^-i) in 0.32 format fractions of pi where pi = 0xffffffff / 2 */
50static const unsigned long atan_table[] = {
51 0x1fffffff, /* +0.785398163 (or pi/4) */
52 0x12e4051d, /* +0.463647609 */
53 0x09fb385b, /* +0.244978663 */
54 0x051111d4, /* +0.124354995 */
55 0x028b0d43, /* +0.062418810 */
56 0x0145d7e1, /* +0.031239833 */
57 0x00a2f61e, /* +0.015623729 */
58 0x00517c55, /* +0.007812341 */
59 0x0028be53, /* +0.003906230 */
60 0x00145f2e, /* +0.001953123 */
61 0x000a2f98, /* +0.000976562 */
62 0x000517cc, /* +0.000488281 */
63 0x00028be6, /* +0.000244141 */
64 0x000145f3, /* +0.000122070 */
65 0x0000a2f9, /* +0.000061035 */
66 0x0000517c, /* +0.000030518 */
67 0x000028be, /* +0.000015259 */
68 0x0000145f, /* +0.000007629 */
69 0x00000a2f, /* +0.000003815 */
70 0x00000517, /* +0.000001907 */
71 0x0000028b, /* +0.000000954 */
72 0x00000145, /* +0.000000477 */
73 0x000000a2, /* +0.000000238 */
74 0x00000051, /* +0.000000119 */
75 0x00000028, /* +0.000000060 */
76 0x00000014, /* +0.000000030 */
77 0x0000000a, /* +0.000000015 */
78 0x00000005, /* +0.000000007 */
79 0x00000002, /* +0.000000004 */
80 0x00000001, /* +0.000000002 */
81 0x00000000, /* +0.000000001 */
82 0x00000000, /* +0.000000000 */
83};
84
85/**
86 * Implements sin and cos using CORDIC rotation.
87 *
88 * @param phase has range from 0 to 0xffffffff, representing 0 and
89 * 2*pi respectively.
90 * @param cos return address for cos
91 * @return sin of phase, value is a signed value from LONG_MIN to LONG_MAX,
92 * representing -1 and 1 respectively.
93 */
94static long fsincos(unsigned long phase, long *cos)
95{
96 int32_t x, x1, y, y1;
97 unsigned long z, z1;
98 int i;
99
100 /* Setup initial vector */
101 x = cordic_circular_gain;
102 y = 0;
103 z = phase;
104
105 /* The phase has to be somewhere between 0..pi for this to work right */
106 if (z < 0xffffffff / 4) {
107 /* z in first quadrant, z += pi/2 to correct */
108 x = -x;
109 z += 0xffffffff / 4;
110 } else if (z < 3 * (0xffffffff / 4)) {
111 /* z in third quadrant, z -= pi/2 to correct */
112 z -= 0xffffffff / 4;
113 } else {
114 /* z in fourth quadrant, z -= 3pi/2 to correct */
115 x = -x;
116 z -= 3 * (0xffffffff / 4);
117 }
118
119 /* Each iteration adds roughly 1-bit of extra precision */
120 for (i = 0; i < 31; i++) {
121 x1 = x >> i;
122 y1 = y >> i;
123 z1 = atan_table[i];
124
125 /* Decided which direction to rotate vector. Pivot point is pi/2 */
126 if (z >= 0xffffffff / 4) {
127 x -= y1;
128 y += x1;
129 z -= z1;
130 } else {
131 x += y1;
132 y -= x1;
133 z += z1;
134 }
135 }
136
137 if (cos)
138 *cos = x;
139
140 return y;
141}
142
143/**
144 * Fixed point square root via Newton-Raphson.
145 * @param a square root argument.
146 * @param fracbits specifies number of fractional bits in argument.
147 * @return Square root of argument in same fixed point format as input.
148 */
149static long fsqrt(long a, unsigned int fracbits)
150{
151 long b = a/2 + (1 << fracbits); /* initial approximation */
152 unsigned n;
153 const unsigned iterations = 8; /* bumped up from 4 as it wasn't
154 nearly enough for 28 fractional bits */
155
156 for (n = 0; n < iterations; ++n)
157 b = (b + (long)(((long long)(a) << fracbits)/b))/2;
158
159 return b;
160}
161
162/* this is the codec entry point */ 45/* this is the codec entry point */
163enum codec_status codec_main(void) 46enum codec_status codec_main(void)
164{ 47{
@@ -238,7 +121,7 @@ next_track:
238 int64_t c; 121 int64_t c;
239 int64_t d; 122 int64_t d;
240 123
241 fsincos((unsigned long)phasemultiple,&z); 124 fp_sincos((unsigned long)phasemultiple,&z);
242 125
243 a = (M_SQRT2*big28)-(z*big28/LONG_MAX); 126 a = (M_SQRT2*big28)-(z*big28/LONG_MAX);
244 127
@@ -247,7 +130,7 @@ next_track:
247 * which is sufficient here, but this is the only reason why I don't 130 * which is sufficient here, but this is the only reason why I don't
248 * use 32 fractional bits everywhere. 131 * use 32 fractional bits everywhere.
249 */ 132 */
250 d = fsqrt((a+b)*(a-b)/big28,28); 133 d = fp_sqrt((a+b)*(a-b)/big28,28);
251 c = (a-d)*big28/b; 134 c = (a-d)*big28/b;
252 135
253 coef1 = (c*8192) >> 28; 136 coef1 = (c*8192) >> 28;
diff --git a/apps/codecs/lib/SOURCES b/apps/codecs/lib/SOURCES
index cbb8e60372..0141af21dc 100644
--- a/apps/codecs/lib/SOURCES
+++ b/apps/codecs/lib/SOURCES
@@ -1,6 +1,6 @@
1#if CONFIG_CODEC == SWCODEC /* software codec platforms */ 1#if CONFIG_CODEC == SWCODEC /* software codec platforms */
2codeclib.c 2codeclib.c
3 3fixedpoint.c
4 4
5mdct2.c 5mdct2.c
6#ifdef CPU_ARM 6#ifdef CPU_ARM
diff --git a/apps/codecs/lib/fixedpoint.c b/apps/codecs/lib/fixedpoint.c
new file mode 100644
index 0000000000..352e246673
--- /dev/null
+++ b/apps/codecs/lib/fixedpoint.c
@@ -0,0 +1 @@
#include "../../fixedpoint.c"
diff --git a/apps/codecs/lib/fixedpoint.h b/apps/codecs/lib/fixedpoint.h
new file mode 100644
index 0000000000..e6ed5208ca
--- /dev/null
+++ b/apps/codecs/lib/fixedpoint.h
@@ -0,0 +1,49 @@
1/***************************************************************************
2 * __________ __ ___.
3 * Open \______ \ ____ ____ | | _\_ |__ _______ ___
4 * Source | _// _ \_/ ___\| |/ /| __ \ / _ \ \/ /
5 * Jukebox | | ( <_> ) \___| < | \_\ ( <_> > < <
6 * Firmware |____|_ /\____/ \___ >__|_ \|___ /\____/__/\_ \
7 * \/ \/ \/ \/ \/
8 * $Id: fixedpoint.h -1 $
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 /** CODECS - FIXED POINT MATH ROUTINES - USAGE
25 *
26 * - x and y arguments are fixed point integers
27 * - fracbits is the number of fractional bits in the argument(s)
28 * - functions return long fixed point integers with the specified number
29 * of fractional bits unless otherwise specified
30 *
31 * Calculate sin and cos of an angle:
32 * fp_sincos(phase, *cos)
33 * where phase is a 32 bit unsigned integer with 0 representing 0
34 * and 0xFFFFFFFF representing 2*pi, and *cos is the address to
35 * a long signed integer. Value returned is a long signed integer
36 * from LONG_MIN to LONG_MAX, representing -1 to 1 respectively.
37 * That is, value is a fixed point integer with 31 fractional bits.
38 *
39 * Take square root of a fixed point number:
40 * fp_sqrt(x, fracbits)
41 *
42 */
43#ifndef _FIXEDPOINT_H_CODECS
44#define _FIXEDPOINT_H_CODECS
45
46long fp_sincos(unsigned long phase, long *cos);
47long fp_sqrt(long a, unsigned int fracbits);
48
49#endif
diff --git a/apps/codecs/spc.c b/apps/codecs/spc.c
index 6ceb704c7c..3a06a228ef 100644
--- a/apps/codecs/spc.c
+++ b/apps/codecs/spc.c
@@ -26,6 +26,7 @@
26/* DSP Based on Brad Martin's OpenSPC DSP emulator */ 26/* DSP Based on Brad Martin's OpenSPC DSP emulator */
27/* tag reading from sexyspc by John Brawn (John_Brawn@yahoo.com) and others */ 27/* tag reading from sexyspc by John Brawn (John_Brawn@yahoo.com) and others */
28#include "codeclib.h" 28#include "codeclib.h"
29#include "../fracmul.h"
29#include "libspc/spc_codec.h" 30#include "libspc/spc_codec.h"
30#include "libspc/spc_profiler.h" 31#include "libspc/spc_profiler.h"
31 32
diff --git a/apps/dsp.c b/apps/dsp.c
index a760865afb..30b4ed357b 100644
--- a/apps/dsp.c
+++ b/apps/dsp.c
@@ -33,6 +33,8 @@
33#include "misc.h" 33#include "misc.h"
34#include "tdspeed.h" 34#include "tdspeed.h"
35#include "buffer.h" 35#include "buffer.h"
36#include "fixedpoint.h"
37#include "fracmul.h"
36 38
37/* 16-bit samples are scaled based on these constants. The shift should be 39/* 16-bit samples are scaled based on these constants. The shift should be
38 * no more than 15. 40 * no more than 15.
@@ -841,7 +843,7 @@ void dsp_set_crossfeed_cross_params(long lf_gain, long hf_gain, long cutoff)
841 * crossfeed shelf filter and should be removed if crossfeed settings are 843 * crossfeed shelf filter and should be removed if crossfeed settings are
842 * ever made incompatible for any other good reason. 844 * ever made incompatible for any other good reason.
843 */ 845 */
844 cutoff = DIV64(cutoff, get_replaygain_int(hf_gain*5), 24); 846 cutoff = fp_div(cutoff, get_replaygain_int(hf_gain*5), 24);
845 filter_shelf_coefs(cutoff, hf_gain, false, c); 847 filter_shelf_coefs(cutoff, hf_gain, false, c);
846 /* Scale coefs by LF gain and shift them to s0.31 format. We have no gains 848 /* Scale coefs by LF gain and shift them to s0.31 format. We have no gains
847 * over 1 and can do this safely 849 * over 1 and can do this safely
diff --git a/apps/dsp.h b/apps/dsp.h
index 8c23c3053d..3d24b24245 100644
--- a/apps/dsp.h
+++ b/apps/dsp.h
@@ -64,86 +64,6 @@ enum {
64 DSP_CALLBACK_SET_STEREO_WIDTH 64 DSP_CALLBACK_SET_STEREO_WIDTH
65}; 65};
66 66
67/* A bunch of fixed point assembler helper macros */
68#if defined(CPU_COLDFIRE)
69/* These macros use the Coldfire EMAC extension and need the MACSR flags set
70 * to fractional mode with no rounding.
71 */
72
73/* Multiply two S.31 fractional integers and return the sign bit and the
74 * 31 most significant bits of the result.
75 */
76#define FRACMUL(x, y) \
77({ \
78 long t; \
79 asm ("mac.l %[a], %[b], %%acc0\n\t" \
80 "movclr.l %%acc0, %[t]\n\t" \
81 : [t] "=r" (t) : [a] "r" (x), [b] "r" (y)); \
82 t; \
83})
84
85/* Multiply two S.31 fractional integers, and return the 32 most significant
86 * bits after a shift left by the constant z. NOTE: Only works for shifts of
87 * 1 to 8 on Coldfire!
88 */
89#define FRACMUL_SHL(x, y, z) \
90({ \
91 long t, t2; \
92 asm ("mac.l %[a], %[b], %%acc0\n\t" \
93 "moveq.l %[d], %[t]\n\t" \
94 "move.l %%accext01, %[t2]\n\t" \
95 "and.l %[mask], %[t2]\n\t" \
96 "lsr.l %[t], %[t2]\n\t" \
97 "movclr.l %%acc0, %[t]\n\t" \
98 "asl.l %[c], %[t]\n\t" \
99 "or.l %[t2], %[t]\n\t" \
100 : [t] "=&d" (t), [t2] "=&d" (t2) \
101 : [a] "r" (x), [b] "r" (y), [mask] "d" (0xff), \
102 [c] "i" ((z)), [d] "i" (8 - (z))); \
103 t; \
104})
105
106#elif defined(CPU_ARM)
107
108/* Multiply two S.31 fractional integers and return the sign bit and the
109 * 31 most significant bits of the result.
110 */
111#define FRACMUL(x, y) \
112({ \
113 long t, t2; \
114 asm ("smull %[t], %[t2], %[a], %[b]\n\t" \
115 "mov %[t2], %[t2], asl #1\n\t" \
116 "orr %[t], %[t2], %[t], lsr #31\n\t" \
117 : [t] "=&r" (t), [t2] "=&r" (t2) \
118 : [a] "r" (x), [b] "r" (y)); \
119 t; \
120})
121
122/* Multiply two S.31 fractional integers, and return the 32 most significant
123 * bits after a shift left by the constant z.
124 */
125#define FRACMUL_SHL(x, y, z) \
126({ \
127 long t, t2; \
128 asm ("smull %[t], %[t2], %[a], %[b]\n\t" \
129 "mov %[t2], %[t2], asl %[c]\n\t" \
130 "orr %[t], %[t2], %[t], lsr %[d]\n\t" \
131 : [t] "=&r" (t), [t2] "=&r" (t2) \
132 : [a] "r" (x), [b] "r" (y), \
133 [c] "M" ((z) + 1), [d] "M" (31 - (z))); \
134 t; \
135})
136
137#else
138
139#define FRACMUL(x, y) (long) (((((long long) (x)) * ((long long) (y))) >> 31))
140#define FRACMUL_SHL(x, y, z) \
141((long)(((((long long) (x)) * ((long long) (y))) >> (31 - (z)))))
142
143#endif
144
145#define DIV64(x, y, z) (long)(((long long)(x) << (z))/(y))
146
147struct dsp_config; 67struct dsp_config;
148 68
149int dsp_process(struct dsp_config *dsp, char *dest, 69int dsp_process(struct dsp_config *dsp, char *dest,
diff --git a/apps/eq.c b/apps/eq.c
index 5977200c9c..6437fed906 100644
--- a/apps/eq.c
+++ b/apps/eq.c
@@ -21,105 +21,11 @@
21 21
22#include <inttypes.h> 22#include <inttypes.h>
23#include "config.h" 23#include "config.h"
24#include "dsp.h" 24#include "fixedpoint.h"
25#include "fracmul.h"
25#include "eq.h" 26#include "eq.h"
26#include "replaygain.h" 27#include "replaygain.h"
27 28
28/* Inverse gain of circular cordic rotation in s0.31 format. */
29static const long cordic_circular_gain = 0xb2458939; /* 0.607252929 */
30
31/* Table of values of atan(2^-i) in 0.32 format fractions of pi where pi = 0xffffffff / 2 */
32static const unsigned long atan_table[] = {
33 0x1fffffff, /* +0.785398163 (or pi/4) */
34 0x12e4051d, /* +0.463647609 */
35 0x09fb385b, /* +0.244978663 */
36 0x051111d4, /* +0.124354995 */
37 0x028b0d43, /* +0.062418810 */
38 0x0145d7e1, /* +0.031239833 */
39 0x00a2f61e, /* +0.015623729 */
40 0x00517c55, /* +0.007812341 */
41 0x0028be53, /* +0.003906230 */
42 0x00145f2e, /* +0.001953123 */
43 0x000a2f98, /* +0.000976562 */
44 0x000517cc, /* +0.000488281 */
45 0x00028be6, /* +0.000244141 */
46 0x000145f3, /* +0.000122070 */
47 0x0000a2f9, /* +0.000061035 */
48 0x0000517c, /* +0.000030518 */
49 0x000028be, /* +0.000015259 */
50 0x0000145f, /* +0.000007629 */
51 0x00000a2f, /* +0.000003815 */
52 0x00000517, /* +0.000001907 */
53 0x0000028b, /* +0.000000954 */
54 0x00000145, /* +0.000000477 */
55 0x000000a2, /* +0.000000238 */
56 0x00000051, /* +0.000000119 */
57 0x00000028, /* +0.000000060 */
58 0x00000014, /* +0.000000030 */
59 0x0000000a, /* +0.000000015 */
60 0x00000005, /* +0.000000007 */
61 0x00000002, /* +0.000000004 */
62 0x00000001, /* +0.000000002 */
63 0x00000000, /* +0.000000001 */
64 0x00000000, /* +0.000000000 */
65};
66
67/**
68 * Implements sin and cos using CORDIC rotation.
69 *
70 * @param phase has range from 0 to 0xffffffff, representing 0 and
71 * 2*pi respectively.
72 * @param cos return address for cos
73 * @return sin of phase, value is a signed value from LONG_MIN to LONG_MAX,
74 * representing -1 and 1 respectively.
75 */
76static long fsincos(unsigned long phase, long *cos) {
77 int32_t x, x1, y, y1;
78 unsigned long z, z1;
79 int i;
80
81 /* Setup initial vector */
82 x = cordic_circular_gain;
83 y = 0;
84 z = phase;
85
86 /* The phase has to be somewhere between 0..pi for this to work right */
87 if (z < 0xffffffff / 4) {
88 /* z in first quadrant, z += pi/2 to correct */
89 x = -x;
90 z += 0xffffffff / 4;
91 } else if (z < 3 * (0xffffffff / 4)) {
92 /* z in third quadrant, z -= pi/2 to correct */
93 z -= 0xffffffff / 4;
94 } else {
95 /* z in fourth quadrant, z -= 3pi/2 to correct */
96 x = -x;
97 z -= 3 * (0xffffffff / 4);
98 }
99
100 /* Each iteration adds roughly 1-bit of extra precision */
101 for (i = 0; i < 31; i++) {
102 x1 = x >> i;
103 y1 = y >> i;
104 z1 = atan_table[i];
105
106 /* Decided which direction to rotate vector. Pivot point is pi/2 */
107 if (z >= 0xffffffff / 4) {
108 x -= y1;
109 y += x1;
110 z -= z1;
111 } else {
112 x += y1;
113 y -= x1;
114 z += z1;
115 }
116 }
117
118 *cos = x;
119
120 return y;
121}
122
123/** 29/**
124 * Calculate first order shelving filter. Filter is not directly usable by the 30 * Calculate first order shelving filter. Filter is not directly usable by the
125 * eq_filter() function. 31 * eq_filter() function.
@@ -135,16 +41,16 @@ void filter_shelf_coefs(unsigned long cutoff, long A, bool low, int32_t *c)
135 int32_t b0, b1, a0, a1; /* s3.28 */ 41 int32_t b0, b1, a0, a1; /* s3.28 */
136 const long g = get_replaygain_int(A*5) << 4; /* 10^(db/40), s3.28 */ 42 const long g = get_replaygain_int(A*5) << 4; /* 10^(db/40), s3.28 */
137 43
138 sin = fsincos(cutoff/2, &cos); 44 sin = fp_sincos(cutoff/2, &cos);
139 if (low) { 45 if (low) {
140 const int32_t sin_div_g = DIV64(sin, g, 25); 46 const int32_t sin_div_g = fp_div(sin, g, 25);
141 cos >>= 3; 47 cos >>= 3;
142 b0 = FRACMUL(sin, g) + cos; /* 0.25 .. 4.10 */ 48 b0 = FRACMUL(sin, g) + cos; /* 0.25 .. 4.10 */
143 b1 = FRACMUL(sin, g) - cos; /* -1 .. 3.98 */ 49 b1 = FRACMUL(sin, g) - cos; /* -1 .. 3.98 */
144 a0 = sin_div_g + cos; /* 0.25 .. 4.10 */ 50 a0 = sin_div_g + cos; /* 0.25 .. 4.10 */
145 a1 = sin_div_g - cos; /* -1 .. 3.98 */ 51 a1 = sin_div_g - cos; /* -1 .. 3.98 */
146 } else { 52 } else {
147 const int32_t cos_div_g = DIV64(cos, g, 25); 53 const int32_t cos_div_g = fp_div(cos, g, 25);
148 sin >>= 3; 54 sin >>= 3;
149 b0 = sin + FRACMUL(cos, g); /* 0.25 .. 4.10 */ 55 b0 = sin + FRACMUL(cos, g); /* 0.25 .. 4.10 */
150 b1 = sin - FRACMUL(cos, g); /* -3.98 .. 1 */ 56 b1 = sin - FRACMUL(cos, g); /* -3.98 .. 1 */
@@ -152,7 +58,7 @@ void filter_shelf_coefs(unsigned long cutoff, long A, bool low, int32_t *c)
152 a1 = sin - cos_div_g; /* -3.98 .. 1 */ 58 a1 = sin - cos_div_g; /* -3.98 .. 1 */
153 } 59 }
154 60
155 const int32_t rcp_a0 = DIV64(1, a0, 57); /* 0.24 .. 3.98, s2.29 */ 61 const int32_t rcp_a0 = fp_div(1, a0, 57); /* 0.24 .. 3.98, s2.29 */
156 *c++ = FRACMUL_SHL(b0, rcp_a0, 1); /* 0.063 .. 15.85 */ 62 *c++ = FRACMUL_SHL(b0, rcp_a0, 1); /* 0.063 .. 15.85 */
157 *c++ = FRACMUL_SHL(b1, rcp_a0, 1); /* -15.85 .. 15.85 */ 63 *c++ = FRACMUL_SHL(b1, rcp_a0, 1); /* -15.85 .. 15.85 */
158 *c++ = -FRACMUL_SHL(a1, rcp_a0, 1); /* -1 .. 1 */ 64 *c++ = -FRACMUL_SHL(a1, rcp_a0, 1); /* -1 .. 1 */
@@ -220,10 +126,10 @@ void eq_pk_coefs(unsigned long cutoff, unsigned long Q, long db, int32_t *c)
220 long cs; 126 long cs;
221 const long one = 1 << 28; /* s3.28 */ 127 const long one = 1 << 28; /* s3.28 */
222 const long A = get_replaygain_int(db*5) << 5; /* 10^(db/40), s2.29 */ 128 const long A = get_replaygain_int(db*5) << 5; /* 10^(db/40), s2.29 */
223 const long alpha = fsincos(cutoff, &cs)/(2*Q)*10 >> 1; /* s1.30 */ 129 const long alpha = fp_sincos(cutoff, &cs)/(2*Q)*10 >> 1; /* s1.30 */
224 int32_t a0, a1, a2; /* these are all s3.28 format */ 130 int32_t a0, a1, a2; /* these are all s3.28 format */
225 int32_t b0, b1, b2; 131 int32_t b0, b1, b2;
226 const long alphadivA = DIV64(alpha, A, 27); 132 const long alphadivA = fp_div(alpha, A, 27);
227 133
228 /* possible numerical ranges are in comments by each coef */ 134 /* possible numerical ranges are in comments by each coef */
229 b0 = one + FRACMUL(alpha, A); /* [1 .. 5] */ 135 b0 = one + FRACMUL(alpha, A); /* [1 .. 5] */
@@ -233,7 +139,7 @@ void eq_pk_coefs(unsigned long cutoff, unsigned long Q, long db, int32_t *c)
233 a2 = one - alphadivA; /* [-3 .. 1] */ 139 a2 = one - alphadivA; /* [-3 .. 1] */
234 140
235 /* range of this is roughly [0.2 .. 1], but we'll never hit 1 completely */ 141 /* range of this is roughly [0.2 .. 1], but we'll never hit 1 completely */
236 const long rcp_a0 = DIV64(1, a0, 59); /* s0.31 */ 142 const long rcp_a0 = fp_div(1, a0, 59); /* s0.31 */
237 *c++ = FRACMUL(b0, rcp_a0); /* [0.25 .. 4] */ 143 *c++ = FRACMUL(b0, rcp_a0); /* [0.25 .. 4] */
238 *c++ = FRACMUL(b1, rcp_a0); /* [-2 .. 2] */ 144 *c++ = FRACMUL(b1, rcp_a0); /* [-2 .. 2] */
239 *c++ = FRACMUL(b2, rcp_a0); /* [-2.4 .. 1] */ 145 *c++ = FRACMUL(b2, rcp_a0); /* [-2.4 .. 1] */
@@ -251,7 +157,7 @@ void eq_ls_coefs(unsigned long cutoff, unsigned long Q, long db, int32_t *c)
251 const long one = 1 << 25; /* s6.25 */ 157 const long one = 1 << 25; /* s6.25 */
252 const long sqrtA = get_replaygain_int(db*5/2) << 2; /* 10^(db/80), s5.26 */ 158 const long sqrtA = get_replaygain_int(db*5/2) << 2; /* 10^(db/80), s5.26 */
253 const long A = FRACMUL_SHL(sqrtA, sqrtA, 8); /* s2.29 */ 159 const long A = FRACMUL_SHL(sqrtA, sqrtA, 8); /* s2.29 */
254 const long alpha = fsincos(cutoff, &cs)/(2*Q)*10 >> 1; /* s1.30 */ 160 const long alpha = fp_sincos(cutoff, &cs)/(2*Q)*10 >> 1; /* s1.30 */
255 const long ap1 = (A >> 4) + one; 161 const long ap1 = (A >> 4) + one;
256 const long am1 = (A >> 4) - one; 162 const long am1 = (A >> 4) - one;
257 const long twosqrtalpha = 2*FRACMUL(sqrtA, alpha); 163 const long twosqrtalpha = 2*FRACMUL(sqrtA, alpha);
@@ -272,7 +178,7 @@ void eq_ls_coefs(unsigned long cutoff, unsigned long Q, long db, int32_t *c)
272 a2 = ap1 + FRACMUL(am1, cs) - twosqrtalpha; 178 a2 = ap1 + FRACMUL(am1, cs) - twosqrtalpha;
273 179
274 /* [0.1 .. 1.99] */ 180 /* [0.1 .. 1.99] */
275 const long rcp_a0 = DIV64(1, a0, 55); /* s1.30 */ 181 const long rcp_a0 = fp_div(1, a0, 55); /* s1.30 */
276 *c++ = FRACMUL_SHL(b0, rcp_a0, 2); /* [0.06 .. 15.9] */ 182 *c++ = FRACMUL_SHL(b0, rcp_a0, 2); /* [0.06 .. 15.9] */
277 *c++ = FRACMUL_SHL(b1, rcp_a0, 2); /* [-2 .. 31.7] */ 183 *c++ = FRACMUL_SHL(b1, rcp_a0, 2); /* [-2 .. 31.7] */
278 *c++ = FRACMUL_SHL(b2, rcp_a0, 2); /* [0 .. 15.9] */ 184 *c++ = FRACMUL_SHL(b2, rcp_a0, 2); /* [0 .. 15.9] */
@@ -290,7 +196,7 @@ void eq_hs_coefs(unsigned long cutoff, unsigned long Q, long db, int32_t *c)
290 const long one = 1 << 25; /* s6.25 */ 196 const long one = 1 << 25; /* s6.25 */
291 const long sqrtA = get_replaygain_int(db*5/2) << 2; /* 10^(db/80), s5.26 */ 197 const long sqrtA = get_replaygain_int(db*5/2) << 2; /* 10^(db/80), s5.26 */
292 const long A = FRACMUL_SHL(sqrtA, sqrtA, 8); /* s2.29 */ 198 const long A = FRACMUL_SHL(sqrtA, sqrtA, 8); /* s2.29 */
293 const long alpha = fsincos(cutoff, &cs)/(2*Q)*10 >> 1; /* s1.30 */ 199 const long alpha = fp_sincos(cutoff, &cs)/(2*Q)*10 >> 1; /* s1.30 */
294 const long ap1 = (A >> 4) + one; 200 const long ap1 = (A >> 4) + one;
295 const long am1 = (A >> 4) - one; 201 const long am1 = (A >> 4) - one;
296 const long twosqrtalpha = 2*FRACMUL(sqrtA, alpha); 202 const long twosqrtalpha = 2*FRACMUL(sqrtA, alpha);
@@ -311,7 +217,7 @@ void eq_hs_coefs(unsigned long cutoff, unsigned long Q, long db, int32_t *c)
311 a2 = ap1 - FRACMUL(am1, cs) - twosqrtalpha; 217 a2 = ap1 - FRACMUL(am1, cs) - twosqrtalpha;
312 218
313 /* [0.1 .. 1.99] */ 219 /* [0.1 .. 1.99] */
314 const long rcp_a0 = DIV64(1, a0, 55); /* s1.30 */ 220 const long rcp_a0 = fp_div(1, a0, 55); /* s1.30 */
315 *c++ = FRACMUL_SHL(b0, rcp_a0, 2); /* [0 .. 16] */ 221 *c++ = FRACMUL_SHL(b0, rcp_a0, 2); /* [0 .. 16] */
316 *c++ = FRACMUL_SHL(b1, rcp_a0, 2); /* [-31.7 .. 2] */ 222 *c++ = FRACMUL_SHL(b1, rcp_a0, 2); /* [-31.7 .. 2] */
317 *c++ = FRACMUL_SHL(b2, rcp_a0, 2); /* [0 .. 16] */ 223 *c++ = FRACMUL_SHL(b2, rcp_a0, 2); /* [0 .. 16] */
diff --git a/apps/eq.h b/apps/eq.h
index 1c3efe50e9..a44e9153ac 100644
--- a/apps/eq.h
+++ b/apps/eq.h
@@ -23,6 +23,7 @@
23#define _EQ_H 23#define _EQ_H
24 24
25#include <inttypes.h> 25#include <inttypes.h>
26#include <stdbool.h>
26 27
27/* These depend on the fixed point formats used by the different filter types 28/* These depend on the fixed point formats used by the different filter types
28 and need to be changed when they change. 29 and need to be changed when they change.
diff --git a/apps/fixedpoint.c b/apps/fixedpoint.c
new file mode 100644
index 0000000000..917f624258
--- /dev/null
+++ b/apps/fixedpoint.c
@@ -0,0 +1,452 @@
1/***************************************************************************
2 * __________ __ ___.
3 * Open \______ \ ____ ____ | | _\_ |__ _______ ___
4 * Source | _// _ \_/ ___\| |/ /| __ \ / _ \ \/ /
5 * Jukebox | | ( <_> ) \___| < | \_\ ( <_> > < <
6 * Firmware |____|_ /\____/ \___ >__|_ \|___ /\____/__/\_ \
7 * \/ \/ \/ \/ \/
8 * $Id: fixedpoint.c -1 $
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. In
156 * general, the closer x is to 1, the quicker the calculation.
157 */
158long fp_sqrt(long x, unsigned int fracbits)
159{
160 long b = x/2 + BIT_N(fracbits); /* initial approximation */
161 long c;
162 unsigned n;
163 const unsigned iterations = 8;
164
165 for (n = 0; n < iterations; ++n)
166 {
167 c = fp_div(x, b, fracbits);
168 if (c == b) break;
169 b = (b + c)/2;
170 }
171
172 return b;
173}
174#endif /* PLUGIN or CODEC */
175
176
177#if defined(PLUGIN)
178/**
179 * Fixed point sinus using a lookup table
180 * don't forget to divide the result by 16384 to get the actual sinus value
181 * @param val sinus argument in degree
182 * @return sin(val)*16384
183 */
184long fp14_sin(int val)
185{
186 val = (val+360)%360;
187 if (val < 181)
188 {
189 if (val < 91)/* phase 0-90 degree */
190 return (long)sin_table[val];
191 else/* phase 91-180 degree */
192 return (long)sin_table[180-val];
193 }
194 else
195 {
196 if (val < 271)/* phase 181-270 degree */
197 return -(long)sin_table[val-180];
198 else/* phase 270-359 degree */
199 return -(long)sin_table[360-val];
200 }
201 return 0;
202}
203
204/**
205 * Fixed point cosinus using a lookup table
206 * don't forget to divide the result by 16384 to get the actual cosinus value
207 * @param val sinus argument in degree
208 * @return cos(val)*16384
209 */
210long fp14_cos(int val)
211{
212 val = (val+360)%360;
213 if (val < 181)
214 {
215 if (val < 91)/* phase 0-90 degree */
216 return (long)sin_table[90-val];
217 else/* phase 91-180 degree */
218 return -(long)sin_table[val-90];
219 }
220 else
221 {
222 if (val < 271)/* phase 181-270 degree */
223 return -(long)sin_table[270-val];
224 else/* phase 270-359 degree */
225 return (long)sin_table[val-270];
226 }
227 return 0;
228}
229
230/**
231 * Fixed-point natural log
232 * taken from http://www.quinapalus.com/efunc.html
233 * "The code assumes integers are at least 32 bits long. The (positive)
234 * argument and the result of the function are both expressed as fixed-point
235 * values with 16 fractional bits, although intermediates are kept with 28
236 * bits of precision to avoid loss of accuracy during shifts."
237 */
238
239long fp16_log(int x) {
240 long t,y;
241
242 y=0xa65af;
243 if(x<0x00008000) x<<=16, y-=0xb1721;
244 if(x<0x00800000) x<<= 8, y-=0x58b91;
245 if(x<0x08000000) x<<= 4, y-=0x2c5c8;
246 if(x<0x20000000) x<<= 2, y-=0x162e4;
247 if(x<0x40000000) x<<= 1, y-=0x0b172;
248 t=x+(x>>1); if((t&0x80000000)==0) x=t,y-=0x067cd;
249 t=x+(x>>2); if((t&0x80000000)==0) x=t,y-=0x03920;
250 t=x+(x>>3); if((t&0x80000000)==0) x=t,y-=0x01e27;
251 t=x+(x>>4); if((t&0x80000000)==0) x=t,y-=0x00f85;
252 t=x+(x>>5); if((t&0x80000000)==0) x=t,y-=0x007e1;
253 t=x+(x>>6); if((t&0x80000000)==0) x=t,y-=0x003f8;
254 t=x+(x>>7); if((t&0x80000000)==0) x=t,y-=0x001fe;
255 x=0x80000000-x;
256 y-=x>>15;
257 return y;
258}
259#endif /* PLUGIN */
260
261
262#if (!defined(PLUGIN) && !defined(CODEC))
263/** MODIFIED FROM replaygain.c */
264/* These math routines have 64-bit internal precision to avoid overflows.
265 * Arguments and return values are 32-bit (long) precision.
266 */
267
268#define FP_MUL64(x, y) (((x) * (y)) >> (fracbits))
269#define FP_DIV64(x, y) (((x) << (fracbits)) / (y))
270
271static long long fp_exp10(long long x, unsigned int fracbits);
272/* static long long fp_log10(long long n, unsigned int fracbits); */
273
274/* constants in fixed point format, 28 fractional bits */
275#define FP28_LN2 (186065279LL) /* ln(2) */
276#define FP28_LN2_INV (387270501LL) /* 1/ln(2) */
277#define FP28_EXP_ZERO (44739243LL) /* 1/6 */
278#define FP28_EXP_ONE (-745654LL) /* -1/360 */
279#define FP28_EXP_TWO (12428LL) /* 1/21600 */
280#define FP28_LN10 (618095479LL) /* ln(10) */
281#define FP28_LOG10OF2 (80807124LL) /* log10(2) */
282
283#define TOL_BITS 2 /* log calculation tolerance */
284
285
286/* The fpexp10 fixed point math routine is based
287 * on oMathFP by Dan Carter (http://orbisstudios.com).
288 */
289
290/** FIXED POINT EXP10
291 * Return 10^x as FP integer. Argument is FP integer.
292 */
293static long long fp_exp10(long long x, unsigned int fracbits)
294{
295 long long k;
296 long long z;
297 long long R;
298 long long xp;
299
300 /* scale constants */
301 const long long fp_one = (1 << fracbits);
302 const long long fp_half = (1 << (fracbits - 1));
303 const long long fp_two = (2 << fracbits);
304 const long long fp_mask = (fp_one - 1);
305 const long long fp_ln2_inv = (FP28_LN2_INV >> (28 - fracbits));
306 const long long fp_ln2 = (FP28_LN2 >> (28 - fracbits));
307 const long long fp_ln10 = (FP28_LN10 >> (28 - fracbits));
308 const long long fp_exp_zero = (FP28_EXP_ZERO >> (28 - fracbits));
309 const long long fp_exp_one = (FP28_EXP_ONE >> (28 - fracbits));
310 const long long fp_exp_two = (FP28_EXP_TWO >> (28 - fracbits));
311
312 /* exp(0) = 1 */
313 if (x == 0)
314 {
315 return fp_one;
316 }
317
318 /* convert from base 10 to base e */
319 x = FP_MUL64(x, fp_ln10);
320
321 /* calculate exp(x) */
322 k = (FP_MUL64(abs(x), fp_ln2_inv) + fp_half) & ~fp_mask;
323
324 if (x < 0)
325 {
326 k = -k;
327 }
328
329 x -= FP_MUL64(k, fp_ln2);
330 z = FP_MUL64(x, x);
331 R = fp_two + FP_MUL64(z, fp_exp_zero + FP_MUL64(z, fp_exp_one
332 + FP_MUL64(z, fp_exp_two)));
333 xp = fp_one + FP_DIV64(FP_MUL64(fp_two, x), R - x);
334
335 if (k < 0)
336 {
337 k = fp_one >> (-k >> fracbits);
338 }
339 else
340 {
341 k = fp_one << (k >> fracbits);
342 }
343
344 return FP_MUL64(k, xp);
345}
346
347
348#if 0 /* useful code, but not currently used */
349/** FIXED POINT LOG10
350 * Return log10(x) as FP integer. Argument is FP integer.
351 */
352static long long fp_log10(long long n, unsigned int fracbits)
353{
354 /* Calculate log2 of argument */
355
356 long long log2, frac;
357 const long long fp_one = (1 << fracbits);
358 const long long fp_two = (2 << fracbits);
359 const long tolerance = (1 << ((fracbits / 2) + 2));
360
361 if (n <=0) return FP_NEGINF;
362 log2 = 0;
363
364 /* integer part */
365 while (n < fp_one)
366 {
367 log2 -= fp_one;
368 n <<= 1;
369 }
370 while (n >= fp_two)
371 {
372 log2 += fp_one;
373 n >>= 1;
374 }
375
376 /* fractional part */
377 frac = fp_one;
378 while (frac > tolerance)
379 {
380 frac >>= 1;
381 n = FP_MUL64(n, n);
382 if (n >= fp_two)
383 {
384 n >>= 1;
385 log2 += frac;
386 }
387 }
388
389 /* convert log2 to log10 */
390 return FP_MUL64(log2, (FP28_LOG10OF2 >> (28 - fracbits)));
391}
392
393
394/** CONVERT FACTOR TO DECIBELS */
395long fp_decibels(unsigned long factor, unsigned int fracbits)
396{
397 long long decibels;
398 long long f = (long long)factor;
399 bool neg;
400
401 /* keep factor in signed long range */
402 if (f >= (1LL << 31))
403 f = (1LL << 31) - 1;
404
405 /* decibels = 20 * log10(factor) */
406 decibels = FP_MUL64((20LL << fracbits), fp_log10(f, fracbits));
407
408 /* keep result in signed long range */
409 if ((neg = (decibels < 0)))
410 decibels = -decibels;
411 if (decibels >= (1LL << 31))
412 return neg ? FP_NEGINF : FP_INF;
413
414 return neg ? (long)-decibels : (long)decibels;
415}
416#endif /* unused code */
417
418
419/** CONVERT DECIBELS TO FACTOR */
420long fp_factor(long decibels, unsigned int fracbits)
421{
422 bool neg;
423 long long factor;
424 long long db = (long long)decibels;
425
426 /* if decibels is 0, factor is 1 */
427 if (db == 0)
428 return (1L << fracbits);
429
430 /* calculate for positive decibels only */
431 if ((neg = (db < 0)))
432 db = -db;
433
434 /* factor = 10 ^ (decibels / 20) */
435 factor = fp_exp10(FP_DIV64(db, (20LL << fracbits)), fracbits);
436
437 /* keep result in signed long range, return 0 if very small */
438 if (factor >= (1LL << 31))
439 {
440 if (neg)
441 return 0;
442 else
443 return FP_INF;
444 }
445
446 /* if negative argument, factor is 1 / result */
447 if (neg)
448 factor = FP_DIV64((1LL << fracbits), factor);
449
450 return (long)factor;
451}
452#endif /* !PLUGIN and !CODEC */
diff --git a/apps/fixedpoint.h b/apps/fixedpoint.h
new file mode 100644
index 0000000000..49292f3417
--- /dev/null
+++ b/apps/fixedpoint.h
@@ -0,0 +1,69 @@
1/***************************************************************************
2 * __________ __ ___.
3 * Open \______ \ ____ ____ | | _\_ |__ _______ ___
4 * Source | _// _ \_/ ___\| |/ /| __ \ / _ \ \/ /
5 * Jukebox | | ( <_> ) \___| < | \_\ ( <_> > < <
6 * Firmware |____|_ /\____/ \___ >__|_ \|___ /\____/__/\_ \
7 * \/ \/ \/ \/ \/
8 * $Id: fixedpoint.h -1 $
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/** APPS - FIXED POINT MATH ROUTINES - USAGE
25 *
26 * - x and y arguments are fixed point integers
27 * - fracbits is the number of fractional bits in the argument(s)
28 * - functions return long fixed point integers with the specified number
29 * of fractional bits unless otherwise specified
30 *
31 * Multiply two fixed point numbers:
32 * fp_mul(x, y, fracbits)
33 *
34 * Divide two fixed point numbers:
35 * fp_div(x, y, fracbits)
36 *
37 * Calculate decibel equivalent of a gain factor:
38 * fp_decibels(factor, fracbits)
39 * where fracbits is in the range 12 to 22 (higher is better),
40 * and factor is a positive fixed point integer.
41 *
42 * Calculate factor equivalent of a decibel value:
43 * fp_factor(decibels, fracbits)
44 * where fracbits is in the range 12 to 22 (lower is better),
45 * and decibels is a fixed point integer.
46 */
47
48#ifndef _FIXEDPOINT_H_APPS
49#define _FIXEDPOINT_H_APPS
50
51#define fp_mul(x, y, z) (long)((((long long)(x)) * ((long long)(y))) >> (z))
52#define fp_div(x, y, z) (long)((((long long)(x)) << (z)) / ((long long)(y)))
53
54
55/** TAKEN FROM ORIGINAL fixedpoint.h */
56long fp_sincos(unsigned long phase, long *cos);
57
58
59/** MODIFIED FROM replaygain.c */
60#define FP_INF (0x7fffffff)
61#define FP_NEGINF -(0x7fffffff)
62
63/* fracbits in range 12 - 22 work well. Higher is better for
64 * calculating dB, lower is better for calculating ratio.
65 */
66/* long fp_decibels(unsigned long factor, unsigned int fracbits); */
67long fp_factor(long decibels, unsigned int fracbits);
68
69#endif
diff --git a/apps/fracmul.h b/apps/fracmul.h
new file mode 100644
index 0000000000..5cc83af624
--- /dev/null
+++ b/apps/fracmul.h
@@ -0,0 +1,93 @@
1#ifndef _FRACMUL_H
2#define _FRACMUL_H
3
4/** FRACTIONAL MULTIPLICATION - TAKEN FROM apps/dsp.h
5 * Multiply two fixed point numbers with 31 fractional bits:
6 * FRACMUL(x, y)
7 *
8 * Multiply two fixed point numbers with 31 fractional bits,
9 * then shift left by z bits:
10 * FRACMUL_SHL(x, y, z)
11 * NOTE: z must be in the range 1-8 on Coldfire targets.
12 */
13
14
15/* A bunch of fixed point assembler helper macros */
16#if defined(CPU_COLDFIRE)
17/* These macros use the Coldfire EMAC extension and need the MACSR flags set
18 * to fractional mode with no rounding.
19 */
20
21/* Multiply two S.31 fractional integers and return the sign bit and the
22 * 31 most significant bits of the result.
23 */
24#define FRACMUL(x, y) \
25({ \
26 long t; \
27 asm ("mac.l %[a], %[b], %%acc0\n\t" \
28 "movclr.l %%acc0, %[t]\n\t" \
29 : [t] "=r" (t) : [a] "r" (x), [b] "r" (y)); \
30 t; \
31})
32
33/* Multiply two S.31 fractional integers, and return the 32 most significant
34 * bits after a shift left by the constant z. NOTE: Only works for shifts of
35 * 1 to 8 on Coldfire!
36 */
37#define FRACMUL_SHL(x, y, z) \
38({ \
39 long t, t2; \
40 asm ("mac.l %[a], %[b], %%acc0\n\t" \
41 "moveq.l %[d], %[t]\n\t" \
42 "move.l %%accext01, %[t2]\n\t" \
43 "and.l %[mask], %[t2]\n\t" \
44 "lsr.l %[t], %[t2]\n\t" \
45 "movclr.l %%acc0, %[t]\n\t" \
46 "asl.l %[c], %[t]\n\t" \
47 "or.l %[t2], %[t]\n\t" \
48 : [t] "=&d" (t), [t2] "=&d" (t2) \
49 : [a] "r" (x), [b] "r" (y), [mask] "d" (0xff), \
50 [c] "i" ((z)), [d] "i" (8 - (z))); \
51 t; \
52})
53
54#elif defined(CPU_ARM)
55
56/* Multiply two S.31 fractional integers and return the sign bit and the
57 * 31 most significant bits of the result.
58 */
59#define FRACMUL(x, y) \
60({ \
61 long t, t2; \
62 asm ("smull %[t], %[t2], %[a], %[b]\n\t" \
63 "mov %[t2], %[t2], asl #1\n\t" \
64 "orr %[t], %[t2], %[t], lsr #31\n\t" \
65 : [t] "=&r" (t), [t2] "=&r" (t2) \
66 : [a] "r" (x), [b] "r" (y)); \
67 t; \
68})
69
70/* Multiply two S.31 fractional integers, and return the 32 most significant
71 * bits after a shift left by the constant z.
72 */
73#define FRACMUL_SHL(x, y, z) \
74({ \
75 long t, t2; \
76 asm ("smull %[t], %[t2], %[a], %[b]\n\t" \
77 "mov %[t2], %[t2], asl %[c]\n\t" \
78 "orr %[t], %[t2], %[t], lsr %[d]\n\t" \
79 : [t] "=&r" (t), [t2] "=&r" (t2) \
80 : [a] "r" (x), [b] "r" (y), \
81 [c] "M" ((z) + 1), [d] "M" (31 - (z))); \
82 t; \
83})
84
85#else
86
87#define FRACMUL(x, y) (long) (((((long long) (x)) * ((long long) (y))) >> 31))
88#define FRACMUL_SHL(x, y, z) \
89((long)(((((long long) (x)) * ((long long) (y))) >> (31 - (z)))))
90
91#endif
92
93#endif
diff --git a/apps/plugins/bounce.c b/apps/plugins/bounce.c
index ee4c3e443c..14bc7dea98 100644
--- a/apps/plugins/bounce.c
+++ b/apps/plugins/bounce.c
@@ -344,7 +344,7 @@ static void init_tables(void)
344 phase = pfrac = 0; 344 phase = pfrac = 0;
345 345
346 for (i = 0; i < TABLE_SIZE; i++) { 346 for (i = 0; i < TABLE_SIZE; i++) {
347 sin = fsincos(phase, NULL); 347 sin = fp_sincos(phase, NULL);
348 xtable[i] = RADIUS_X + sin / DIV_X; 348 xtable[i] = RADIUS_X + sin / DIV_X;
349 ytable[i] = RADIUS_Y + sin / DIV_Y; 349 ytable[i] = RADIUS_Y + sin / DIV_Y;
350 350
@@ -411,7 +411,7 @@ static void init_clock(void)
411 phase = pfrac = 0; 411 phase = pfrac = 0;
412 412
413 for (i = 0; i < 60; i++) { 413 for (i = 0; i < 60; i++) {
414 sin = fsincos(phase, &cos); 414 sin = fp_sincos(phase, &cos);
415 xminute[i] = LCD_WIDTH/2 + sin / DIV_MX; 415 xminute[i] = LCD_WIDTH/2 + sin / DIV_MX;
416 yminute[i] = LCD_HEIGHT/2 - cos / DIV_MY; 416 yminute[i] = LCD_HEIGHT/2 - cos / DIV_MY;
417 xhour[i] = LCD_WIDTH/2 + sin / DIV_HX; 417 xhour[i] = LCD_WIDTH/2 + sin / DIV_HX;
diff --git a/apps/plugins/bubbles.c b/apps/plugins/bubbles.c
index 4146b45b36..44d172c4ee 100644
--- a/apps/plugins/bubbles.c
+++ b/apps/plugins/bubbles.c
@@ -1469,17 +1469,17 @@ static void bubbles_drawboard(struct game_context* bb) {
1469 ROW_HEIGHT*(BB_HEIGHT-2)+BUBBLE_HEIGHT); 1469 ROW_HEIGHT*(BB_HEIGHT-2)+BUBBLE_HEIGHT);
1470 1470
1471 /* draw arrow */ 1471 /* draw arrow */
1472 tipx = SHOTX+BUBBLE_WIDTH/2+(((sin_int(bb->angle)>>4)*BUBBLE_WIDTH*3/2)>>10); 1472 tipx = SHOTX+BUBBLE_WIDTH/2+(((fp14_sin(bb->angle)>>4)*BUBBLE_WIDTH*3/2)>>10);
1473 tipy = SHOTY+BUBBLE_HEIGHT/2-(((cos_int(bb->angle)>>4)*BUBBLE_HEIGHT*3/2)>>10); 1473 tipy = SHOTY+BUBBLE_HEIGHT/2-(((fp14_cos(bb->angle)>>4)*BUBBLE_HEIGHT*3/2)>>10);
1474 1474
1475 rb->lcd_drawline(SHOTX+BUBBLE_WIDTH/2+(((sin_int(bb->angle)>>4)*BUBBLE_WIDTH/2)>>10), 1475 rb->lcd_drawline(SHOTX+BUBBLE_WIDTH/2+(((fp14_sin(bb->angle)>>4)*BUBBLE_WIDTH/2)>>10),
1476 SHOTY+BUBBLE_HEIGHT/2-(((cos_int(bb->angle)>>4)*BUBBLE_HEIGHT/2)>>10), 1476 SHOTY+BUBBLE_HEIGHT/2-(((fp14_cos(bb->angle)>>4)*BUBBLE_HEIGHT/2)>>10),
1477 tipx, tipy); 1477 tipx, tipy);
1478 xlcd_filltriangle(tipx, tipy, 1478 xlcd_filltriangle(tipx, tipy,
1479 tipx+(((sin_int(bb->angle-135)>>4)*BUBBLE_WIDTH/3)>>10), 1479 tipx+(((fp14_sin(bb->angle-135)>>4)*BUBBLE_WIDTH/3)>>10),
1480 tipy-(((cos_int(bb->angle-135)>>4)*BUBBLE_HEIGHT/3)>>10), 1480 tipy-(((fp14_cos(bb->angle-135)>>4)*BUBBLE_HEIGHT/3)>>10),
1481 tipx+(((sin_int(bb->angle+135)>>4)*BUBBLE_WIDTH/3)>>10), 1481 tipx+(((fp14_sin(bb->angle+135)>>4)*BUBBLE_WIDTH/3)>>10),
1482 tipy-(((cos_int(bb->angle+135)>>4)*BUBBLE_HEIGHT/3)>>10)); 1482 tipy-(((fp14_cos(bb->angle+135)>>4)*BUBBLE_HEIGHT/3)>>10));
1483 1483
1484 /* draw text */ 1484 /* draw text */
1485 rb->lcd_getstringsize(level, &w, &h); 1485 rb->lcd_getstringsize(level, &w, &h);
@@ -1524,8 +1524,8 @@ static int bubbles_fire(struct game_context* bb) {
1524 1524
1525 /* get current bubble */ 1525 /* get current bubble */
1526 bubblecur = bb->queue[bb->nextinq]; 1526 bubblecur = bb->queue[bb->nextinq];
1527 shotxinc = ((sin_int(bb->angle)>>4)*BUBBLE_WIDTH)/3; 1527 shotxinc = ((fp14_sin(bb->angle)>>4)*BUBBLE_WIDTH)/3;
1528 shotyinc = ((-1*(cos_int(bb->angle)>>4))*BUBBLE_HEIGHT)/3; 1528 shotyinc = ((-1*(fp14_cos(bb->angle)>>4))*BUBBLE_HEIGHT)/3;
1529 shotxofs = shotyofs = 0; 1529 shotxofs = shotyofs = 0;
1530 1530
1531 /* advance the queue */ 1531 /* advance the queue */
diff --git a/apps/plugins/clock/clock_draw_analog.c b/apps/plugins/clock/clock_draw_analog.c
index c41ec3b24c..9efe3623a3 100644
--- a/apps/plugins/clock/clock_draw_analog.c
+++ b/apps/plugins/clock/clock_draw_analog.c
@@ -41,11 +41,11 @@ void polar_to_cartesian(int a, int r, int* x, int* y)
41{ 41{
42#if CONFIG_LCD == LCD_SSD1815 42#if CONFIG_LCD == LCD_SSD1815
43 /* Correct non-square pixel aspect of archos recorder LCD */ 43 /* Correct non-square pixel aspect of archos recorder LCD */
44 *x = (sin_int(a) * 5 / 4 * r) >> 14; 44 *x = (fp14_sin(a) * 5 / 4 * r) >> 14;
45#else 45#else
46 *x = (sin_int(a) * r) >> 14; 46 *x = (fp14_sin(a) * r) >> 14;
47#endif 47#endif
48 *y = (sin_int(a-90) * r) >> 14; 48 *y = (fp14_sin(a-90) * r) >> 14;
49} 49}
50 50
51void polar_to_cartesian_screen_centered(struct screen * display, 51void polar_to_cartesian_screen_centered(struct screen * display,
diff --git a/apps/plugins/cube.c b/apps/plugins/cube.c
index 55219e5a5e..c770214700 100644
--- a/apps/plugins/cube.c
+++ b/apps/plugins/cube.c
@@ -433,12 +433,12 @@ static void cube_rotate(int xa, int ya, int za)
433 /* Just to prevent unnecessary lookups */ 433 /* Just to prevent unnecessary lookups */
434 long sxa, cxa, sya, cya, sza, cza; 434 long sxa, cxa, sya, cya, sza, cza;
435 435
436 sxa = sin_int(xa); 436 sxa = fp14_sin(xa);
437 cxa = cos_int(xa); 437 cxa = fp14_cos(xa);
438 sya = sin_int(ya); 438 sya = fp14_sin(ya);
439 cya = cos_int(ya); 439 cya = fp14_cos(ya);
440 sza = sin_int(za); 440 sza = fp14_sin(za);
441 cza = cos_int(za); 441 cza = fp14_cos(za);
442 442
443 /* calculate overall translation matrix */ 443 /* calculate overall translation matrix */
444 matrice[0][0] = (cza * cya) >> 14; 444 matrice[0][0] = (cza * cya) >> 14;
diff --git a/apps/plugins/lib/fixedpoint.c b/apps/plugins/lib/fixedpoint.c
index 0ae2cded69..352e246673 100644
--- a/apps/plugins/lib/fixedpoint.c
+++ b/apps/plugins/lib/fixedpoint.c
@@ -1,238 +1 @@
1/*************************************************************************** #include "../../fixedpoint.c"
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 <inttypes.h>
25#include "plugin.h"
26#include "fixedpoint.h"
27
28/* Inverse gain of circular cordic rotation in s0.31 format. */
29static const long cordic_circular_gain = 0xb2458939; /* 0.607252929 */
30
31/* Table of values of atan(2^-i) in 0.32 format fractions of pi where pi = 0xffffffff / 2 */
32static const unsigned long atan_table[] = {
33 0x1fffffff, /* +0.785398163 (or pi/4) */
34 0x12e4051d, /* +0.463647609 */
35 0x09fb385b, /* +0.244978663 */
36 0x051111d4, /* +0.124354995 */
37 0x028b0d43, /* +0.062418810 */
38 0x0145d7e1, /* +0.031239833 */
39 0x00a2f61e, /* +0.015623729 */
40 0x00517c55, /* +0.007812341 */
41 0x0028be53, /* +0.003906230 */
42 0x00145f2e, /* +0.001953123 */
43 0x000a2f98, /* +0.000976562 */
44 0x000517cc, /* +0.000488281 */
45 0x00028be6, /* +0.000244141 */
46 0x000145f3, /* +0.000122070 */
47 0x0000a2f9, /* +0.000061035 */
48 0x0000517c, /* +0.000030518 */
49 0x000028be, /* +0.000015259 */
50 0x0000145f, /* +0.000007629 */
51 0x00000a2f, /* +0.000003815 */
52 0x00000517, /* +0.000001907 */
53 0x0000028b, /* +0.000000954 */
54 0x00000145, /* +0.000000477 */
55 0x000000a2, /* +0.000000238 */
56 0x00000051, /* +0.000000119 */
57 0x00000028, /* +0.000000060 */
58 0x00000014, /* +0.000000030 */
59 0x0000000a, /* +0.000000015 */
60 0x00000005, /* +0.000000007 */
61 0x00000002, /* +0.000000004 */
62 0x00000001, /* +0.000000002 */
63 0x00000000, /* +0.000000001 */
64 0x00000000, /* +0.000000000 */
65};
66
67/* Precalculated sine and cosine * 16384 (2^14) (fixed point 18.14) */
68static const short sin_table[91] =
69{
70 0, 285, 571, 857, 1142, 1427, 1712, 1996, 2280, 2563,
71 2845, 3126, 3406, 3685, 3963, 4240, 4516, 4790, 5062, 5334,
72 5603, 5871, 6137, 6401, 6663, 6924, 7182, 7438, 7691, 7943,
73 8191, 8438, 8682, 8923, 9161, 9397, 9630, 9860, 10086, 10310,
74 10531, 10748, 10963, 11173, 11381, 11585, 11785, 11982, 12175, 12365,
75 12550, 12732, 12910, 13084, 13254, 13420, 13582, 13740, 13894, 14043,
76 14188, 14329, 14466, 14598, 14725, 14848, 14967, 15081, 15190, 15295,
77 15395, 15491, 15582, 15668, 15749, 15825, 15897, 15964, 16025, 16082,
78 16135, 16182, 16224, 16261, 16294, 16321, 16344, 16361, 16374, 16381,
79 16384
80};
81
82/**
83 * Implements sin and cos using CORDIC rotation.
84 *
85 * @param phase has range from 0 to 0xffffffff, representing 0 and
86 * 2*pi respectively.
87 * @param cos return address for cos
88 * @return sin of phase, value is a signed value from LONG_MIN to LONG_MAX,
89 * representing -1 and 1 respectively.
90 */
91long fsincos(unsigned long phase, long *cos)
92{
93 int32_t x, x1, y, y1;
94 unsigned long z, z1;
95 int i;
96
97 /* Setup initial vector */
98 x = cordic_circular_gain;
99 y = 0;
100 z = phase;
101
102 /* The phase has to be somewhere between 0..pi for this to work right */
103 if (z < 0xffffffff / 4) {
104 /* z in first quadrant, z += pi/2 to correct */
105 x = -x;
106 z += 0xffffffff / 4;
107 } else if (z < 3 * (0xffffffff / 4)) {
108 /* z in third quadrant, z -= pi/2 to correct */
109 z -= 0xffffffff / 4;
110 } else {
111 /* z in fourth quadrant, z -= 3pi/2 to correct */
112 x = -x;
113 z -= 3 * (0xffffffff / 4);
114 }
115
116 /* Each iteration adds roughly 1-bit of extra precision */
117 for (i = 0; i < 31; i++) {
118 x1 = x >> i;
119 y1 = y >> i;
120 z1 = atan_table[i];
121
122 /* Decided which direction to rotate vector. Pivot point is pi/2 */
123 if (z >= 0xffffffff / 4) {
124 x -= y1;
125 y += x1;
126 z -= z1;
127 } else {
128 x += y1;
129 y -= x1;
130 z += z1;
131 }
132 }
133
134 if (cos)
135 *cos = x;
136
137 return y;
138}
139
140/**
141 * Fixed point square root via Newton-Raphson.
142 * @param a square root argument.
143 * @param fracbits specifies number of fractional bits in argument.
144 * @return Square root of argument in same fixed point format as input.
145 */
146long fsqrt(long a, unsigned int fracbits)
147{
148 long b = a/2 + BIT_N(fracbits); /* initial approximation */
149 unsigned n;
150 const unsigned iterations = 4;
151
152 for (n = 0; n < iterations; ++n)
153 b = (b + (long)(((long long)(a) << fracbits)/b))/2;
154
155 return b;
156}
157
158/**
159 * Fixed point sinus using a lookup table
160 * don't forget to divide the result by 16384 to get the actual sinus value
161 * @param val sinus argument in degree
162 * @return sin(val)*16384
163 */
164long sin_int(int val)
165{
166 val = (val+360)%360;
167 if (val < 181)
168 {
169 if (val < 91)/* phase 0-90 degree */
170 return (long)sin_table[val];
171 else/* phase 91-180 degree */
172 return (long)sin_table[180-val];
173 }
174 else
175 {
176 if (val < 271)/* phase 181-270 degree */
177 return -(long)sin_table[val-180];
178 else/* phase 270-359 degree */
179 return -(long)sin_table[360-val];
180 }
181 return 0;
182}
183
184/**
185 * Fixed point cosinus using a lookup table
186 * don't forget to divide the result by 16384 to get the actual cosinus value
187 * @param val sinus argument in degree
188 * @return cos(val)*16384
189 */
190long cos_int(int val)
191{
192 val = (val+360)%360;
193 if (val < 181)
194 {
195 if (val < 91)/* phase 0-90 degree */
196 return (long)sin_table[90-val];
197 else/* phase 91-180 degree */
198 return -(long)sin_table[val-90];
199 }
200 else
201 {
202 if (val < 271)/* phase 181-270 degree */
203 return -(long)sin_table[270-val];
204 else/* phase 270-359 degree */
205 return (long)sin_table[val-270];
206 }
207 return 0;
208}
209
210/**
211 * Fixed-point natural log
212 * taken from http://www.quinapalus.com/efunc.html
213 * "The code assumes integers are at least 32 bits long. The (positive)
214 * argument and the result of the function are both expressed as fixed-point
215 * values with 16 fractional bits, although intermediates are kept with 28
216 * bits of precision to avoid loss of accuracy during shifts."
217 */
218
219long flog(int x) {
220 long t,y;
221
222 y=0xa65af;
223 if(x<0x00008000) x<<=16, y-=0xb1721;
224 if(x<0x00800000) x<<= 8, y-=0x58b91;
225 if(x<0x08000000) x<<= 4, y-=0x2c5c8;
226 if(x<0x20000000) x<<= 2, y-=0x162e4;
227 if(x<0x40000000) x<<= 1, y-=0x0b172;
228 t=x+(x>>1); if((t&0x80000000)==0) x=t,y-=0x067cd;
229 t=x+(x>>2); if((t&0x80000000)==0) x=t,y-=0x03920;
230 t=x+(x>>3); if((t&0x80000000)==0) x=t,y-=0x01e27;
231 t=x+(x>>4); if((t&0x80000000)==0) x=t,y-=0x00f85;
232 t=x+(x>>5); if((t&0x80000000)==0) x=t,y-=0x007e1;
233 t=x+(x>>6); if((t&0x80000000)==0) x=t,y-=0x003f8;
234 t=x+(x>>7); if((t&0x80000000)==0) x=t,y-=0x001fe;
235 x=0x80000000-x;
236 y-=x>>15;
237 return y;
238}
diff --git a/apps/plugins/lib/fixedpoint.h b/apps/plugins/lib/fixedpoint.h
index dfabbad8dc..ef50dd0085 100644
--- a/apps/plugins/lib/fixedpoint.h
+++ b/apps/plugins/lib/fixedpoint.h
@@ -21,11 +21,44 @@
21 * 21 *
22 ****************************************************************************/ 22 ****************************************************************************/
23 23
24long fsincos(unsigned long phase, long *cos); 24/** PLUGINS - FIXED POINT MATH ROUTINES - USAGE
25long fsqrt(long a, unsigned int fracbits); 25 *
26long cos_int(int val); 26 * - x and y arguments are fixed point integers
27long sin_int(int val); 27 * - fracbits is the number of fractional bits in the argument(s)
28long flog(int x); 28 * - functions return long fixed point integers with the specified number
29 * of fractional bits unless otherwise specified
30 *
31 * Calculate sin and cos of an angle:
32 * fp_sincos(phase, *cos)
33 * where phase is a 32 bit unsigned integer with 0 representing 0
34 * and 0xFFFFFFFF representing 2*pi, and *cos is the address to
35 * a long signed integer. Value returned is a long signed integer
36 * from LONG_MIN to LONG_MAX, representing -1 to 1 respectively.
37 * That is, value is a fixed point integer with 31 fractional bits.
38 *
39 * Take square root of a fixed point number:
40 * fp_sqrt(x, fracbits)
41 *
42 * Calculate sin or cos of an angle (very fast, from a table):
43 * fp14_sin(angle)
44 * fp14_cos(angle)
45 * where angle is a non-fixed point integer in degrees. Value
46 * returned is a fixed point integer with 14 fractional bits.
47 *
48 * Calculate the natural log of a positive fixed point integer
49 * fp16_log(x)
50 * where x and the value returned are fixed point integers
51 * with 16 fractional bits.
52 */
53
54#ifndef _FIXEDPOINT_H_PLUGINS
55#define _FIXEDPOINT_H_PLUGINS
56
57long fp_sincos(unsigned long phase, long *cos);
58long fp_sqrt(long a, unsigned int fracbits);
59long fp14_cos(int val);
60long fp14_sin(int val);
61long fp16_log(int x);
29 62
30/* fast unsigned multiplication (16x16bit->32bit or 32x32bit->32bit, 63/* fast unsigned multiplication (16x16bit->32bit or 32x32bit->32bit,
31 * whichever is faster for the architecture) */ 64 * whichever is faster for the architecture) */
@@ -34,3 +67,5 @@ long flog(int x);
34#else /* SH1, coldfire */ 67#else /* SH1, coldfire */
35#define FMULU(a, b) ((uint32_t) (((uint16_t) (a)) * ((uint16_t) (b)))) 68#define FMULU(a, b) ((uint32_t) (((uint16_t) (a)) * ((uint16_t) (b))))
36#endif 69#endif
70
71#endif
diff --git a/apps/plugins/plasma.c b/apps/plugins/plasma.c
index 2a3e43e6b8..00287eb0b8 100644
--- a/apps/plugins/plasma.c
+++ b/apps/plugins/plasma.c
@@ -198,7 +198,7 @@ static void wave_table_generate(void)
198 for (i=0;i<256;++i) 198 for (i=0;i<256;++i)
199 { 199 {
200 wave_array[i] = (unsigned char)((WAV_AMP 200 wave_array[i] = (unsigned char)((WAV_AMP
201 * (sin_int((i * 360 * plasma_frequency) / 256))) / 16384); 201 * (fp14_sin((i * 360 * plasma_frequency) / 256))) / 16384);
202 } 202 }
203} 203}
204 204
diff --git a/apps/plugins/vu_meter.c b/apps/plugins/vu_meter.c
index 16aac3a011..74c3b1cf97 100644
--- a/apps/plugins/vu_meter.c
+++ b/apps/plugins/vu_meter.c
@@ -415,7 +415,7 @@ void calc_scales(void)
415 for (i=1; i <= half_width; i++) 415 for (i=1; i <= half_width; i++)
416 { 416 {
417 /* analog scale */ 417 /* analog scale */
418 y = (half_width/5)*flog(i*fx_log_factor); 418 y = (half_width/5)*fp16_log(i*fx_log_factor);
419 419
420 /* better way of checking for negative values? */ 420 /* better way of checking for negative values? */
421 z = y>>16; 421 z = y>>16;
@@ -431,7 +431,7 @@ void calc_scales(void)
431 k = nh2 - ( j * j ); 431 k = nh2 - ( j * j );
432 432
433 /* fsqrt+1 seems to give a closer approximation */ 433 /* fsqrt+1 seems to give a closer approximation */
434 y_values[i-1] = LCD_HEIGHT - (fsqrt(k, 16)>>8) - 1; 434 y_values[i-1] = LCD_HEIGHT - (fp_sqrt(k, 16)>>8) - 1;
435 rb->yield(); 435 rb->yield();
436 } 436 }
437} 437}
diff --git a/apps/replaygain.c b/apps/replaygain.c
index 90944f91d0..b398afc294 100644
--- a/apps/replaygain.c
+++ b/apps/replaygain.c
@@ -30,188 +30,11 @@
30#include "metadata.h" 30#include "metadata.h"
31#include "debug.h" 31#include "debug.h"
32#include "replaygain.h" 32#include "replaygain.h"
33 33#include "fixedpoint.h"
34/* The fixed point math routines (with the exception of fp_atof) are based
35 * on oMathFP by Dan Carter (http://orbisstudios.com).
36 */
37
38/* 12 bits of precision gives fairly accurate result, but still allows a
39 * compact implementation. The math code supports up to 13...
40 */
41 34
42#define FP_BITS (12) 35#define FP_BITS (12)
43#define FP_MASK ((1 << FP_BITS) - 1)
44#define FP_ONE (1 << FP_BITS) 36#define FP_ONE (1 << FP_BITS)
45#define FP_TWO (2 << FP_BITS)
46#define FP_HALF (1 << (FP_BITS - 1))
47#define FP_LN2 ( 45426 >> (16 - FP_BITS))
48#define FP_LN2_INV ( 94548 >> (16 - FP_BITS))
49#define FP_EXP_ZERO ( 10922 >> (16 - FP_BITS))
50#define FP_EXP_ONE ( -182 >> (16 - FP_BITS))
51#define FP_EXP_TWO ( 4 >> (16 - FP_BITS))
52#define FP_INF (0x7fffffff)
53#define FP_LN10 (150902 >> (16 - FP_BITS))
54
55#define FP_MAX_DIGITS (4)
56#define FP_MAX_DIGITS_INT (10000)
57
58#define FP_FAST_MUL_DIV
59
60#ifdef FP_FAST_MUL_DIV
61
62/* These macros can easily overflow, but they are good enough for our uses,
63 * and saves some code.
64 */
65#define fp_mul(x, y) (((x) * (y)) >> FP_BITS)
66#define fp_div(x, y) (((x) << FP_BITS) / (y))
67
68#else
69
70static long fp_mul(long x, long y)
71{
72 long x_neg = 0;
73 long y_neg = 0;
74 long rc;
75
76 if ((x == 0) || (y == 0))
77 {
78 return 0;
79 }
80
81 if (x < 0)
82 {
83 x_neg = 1;
84 x = -x;
85 }
86
87 if (y < 0)
88 {
89 y_neg = 1;
90 y = -y;
91 }
92
93 rc = (((x >> FP_BITS) * (y >> FP_BITS)) << FP_BITS)
94 + (((x & FP_MASK) * (y & FP_MASK)) >> FP_BITS)
95 + ((x & FP_MASK) * (y >> FP_BITS))
96 + ((x >> FP_BITS) * (y & FP_MASK));
97
98 if ((x_neg ^ y_neg) == 1)
99 {
100 rc = -rc;
101 }
102
103 return rc;
104}
105
106static long fp_div(long x, long y)
107{
108 long x_neg = 0;
109 long y_neg = 0;
110 long shifty;
111 long rc;
112 int msb = 0;
113 int lsb = 0;
114
115 if (x == 0)
116 {
117 return 0;
118 }
119
120 if (y == 0)
121 {
122 return (x < 0) ? -FP_INF : FP_INF;
123 }
124
125 if (x < 0)
126 {
127 x_neg = 1;
128 x = -x;
129 }
130
131 if (y < 0)
132 {
133 y_neg = 1;
134 y = -y;
135 }
136
137 while ((x & BIT_N(30 - msb)) == 0)
138 {
139 msb++;
140 }
141
142 while ((y & BIT_N(lsb)) == 0)
143 {
144 lsb++;
145 }
146
147 shifty = FP_BITS - (msb + lsb);
148 rc = ((x << msb) / (y >> lsb));
149 37
150 if (shifty > 0)
151 {
152 rc <<= shifty;
153 }
154 else
155 {
156 rc >>= -shifty;
157 }
158
159 if ((x_neg ^ y_neg) == 1)
160 {
161 rc = -rc;
162 }
163
164 return rc;
165}
166
167#endif /* FP_FAST_MUL_DIV */
168
169static long fp_exp(long x)
170{
171 long k;
172 long z;
173 long R;
174 long xp;
175
176 if (x == 0)
177 {
178 return FP_ONE;
179 }
180
181 k = (fp_mul(abs(x), FP_LN2_INV) + FP_HALF) & ~FP_MASK;
182
183 if (x < 0)
184 {
185 k = -k;
186 }
187
188 x -= fp_mul(k, FP_LN2);
189 z = fp_mul(x, x);
190 R = FP_TWO + fp_mul(z, FP_EXP_ZERO + fp_mul(z, FP_EXP_ONE
191 + fp_mul(z, FP_EXP_TWO)));
192 xp = FP_ONE + fp_div(fp_mul(FP_TWO, x), R - x);
193
194 if (k < 0)
195 {
196 k = FP_ONE >> (-k >> FP_BITS);
197 }
198 else
199 {
200 k = FP_ONE << (k >> FP_BITS);
201 }
202
203 return fp_mul(k, xp);
204}
205
206static long fp_exp10(long x)
207{
208 if (x == 0)
209 {
210 return FP_ONE;
211 }
212
213 return fp_exp(fp_mul(FP_LN10, x));
214}
215 38
216static long fp_atof(const char* s, int precision) 39static long fp_atof(const char* s, int precision)
217{ 40{
@@ -300,7 +123,7 @@ static long convert_gain(long gain)
300 gain = 17 * FP_ONE; 123 gain = 17 * FP_ONE;
301 } 124 }
302 125
303 gain = fp_exp10(gain / 20) << (24 - FP_BITS); 126 gain = fp_factor(gain, FP_BITS) << (24 - FP_BITS);
304 127
305 return gain; 128 return gain;
306} 129}