1 /* Copyright (c) 2013 The Chromium OS Authors. All rights reserved.
2 * Use of this source code is governed by a BSD-style license that can be
3 * found in the LICENSE file.
4 */
5
6 #ifndef DRC_MATH_H_
7 #define DRC_MATH_H_
8
9 #ifdef __cplusplus
10 extern "C" {
11 #endif
12
13 #include <stddef.h>
14 #include <math.h>
15
16 /* Uncomment to use the slow but accurate functions. */
17 /* #define SLOW_DB_TO_LINEAR */
18 /* #define SLOW_LINEAR_TO_DB */
19 /* #define SLOW_WARP_SIN */
20 /* #define SLOW_KNEE_EXP */
21 /* #define SLOW_FREXPF */
22
23 #define PI_FLOAT 3.141592653589793f
24 #define PI_OVER_TWO_FLOAT 1.57079632679489661923f
25 #define TWO_OVER_PI_FLOAT 0.63661977236758134f
26 #define NEG_TWO_DB 0.7943282347242815f /* -2dB = 10^(-2/20) */
27
28 #ifndef max
29 #define max(a, b) ({ __typeof__(a) _a = (a); \
30 __typeof__(b) _b = (b); \
31 _a > _b ? _a : _b; })
32 #endif
33
34 #ifndef min
35 #define min(a, b) ({ __typeof__(a) _a = (a); \
36 __typeof__(b) _b = (b); \
37 _a < _b ? _a : _b; })
38 #endif
39
40 #ifndef M_PI
41 #define M_PI 3.14159265358979323846
42 #endif
43
44 #define PURE __attribute__ ((pure))
45 static inline float decibels_to_linear(float decibels) PURE;
46 static inline float linear_to_decibels(float linear) PURE;
47 static inline float warp_sinf(float x) PURE;
48 static inline float warp_asinf(float x) PURE;
49 static inline float knee_expf(float input) PURE;
50
51 extern float db_to_linear[201]; /* from -100dB to 100dB */
52
53 void drc_math_init();
54
55 union ieee754_float {
56 float f;
57
58 /* This is the IEEE 754 single-precision format. */
59 struct
60 {
61 /* Little endian. */
62 unsigned int mantissa:23;
63 unsigned int exponent:8;
64 unsigned int negative:1;
65 } ieee;
66 };
67
68 /* Rounds the input number to the nearest integer */
69 #ifdef __arm__
round_int(float x)70 static inline float round_int(float x)
71 {
72 return x < 0 ? (int)(x - 0.5f) : (int)(x + 0.5f);
73 }
74 #else
75 #define round_int rintf /* glibc will use roundss if SSE4.1 is available */
76 #endif
77
decibels_to_linear(float decibels)78 static inline float decibels_to_linear(float decibels)
79 {
80 #ifdef SLOW_DB_TO_LINEAR
81 /* 10^(x/20) = e^(x * log(10^(1/20))) */
82 return expf(0.1151292546497022f * decibels);
83 #else
84 float x;
85 float fi;
86 int i;
87
88 fi = round_int(decibels);
89 x = decibels - fi;
90 i = (int)fi;
91 i = max(min(i, 100), -100);
92
93 /* Coefficients obtained from:
94 * fpminimax(10^(x/20), [|1,2,3|], [|SG...|], [-0.5;0.5], 1, absolute);
95 * max error ~= 7.897e-8
96 */
97 const float A3 = 2.54408805631101131439208984375e-4f;
98 const float A2 = 6.628888659179210662841796875e-3f;
99 const float A1 = 0.11512924730777740478515625f;
100 const float A0 = 1.0f;
101
102 float x2 = x * x;
103 return ((A3 * x + A2)*x2 + (A1 * x + A0)) * db_to_linear[i+100];
104 #endif
105 }
106
frexpf_fast(float x,int * e)107 static inline float frexpf_fast(float x, int *e)
108 {
109 #ifdef SLOW_FREXPF
110 return frexpf(x, e);
111 #else
112 union ieee754_float u;
113 u.f = x;
114 int exp = u.ieee.exponent;
115 if (exp == 0xff)
116 return NAN;
117 *e = exp - 126;
118 u.ieee.exponent = 126;
119 return u.f;
120 #endif
121 }
122
linear_to_decibels(float linear)123 static inline float linear_to_decibels(float linear)
124 {
125 /* For negative or zero, just return a very small dB value. */
126 if (linear <= 0)
127 return -1000;
128
129 #ifdef SLOW_LINEAR_TO_DB
130 /* 20 * log10(x) = 20 / log(10) * log(x) */
131 return 8.6858896380650366f * logf(linear);
132 #else
133 int e = 0;
134 float x = frexpf_fast(linear, &e);
135 float exp = e;
136
137 if (x > 0.707106781186548f) {
138 x *= 0.707106781186548f;
139 exp += 0.5f;
140 }
141
142 /* Coefficients obtained from:
143 * fpminimax(log10(x), 5, [|SG...|], [1/2;sqrt(2)/2], absolute);
144 * max err ~= 6.088e-8
145 */
146 const float A5 = 1.131880283355712890625f;
147 const float A4 = -4.258677959442138671875f;
148 const float A3 = 6.81631565093994140625f;
149 const float A2 = -6.1185703277587890625f;
150 const float A1 = 3.6505267620086669921875f;
151 const float A0 = -1.217894077301025390625f;
152
153 float x2 = x * x;
154 float x4 = x2 * x2;
155 return ((A5 * x + A4)*x4 + (A3 * x + A2)*x2 + (A1 * x + A0)) * 20.0f
156 + exp * 6.0205999132796239f;
157 #endif
158 }
159
160
warp_sinf(float x)161 static inline float warp_sinf(float x)
162 {
163 #ifdef SLOW_WARP_SIN
164 return sinf(PI_OVER_TWO_FLOAT * x);
165 #else
166 /* Coefficients obtained from:
167 * fpminimax(sin(x*pi/2), [|1,3,5,7|], [|SG...|], [-1e-30;1], absolute)
168 * max err ~= 5.901e-7
169 */
170 const float A7 = -4.3330336920917034149169921875e-3f;
171 const float A5 = 7.9434238374233245849609375e-2f;
172 const float A3 = -0.645892798900604248046875f;
173 const float A1 = 1.5707910060882568359375f;
174
175 float x2 = x * x;
176 float x4 = x2 * x2;
177 return x * ((A7 * x2 + A5) * x4 + (A3 * x2 + A1));
178 #endif
179 }
180
warp_asinf(float x)181 static inline float warp_asinf(float x)
182 {
183 return asinf(x) * TWO_OVER_PI_FLOAT;
184 }
185
knee_expf(float input)186 static inline float knee_expf(float input)
187 {
188 #ifdef SLOW_KNEE_EXP
189 return expf(input);
190 #else
191 /* exp(x) = decibels_to_linear(20*log10(e)*x) */
192 return decibels_to_linear(8.685889638065044f * input);
193 #endif
194 }
195
196 /* Returns 1 for nan or inf, 0 otherwise. This is faster than the alternative
197 * return x != 0 && !isnormal(x);
198 */
isbadf(float x)199 static inline int isbadf(float x)
200 {
201 union ieee754_float u;
202 u.f = x;
203 return u.ieee.exponent == 0xff;
204 }
205
206 #ifdef __cplusplus
207 } /* extern "C" */
208 #endif
209
210 #endif /* DRC_MATH_H_ */
211