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