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