10f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid/* Copyright (c) 2013 The Chromium OS Authors. All rights reserved. 20f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid * Use of this source code is governed by a BSD-style license that can be 30f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid * found in the LICENSE file. 40f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid */ 50f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid 60f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid#ifndef DRC_MATH_H_ 70f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid#define DRC_MATH_H_ 80f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid 90f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid#ifdef __cplusplus 100f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reidextern "C" { 110f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid#endif 120f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid 130f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid#include <stddef.h> 140f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid#include <math.h> 150f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid 160f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid/* Uncomment to use the slow but accurate functions. */ 170f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid/* #define SLOW_DB_TO_LINEAR */ 180f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid/* #define SLOW_LINEAR_TO_DB */ 190f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid/* #define SLOW_WARP_SIN */ 200f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid/* #define SLOW_KNEE_EXP */ 210f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid/* #define SLOW_FREXPF */ 220f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid 230f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid#define PI_FLOAT 3.141592653589793f 240f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid#define PI_OVER_TWO_FLOAT 1.57079632679489661923f 250f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid#define TWO_OVER_PI_FLOAT 0.63661977236758134f 260f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid#define NEG_TWO_DB 0.7943282347242815f /* -2dB = 10^(-2/20) */ 270f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid 280f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid#ifndef max 290f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid#define max(a, b) ({ __typeof__(a) _a = (a); \ 300f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid __typeof__(b) _b = (b); \ 310f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid _a > _b ? _a : _b; }) 320f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid#endif 330f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid 340f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid#ifndef min 350f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid#define min(a, b) ({ __typeof__(a) _a = (a); \ 360f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid __typeof__(b) _b = (b); \ 370f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid _a < _b ? _a : _b; }) 380f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid#endif 390f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid 400f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid#ifndef M_PI 410f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid#define M_PI 3.14159265358979323846 420f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid#endif 430f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid 440f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid#define PURE __attribute__ ((pure)) 450f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reidstatic inline float decibels_to_linear(float decibels) PURE; 460f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reidstatic inline float linear_to_decibels(float linear) PURE; 470f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reidstatic inline float warp_sinf(float x) PURE; 480f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reidstatic inline float warp_asinf(float x) PURE; 490f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reidstatic inline float knee_expf(float input) PURE; 500f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid 510f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reidextern float db_to_linear[201]; /* from -100dB to 100dB */ 520f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid 530f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reidvoid drc_math_init(); 540f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid 550f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reidunion ieee754_float { 560f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid float f; 570f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid 580f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid /* This is the IEEE 754 single-precision format. */ 590f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid struct 600f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid { 610f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid /* Little endian. */ 620f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid unsigned int mantissa:23; 630f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid unsigned int exponent:8; 640f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid unsigned int negative:1; 650f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid } ieee; 660f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid}; 670f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid 680f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid/* Rounds the input number to the nearest integer */ 690f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid#ifdef __arm__ 700f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reidstatic inline float round_int(float x) 710f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid{ 720f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid return x < 0 ? (int)(x - 0.5f) : (int)(x + 0.5f); 730f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid} 740f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid#else 750f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid#define round_int rintf /* glibc will use roundss if SSE4.1 is available */ 760f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid#endif 770f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid 780f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reidstatic inline float decibels_to_linear(float decibels) 790f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid{ 800f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid#ifdef SLOW_DB_TO_LINEAR 810f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid /* 10^(x/20) = e^(x * log(10^(1/20))) */ 820f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid return expf(0.1151292546497022f * decibels); 830f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid#else 840f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid float x; 850f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid float fi; 860f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid int i; 870f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid 880f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid fi = round_int(decibels); 890f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid x = decibels - fi; 900f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid i = (int)fi; 910f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid i = max(min(i, 100), -100); 920f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid 930f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid /* Coefficients obtained from: 940f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid * fpminimax(10^(x/20), [|1,2,3|], [|SG...|], [-0.5;0.5], 1, absolute); 950f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid * max error ~= 7.897e-8 960f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid */ 970f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid const float A3 = 2.54408805631101131439208984375e-4f; 980f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid const float A2 = 6.628888659179210662841796875e-3f; 990f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid const float A1 = 0.11512924730777740478515625f; 1000f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid const float A0 = 1.0f; 1010f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid 1020f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid float x2 = x * x; 1030f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid return ((A3 * x + A2)*x2 + (A1 * x + A0)) * db_to_linear[i+100]; 1040f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid#endif 1050f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid} 1060f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid 1070f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reidstatic inline float frexpf_fast(float x, int *e) 1080f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid{ 1090f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid#ifdef SLOW_FREXPF 1100f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid return frexpf(x, e); 1110f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid#else 1120f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid union ieee754_float u; 1130f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid u.f = x; 1140f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid int exp = u.ieee.exponent; 1150f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid if (exp == 0xff) 1160f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid return NAN; 1170f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid *e = exp - 126; 1180f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid u.ieee.exponent = 126; 1190f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid return u.f; 1200f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid#endif 1210f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid} 1220f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid 1230f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reidstatic inline float linear_to_decibels(float linear) 1240f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid{ 1250f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid /* For negative or zero, just return a very small dB value. */ 1260f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid if (linear <= 0) 1270f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid return -1000; 1280f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid 1290f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid#ifdef SLOW_LINEAR_TO_DB 1300f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid /* 20 * log10(x) = 20 / log(10) * log(x) */ 1310f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid return 8.6858896380650366f * logf(linear); 1320f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid#else 1330f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid int e = 0; 1340f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid float x = frexpf_fast(linear, &e); 1350f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid float exp = e; 1360f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid 1370f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid if (x > 0.707106781186548f) { 1380f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid x *= 0.707106781186548f; 1390f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid exp += 0.5f; 1400f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid } 1410f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid 1420f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid /* Coefficients obtained from: 1430f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid * fpminimax(log10(x), 5, [|SG...|], [1/2;sqrt(2)/2], absolute); 1440f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid * max err ~= 6.088e-8 1450f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid */ 1460f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid const float A5 = 1.131880283355712890625f; 1470f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid const float A4 = -4.258677959442138671875f; 1480f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid const float A3 = 6.81631565093994140625f; 1490f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid const float A2 = -6.1185703277587890625f; 1500f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid const float A1 = 3.6505267620086669921875f; 1510f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid const float A0 = -1.217894077301025390625f; 1520f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid 1530f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid float x2 = x * x; 1540f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid float x4 = x2 * x2; 1550f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid return ((A5 * x + A4)*x4 + (A3 * x + A2)*x2 + (A1 * x + A0)) * 20.0f 1560f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid + exp * 6.0205999132796239f; 1570f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid#endif 1580f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid} 1590f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid 1600f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid 1610f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reidstatic inline float warp_sinf(float x) 1620f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid{ 1630f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid#ifdef SLOW_WARP_SIN 1640f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid return sinf(PI_OVER_TWO_FLOAT * x); 1650f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid#else 1660f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid /* Coefficients obtained from: 1670f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid * fpminimax(sin(x*pi/2), [|1,3,5,7|], [|SG...|], [-1e-30;1], absolute) 1680f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid * max err ~= 5.901e-7 1690f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid */ 1700f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid const float A7 = -4.3330336920917034149169921875e-3f; 1710f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid const float A5 = 7.9434238374233245849609375e-2f; 1720f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid const float A3 = -0.645892798900604248046875f; 1730f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid const float A1 = 1.5707910060882568359375f; 1740f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid 1750f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid float x2 = x * x; 1760f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid float x4 = x2 * x2; 1770f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid return x * ((A7 * x2 + A5) * x4 + (A3 * x2 + A1)); 1780f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid#endif 1790f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid} 1800f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid 1810f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reidstatic inline float warp_asinf(float x) 1820f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid{ 1830f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid return asinf(x) * TWO_OVER_PI_FLOAT; 1840f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid} 1850f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid 1860f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reidstatic inline float knee_expf(float input) 1870f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid{ 1880f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid#ifdef SLOW_KNEE_EXP 1890f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid return expf(input); 1900f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid#else 1910f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid /* exp(x) = decibels_to_linear(20*log10(e)*x) */ 1920f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid return decibels_to_linear(8.685889638065044f * input); 1930f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid#endif 1940f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid} 1950f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid 1960f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid/* Returns 1 for nan or inf, 0 otherwise. This is faster than the alternative 1970f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid * return x != 0 && !isnormal(x); 1980f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid */ 1990f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reidstatic inline int isbadf(float x) 2000f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid{ 2010f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid union ieee754_float u; 2020f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid u.f = x; 2030f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid return u.ieee.exponent == 0xff; 2040f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid} 2050f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid 2060f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid#ifdef __cplusplus 2070f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid} /* extern "C" */ 2080f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid#endif 2090f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid 2100f0250476fcb27f45d4f9f98e5769951a5ea5712Dylan Reid#endif /* DRC_MATH_H_ */ 211