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