summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorMaurus Cuelenaere <mcuelenaere@gmail.com>2009-07-04 12:26:45 +0000
committerMaurus Cuelenaere <mcuelenaere@gmail.com>2009-07-04 12:26:45 +0000
commit4710a3280025b0ba8ffb6e8183578a5df48257fa (patch)
tree3dd82b90ab668c18109d0885cd3114449efaddf1
parent69c73e8bd6c8c6ac79c6538cb0ad4686b9d1d920 (diff)
downloadrockbox-4710a3280025b0ba8ffb6e8183578a5df48257fa.tar.gz
rockbox-4710a3280025b0ba8ffb6e8183578a5df48257fa.zip
Consolidate all fixed point math routines in one library (FS#10400) by Jeffrey Goode
git-svn-id: svn://svn.rockbox.org/rockbox/trunk@21633 a1c6a512-1295-4272-9138-f99709370657
-rw-r--r--apps/SOURCES1
-rw-r--r--apps/codecs/adx.c119
-rw-r--r--apps/codecs/lib/SOURCES2
-rw-r--r--apps/codecs/lib/fixedpoint.h126
-rw-r--r--apps/codecs/spc.c1
-rw-r--r--apps/dsp.c1
-rw-r--r--apps/dsp.h80
-rw-r--r--apps/eq.c97
-rw-r--r--apps/eq.h1
-rw-r--r--apps/fixedpoint.c440
-rw-r--r--apps/fixedpoint.h197
-rw-r--r--apps/plugins/lib/SOURCES2
-rw-r--r--apps/plugins/lib/fixedpoint.c238
-rw-r--r--apps/replaygain.c181
14 files changed, 773 insertions, 713 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..e23b3d4f80 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 "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{
diff --git a/apps/codecs/lib/SOURCES b/apps/codecs/lib/SOURCES
index cbb8e60372..a1730f656a 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 3../../fixedpoint.c
4 4
5mdct2.c 5mdct2.c
6#ifdef CPU_ARM 6#ifdef CPU_ARM
diff --git a/apps/codecs/lib/fixedpoint.h b/apps/codecs/lib/fixedpoint.h
new file mode 100644
index 0000000000..54ece27080
--- /dev/null
+++ b/apps/codecs/lib/fixedpoint.h
@@ -0,0 +1,126 @@
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#ifndef _FIXEDPOINT_H
25#define _FIXEDPOINT_H
26
27#include <inttypes.h>
28
29/** TAKEN FROM apps/dsp.h */
30/* A bunch of fixed point assembler helper macros */
31#if defined(CPU_COLDFIRE)
32/* These macros use the Coldfire EMAC extension and need the MACSR flags set
33 * to fractional mode with no rounding.
34 */
35
36/* Multiply two S.31 fractional integers and return the sign bit and the
37 * 31 most significant bits of the result.
38 */
39#define FRACMUL(x, y) \
40({ \
41 long t; \
42 asm ("mac.l %[a], %[b], %%acc0\n\t" \
43 "movclr.l %%acc0, %[t]\n\t" \
44 : [t] "=r" (t) : [a] "r" (x), [b] "r" (y)); \
45 t; \
46})
47
48/* Multiply two S.31 fractional integers, and return the 32 most significant
49 * bits after a shift left by the constant z. NOTE: Only works for shifts of
50 * 1 to 8 on Coldfire!
51 */
52#define FRACMUL_SHL(x, y, z) \
53({ \
54 long t, t2; \
55 asm ("mac.l %[a], %[b], %%acc0\n\t" \
56 "moveq.l %[d], %[t]\n\t" \
57 "move.l %%accext01, %[t2]\n\t" \
58 "and.l %[mask], %[t2]\n\t" \
59 "lsr.l %[t], %[t2]\n\t" \
60 "movclr.l %%acc0, %[t]\n\t" \
61 "asl.l %[c], %[t]\n\t" \
62 "or.l %[t2], %[t]\n\t" \
63 : [t] "=&d" (t), [t2] "=&d" (t2) \
64 : [a] "r" (x), [b] "r" (y), [mask] "d" (0xff), \
65 [c] "i" ((z)), [d] "i" (8 - (z))); \
66 t; \
67})
68
69#elif defined(CPU_ARM)
70
71/* Multiply two S.31 fractional integers and return the sign bit and the
72 * 31 most significant bits of the result.
73 */
74#define FRACMUL(x, y) \
75({ \
76 long t, t2; \
77 asm ("smull %[t], %[t2], %[a], %[b]\n\t" \
78 "mov %[t2], %[t2], asl #1\n\t" \
79 "orr %[t], %[t2], %[t], lsr #31\n\t" \
80 : [t] "=&r" (t), [t2] "=&r" (t2) \
81 : [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.
87 */
88#define FRACMUL_SHL(x, y, z) \
89({ \
90 long t, t2; \
91 asm ("smull %[t], %[t2], %[a], %[b]\n\t" \
92 "mov %[t2], %[t2], asl %[c]\n\t" \
93 "orr %[t], %[t2], %[t], lsr %[d]\n\t" \
94 : [t] "=&r" (t), [t2] "=&r" (t2) \
95 : [a] "r" (x), [b] "r" (y), \
96 [c] "M" ((z) + 1), [d] "M" (31 - (z))); \
97 t; \
98})
99
100#else
101
102#define FRACMUL(x, y) (long) (((((long long) (x)) * ((long long) (y))) >> 31))
103#define FRACMUL_SHL(x, y, z) \
104((long)(((((long long) (x)) * ((long long) (y))) >> (31 - (z)))))
105
106#endif
107
108#define DIV64(x, y, z) (long)(((long long)(x) << (z))/(y))
109
110
111/** TAKEN FROM ORIGINAL fixedpoint.h */
112/* fast unsigned multiplication (16x16bit->32bit or 32x32bit->32bit,
113 * whichever is faster for the architecture) */
114#ifdef CPU_ARM
115#define FMULU(a, b) ((uint32_t) (((uint32_t) (a)) * ((uint32_t) (b))))
116#else /* SH1, coldfire */
117#define FMULU(a, b) ((uint32_t) (((uint16_t) (a)) * ((uint16_t) (b))))
118#endif
119
120long fsincos(unsigned long phase, long *cos);
121long fsqrt(long a, unsigned int fracbits);
122long cos_int(int val);
123long sin_int(int val);
124long flog(int x);
125
126#endif
diff --git a/apps/codecs/spc.c b/apps/codecs/spc.c
index 6ceb704c7c..d5313bfa47 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 "fixedpoint.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..66469304b0 100644
--- a/apps/dsp.c
+++ b/apps/dsp.c
@@ -33,6 +33,7 @@
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"
36 37
37/* 16-bit samples are scaled based on these constants. The shift should be 38/* 16-bit samples are scaled based on these constants. The shift should be
38 * no more than 15. 39 * no more than 15.
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..7b7fba341d 100644
--- a/apps/eq.c
+++ b/apps/eq.c
@@ -21,105 +21,10 @@
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 "eq.h" 25#include "eq.h"
26#include "replaygain.h" 26#include "replaygain.h"
27 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/**
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/** 28/**
124 * Calculate first order shelving filter. Filter is not directly usable by the 29 * Calculate first order shelving filter. Filter is not directly usable by the
125 * eq_filter() function. 30 * eq_filter() function.
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..b65070e348
--- /dev/null
+++ b/apps/fixedpoint.c
@@ -0,0 +1,440 @@
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
28#ifndef BIT_N
29#define BIT_N(n) (1U << (n))
30#endif
31
32/** TAKEN FROM ORIGINAL fixedpoint.h */
33/* Inverse gain of circular cordic rotation in s0.31 format. */
34static const long cordic_circular_gain = 0xb2458939; /* 0.607252929 */
35
36/* Table of values of atan(2^-i) in 0.32 format fractions of pi where pi = 0xffffffff / 2 */
37static 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 */
70};
71
72/* Precalculated sine and cosine * 16384 (2^14) (fixed point 18.14) */
73static const short sin_table[91] =
74{
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,
84 16384
85};
86
87/**
88 * Implements sin and cos using CORDIC rotation.
89 *
90 * @param phase has range from 0 to 0xffffffff, representing 0 and
91 * 2*pi respectively.
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.
95 */
96long fsincos(unsigned long phase, long *cos)
97{
98 int32_t x, x1, y, y1;
99 unsigned long z, z1;
100 int i;
101
102 /* Setup initial vector */
103 x = cordic_circular_gain;
104 y = 0;
105 z = phase;
106
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 */
110 x = -x;
111 z += 0xffffffff / 4;
112 } else if (z < 3 * (0xffffffff / 4)) {
113 /* z in third quadrant, z -= pi/2 to correct */
114 z -= 0xffffffff / 4;
115 } else {
116 /* z in fourth quadrant, z -= 3pi/2 to correct */
117 x = -x;
118 z -= 3 * (0xffffffff / 4);
119 }
120
121 /* Each iteration adds roughly 1-bit of extra precision */
122 for (i = 0; i < 31; i++) {
123 x1 = x >> i;
124 y1 = y >> i;
125 z1 = atan_table[i];
126
127 /* Decided which direction to rotate vector. Pivot point is pi/2 */
128 if (z >= 0xffffffff / 4) {
129 x -= y1;
130 y += x1;
131 z -= z1;
132 } else {
133 x += y1;
134 y -= x1;
135 z += z1;
136 }
137 }
138
139 if (cos)
140 *cos = x;
141
142 return y;
143}
144
145/**
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.
150 *
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.
154 */
155long fsqrt(long x, unsigned int fracbits)
156{
157 long b = x/2 + BIT_N(fracbits); /* initial approximation */
158 long c;
159 unsigned n;
160 const unsigned iterations = 8;
161
162 for (n = 0; n < iterations; ++n)
163 {
164 c = DIV64(x, b, fracbits);
165 if (c == b) break;
166 b = (b + c)/2;
167 }
168
169 return b;
170}
171
172/**
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
177 */
178long sin_int(int val)
179{
180 val = (val+360)%360;
181 if (val < 181)
182 {
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];
187 }
188 else
189 {
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];
194 }
195 return 0;
196}
197
198/**
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
203 */
204long cos_int(int val)
205{
206 val = (val+360)%360;
207 if (val < 181)
208 {
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];
213 }
214 else
215 {
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];
220 }
221 return 0;
222}
223
224/**
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."
231 */
232
233long flog(int x) {
234 long t,y;
235
236 y=0xa65af;
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;
249 x=0x80000000-x;
250 y-=x>>15;
251 return y;
252}
253
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.
257 */
258
259#define FP_MUL64(x, y) (((x) * (y)) >> (fracbits))
260#define FP_DIV64(x, y) (((x) << (fracbits)) / (y))
261
262static long long fp_exp10(long long x, unsigned int fracbits);
263static long long fp_log10(long long n, unsigned int fracbits);
264
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) */
273
274#define TOL_BITS 2 /* log calculation tolerance */
275
276
277/* The fpexp10 fixed point math routine is based
278 * on oMathFP by Dan Carter (http://orbisstudios.com).
279 */
280
281/** FIXED POINT EXP10
282 * Return 10^x as FP integer. Argument is FP integer.
283 */
284static long long fp_exp10(long long x, unsigned int fracbits)
285{
286 long long k;
287 long long z;
288 long long R;
289 long long xp;
290
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));
302
303 /* exp(0) = 1 */
304 if (x == 0)
305 {
306 return fp_one;
307 }
308
309 /* convert from base 10 to base e */
310 x = FP_MUL64(x, fp_ln10);
311
312 /* calculate exp(x) */
313 k = (FP_MUL64(abs(x), fp_ln2_inv) + fp_half) & ~fp_mask;
314
315 if (x < 0)
316 {
317 k = -k;
318 }
319
320 x -= FP_MUL64(k, fp_ln2);
321 z = FP_MUL64(x, x);
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);
325
326 if (k < 0)
327 {
328 k = fp_one >> (-k >> fracbits);
329 }
330 else
331 {
332 k = fp_one << (k >> fracbits);
333 }
334
335 return FP_MUL64(k, xp);
336}
337
338
339/** FIXED POINT LOG10
340 * Return log10(x) as FP integer. Argument is FP integer.
341 */
342static long long fp_log10(long long n, unsigned int fracbits)
343{
344 /* Calculate log2 of argument */
345
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));
350
351 if (n <=0) return FP_NEGINF;
352 log2 = 0;
353
354 /* integer part */
355 while (n < fp_one)
356 {
357 log2 -= fp_one;
358 n <<= 1;
359 }
360 while (n >= fp_two)
361 {
362 log2 += fp_one;
363 n >>= 1;
364 }
365
366 /* fractional part */
367 frac = fp_one;
368 while (frac > tolerance)
369 {
370 frac >>= 1;
371 n = FP_MUL64(n, n);
372 if (n >= fp_two)
373 {
374 n >>= 1;
375 log2 += frac;
376 }
377 }
378
379 /* convert log2 to log10 */
380 return FP_MUL64(log2, (FP28_LOG10OF2 >> (28 - fracbits)));
381}
382
383
384/** CONVERT FACTOR TO DECIBELS */
385long fp_decibels(unsigned long factor, unsigned int fracbits)
386{
387 long long decibels;
388 long long f = (long long)factor;
389 bool neg;
390
391 /* keep factor in signed long range */
392 if (f >= (1LL << 31))
393 f = (1LL << 31) - 1;
394
395 /* decibels = 20 * log10(factor) */
396 decibels = FP_MUL64((20LL << fracbits), fp_log10(f, fracbits));
397
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;
403
404 return neg ? (long)-decibels : (long)decibels;
405}
406
407
408/** CONVERT DECIBELS TO FACTOR */
409long fp_factor(long decibels, unsigned int fracbits)
410{
411 bool neg;
412 long long factor;
413 long long db = (long long)decibels;
414
415 /* if decibels is 0, factor is 1 */
416 if (db == 0)
417 return (1L << fracbits);
418
419 /* calculate for positive decibels only */
420 if ((neg = (db < 0)))
421 db = -db;
422
423 /* factor = 10 ^ (decibels / 20) */
424 factor = fp_exp10(FP_DIV64(db, (20LL << fracbits)), fracbits);
425
426 /* keep result in signed long range, return 0 if very small */
427 if (factor >= (1LL << 31))
428 {
429 if (neg)
430 return 0;
431 else
432 return FP_INF;
433 }
434
435 /* if negative argument, factor is 1 / result */
436 if (neg)
437 factor = FP_DIV64((1LL << fracbits), factor);
438
439 return (long)factor;
440}
diff --git a/apps/fixedpoint.h b/apps/fixedpoint.h
new file mode 100644
index 0000000000..a3ca6ee6ed
--- /dev/null
+++ b/apps/fixedpoint.h
@@ -0,0 +1,197 @@
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/** 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 * Shortcut: Multiply two fixed point numbers with 31 fractional bits:
35 * fp31_mul(x, y)
36 *
37 * Shortcut: Multiply two fixed point numbers with 31 fractional bits,
38 * then shift left by z bits:
39 * fp31_mulshl(x, y, z)
40 * NOTE: z must be in the range 1-8 on Coldfire targets.
41 *
42 * Divide two fixed point numbers:
43 * fp_div(x, y, fracbits)
44 *
45 * Take square root of a fixed point number:
46 * fp_sqrt(x, fracbits)
47 *
48 * Calculate sin and cos of an angle:
49 * fp_sincos(phase, *cos)
50 * where phase is a 32 bit unsigned integer with 0 representing 0
51 * and 0xFFFFFFFF representing 2*pi, and *cos is the address to
52 * a long signed integer. Value returned is a long signed integer
53 * from LONG_MIN to LONG_MAX, representing -1 to 1 respectively.
54 * That is, value is a fixed point integer with 31 fractional bits.
55 *
56 * Calculate sin or cos of an angle (very fast, from a table):
57 * fp14_sin(angle)
58 * fp14_cos(angle)
59 * where angle is a non-fixed point integer in degrees. Value
60 * returned is a fixed point integer with 14 fractional bits.
61 *
62 * Calculate decibel equivalent of a gain factor:
63 * fp_decibels(factor, fracbits)
64 * where fracbits is in the range 12 to 22 (higher is better),
65 * and factor is a positive fixed point integer.
66 *
67 * Calculate factor equivalent of a decibel value:
68 * fp_factor(decibels, fracbits)
69 * where fracbits is in the range 12 to 22 (lower is better),
70 * and decibels is a fixed point integer.
71 */
72
73#ifndef _FIXEDPOINT_H
74#define _FIXEDPOINT_H
75
76#include <inttypes.h>
77
78/* Redefine function names, making sure legacy code is usable */
79#define fp31_mul(x, y) FRACMUL(x, y)
80#define fp31_mulshl(x, y, z) FRACMUL_SHL(x, y, z)
81#define fp_div(x, y, z) DIV64(x, y, z)
82#define fp_sqrt(x, y) fsqrt(x, y)
83#define fp_sincos(x, y) fsincos(x, y)
84#define fp14_sin(x) sin_int(x)
85#define fp14_cos(x) cos_int(x)
86#define fp16_log(x) flog(x)
87
88
89#define fp_mul(x, y, z) (long)((((long long)(x)) * ((long long)(y))) >> (z))
90#define DIV64(x, y, z) (long)((((long long)(x)) << (z)) / ((long long)(y)))
91
92/** TAKEN FROM apps/dsp.h */
93/* A bunch of fixed point assembler helper macros */
94#if defined(CPU_COLDFIRE)
95/* These macros use the Coldfire EMAC extension and need the MACSR flags set
96 * to fractional mode with no rounding.
97 */
98
99/* Multiply two S.31 fractional integers and return the sign bit and the
100 * 31 most significant bits of the result.
101 */
102#define FRACMUL(x, y) \
103({ \
104 long t; \
105 asm ("mac.l %[a], %[b], %%acc0\n\t" \
106 "movclr.l %%acc0, %[t]\n\t" \
107 : [t] "=r" (t) : [a] "r" (x), [b] "r" (y)); \
108 t; \
109})
110
111/* Multiply two S.31 fractional integers, and return the 32 most significant
112 * bits after a shift left by the constant z. NOTE: Only works for shifts of
113 * 1 to 8 on Coldfire!
114 */
115#define FRACMUL_SHL(x, y, z) \
116({ \
117 long t, t2; \
118 asm ("mac.l %[a], %[b], %%acc0\n\t" \
119 "moveq.l %[d], %[t]\n\t" \
120 "move.l %%accext01, %[t2]\n\t" \
121 "and.l %[mask], %[t2]\n\t" \
122 "lsr.l %[t], %[t2]\n\t" \
123 "movclr.l %%acc0, %[t]\n\t" \
124 "asl.l %[c], %[t]\n\t" \
125 "or.l %[t2], %[t]\n\t" \
126 : [t] "=&d" (t), [t2] "=&d" (t2) \
127 : [a] "r" (x), [b] "r" (y), [mask] "d" (0xff), \
128 [c] "i" ((z)), [d] "i" (8 - (z))); \
129 t; \
130})
131
132#elif defined(CPU_ARM)
133
134/* Multiply two S.31 fractional integers and return the sign bit and the
135 * 31 most significant bits of the result.
136 */
137#define FRACMUL(x, y) \
138({ \
139 long t, t2; \
140 asm ("smull %[t], %[t2], %[a], %[b]\n\t" \
141 "mov %[t2], %[t2], asl #1\n\t" \
142 "orr %[t], %[t2], %[t], lsr #31\n\t" \
143 : [t] "=&r" (t), [t2] "=&r" (t2) \
144 : [a] "r" (x), [b] "r" (y)); \
145 t; \
146})
147
148/* Multiply two S.31 fractional integers, and return the 32 most significant
149 * bits after a shift left by the constant z.
150 */
151#define FRACMUL_SHL(x, y, z) \
152({ \
153 long t, t2; \
154 asm ("smull %[t], %[t2], %[a], %[b]\n\t" \
155 "mov %[t2], %[t2], asl %[c]\n\t" \
156 "orr %[t], %[t2], %[t], lsr %[d]\n\t" \
157 : [t] "=&r" (t), [t2] "=&r" (t2) \
158 : [a] "r" (x), [b] "r" (y), \
159 [c] "M" ((z) + 1), [d] "M" (31 - (z))); \
160 t; \
161})
162
163#else
164
165#define FRACMUL(x, y) (long) (((((long long) (x)) * ((long long) (y))) >> 31))
166#define FRACMUL_SHL(x, y, z) \
167((long)(((((long long) (x)) * ((long long) (y))) >> (31 - (z)))))
168
169#endif
170
171/** TAKEN FROM ORIGINAL fixedpoint.h */
172/* fast unsigned multiplication (16x16bit->32bit or 32x32bit->32bit,
173 * whichever is faster for the architecture) */
174#ifdef CPU_ARM
175#define FMULU(a, b) ((uint32_t) (((uint32_t) (a)) * ((uint32_t) (b))))
176#else /* SH1, coldfire */
177#define FMULU(a, b) ((uint32_t) (((uint16_t) (a)) * ((uint16_t) (b))))
178#endif
179
180long fsincos(unsigned long phase, long *cos);
181long fsqrt(long x, unsigned int fracbits);
182long sin_int(int val);
183long cos_int(int val);
184long flog(int x);
185
186
187/** MODIFIED FROM replaygain.c */
188#define FP_INF (0x7fffffff)
189#define FP_NEGINF -(0x7fffffff)
190
191/* fracbits in range 12 - 22 work well. Higher is better for
192 * calculating dB, lower is better for calculating ratio.
193 */
194long fp_decibels(unsigned long factor, unsigned int fracbits);
195long fp_factor(long decibels, unsigned int fracbits);
196
197#endif
diff --git a/apps/plugins/lib/SOURCES b/apps/plugins/lib/SOURCES
index 02adb7089c..bcce3f2969 100644
--- a/apps/plugins/lib/SOURCES
+++ b/apps/plugins/lib/SOURCES
@@ -1,7 +1,7 @@
1gcc-support.c 1gcc-support.c
2jhash.c 2jhash.c
3configfile.c 3configfile.c
4fixedpoint.c 4../../fixedpoint.c
5playback_control.c 5playback_control.c
6rgb_hsv.c 6rgb_hsv.c
7buflib.c 7buflib.c
diff --git a/apps/plugins/lib/fixedpoint.c b/apps/plugins/lib/fixedpoint.c
index 0ae2cded69..e69de29bb2 100644
--- a/apps/plugins/lib/fixedpoint.c
+++ b/apps/plugins/lib/fixedpoint.c
@@ -1,238 +0,0 @@
1/***************************************************************************
2 * __________ __ ___.
3 * Open \______ \ ____ ____ | | _\_ |__ _______ ___
4 * Source | _// _ \_/ ___\| |/ /| __ \ / _ \ \/ /
5 * Jukebox | | ( <_> ) \___| < | \_\ ( <_> > < <
6 * Firmware |____|_ /\____/ \___ >__|_ \|___ /\____/__/\_ \
7 * \/ \/ \/ \/ \/
8 * $Id$
9 *
10 * Copyright (C) 2006 Jens Arnold
11 *
12 * Fixed point library for plugins
13 *
14 * This program is free software; you can redistribute it and/or
15 * modify it under the terms of the GNU General Public License
16 * as published by the Free Software Foundation; either version 2
17 * of the License, or (at your option) any later version.
18 *
19 * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY OF ANY
20 * KIND, either express or implied.
21 *
22 ****************************************************************************/
23
24#include <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/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}