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
10extern "C" {
11#endif
12
13#include <stddef.h>
14#include <math.h>
15
16#ifndef __BIONIC__
17#include <ieee754.h>
18#else
19union 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))
61static inline float decibels_to_linear(float decibels) PURE;
62static inline float linear_to_decibels(float linear) PURE;
63static inline float warp_sinf(float x) PURE;
64static inline float warp_asinf(float x) PURE;
65static inline float knee_expf(float input) PURE;
66
67extern float db_to_linear[201]; /* from -100dB to 100dB */
68
69void drc_math_init();
70
71/* Rounds the input number to the nearest integer */
72#if defined(__arm__)
73static 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
81static 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
110static 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
126static 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
164static 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
184static inline float warp_asinf(float x)
185{
186	return asinf(x) * TWO_OVER_PI_FLOAT;
187}
188
189static 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 */
202static 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