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__) 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 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 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 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 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 184 static inline float warp_asinf(float x) 185 { 186 return asinf(x) * TWO_OVER_PI_FLOAT; 187 } 188 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 */ 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
