DspMath.h

Go to the documentation of this file.

// DSPark -- Professional Audio DSP Framework
// Copyright (c) 2026 Cristian Moresi -- MIT License
 
#pragma once
 
#include <algorithm>
#include <cmath>
#include <concepts>
#include <numbers>
 
namespace dspark {
 
// ============================================================================
// Concepts
// ============================================================================
 
template <typename T>
concept FloatType = std::floating_point<T>;
 
// ============================================================================
// Constants
// ============================================================================
 
template <FloatType T> inline constexpr T pi       = std::numbers::pi_v<T>;
 
template <FloatType T> inline constexpr T twoPi    = T(2) * std::numbers::pi_v<T>;
 
template <FloatType T> inline constexpr T invTwoPi = T(1) / twoPi<T>;
 
template <FloatType T> inline constexpr T halfPi   = std::numbers::pi_v<T> / T(2);
 
template <FloatType T> inline constexpr T sqrt2    = std::numbers::sqrt2_v<T>;
 
template <FloatType T> inline constexpr T invSqrt2 = T(1) / std::numbers::sqrt2_v<T>;
 
// ============================================================================
// Decibel / Gain Conversions
// ============================================================================
 
template <FloatType T>
[[nodiscard]] inline T decibelsToGain(T dB, T minusInfinityDb = T(-100)) noexcept
{
    return dB <= minusInfinityDb ? T(0) : std::pow(T(10), dB / T(20));
}
 
template <FloatType T>
[[nodiscard]] inline T gainToDecibels(T gain, T minusInfinityDb = T(-100)) noexcept
{
    return gain > T(0) ? std::max(minusInfinityDb, T(20) * std::log10(gain))
                       : minusInfinityDb;
}
 
// ============================================================================
// Interpolation and Mapping
// ============================================================================
 
template <FloatType T>
[[nodiscard]] inline T mapRange(T value, T inMin, T inMax, T outMin, T outMax) noexcept
{
    if (inMin == inMax) return outMin; // Safety: prevent NaN/Inf in audio pipeline
    return outMin + (outMax - outMin) * ((value - inMin) / (inMax - inMin));
}
 
// ============================================================================
// Fast Approximations
// ============================================================================
 
template <FloatType T>
[[nodiscard]] inline T fastTanh(T x) noexcept
{
    // std::clamp translates to SIMD min/max intrinsics (no branching).
    // Maintains continuity avoiding step-aliasing at the bounds.
    x = std::clamp(x, T(-3), T(3));
    
    const auto x2 = x * x;
    const auto x4 = x2 * x2;
    // Pade [5,4] approximant: max error < 0.05% for |x| <= 3
    return x * (T(945) + T(105) * x2 + x4) / (T(945) + T(420) * x2 + T(15) * x4);
}
 
template <FloatType T>
[[nodiscard]] inline T fastPow10(T x) noexcept
{
    // log2(10) folded into one constant: rounding it once is more accurate
    // than multiplying log2(e) * ln(10) at runtime, and it saves a multiply
    // (compilers may not reassociate x * a * b without fast-math).
    constexpr T kLog2Of10 =
        static_cast<T>(std::numbers::ln10_v<long double> / std::numbers::ln2_v<long double>);
    return std::exp2(x * kLog2Of10);
}
 
template <FloatType T>
[[nodiscard]] inline T fastExp(T x) noexcept
{
    return std::exp2(x * std::numbers::log2e_v<T>);
}
 
template <FloatType T>
[[nodiscard]] inline T fastTan(T x) noexcept
{
    constexpr T limit = T(1.520);  // ~pi/2 - 0.05
    if (std::abs(x) > limit) return std::tan(x);
    const T x2 = x * x;
    const T x4 = x2 * x2;
    // Pade [5,4]: tan(x) ~= x * (945 - 105*x^2 + x^4) / (945 - 420*x^2 + 15*x^4)
    return x * (T(945) - T(105) * x2 + x4) / (T(945) - T(420) * x2 + T(15) * x4);
}
 
template <FloatType T>
[[nodiscard]] inline T fastSin(T x) noexcept
{
    // Range-reduce to [-pi, pi] with a two-term Cody-Waite split of 2*pi so the
    // subtraction stays accurate in float even for arguments several periods out.
    const T k = std::floor(x * invTwoPi<T> + T(0.5));
    x -= k * T(6.28125);                  // high part (exactly representable)
    x -= k * T(0.0019353071795864769253); // low part of 2*pi
 
    // Fold into [-pi/2, pi/2] where the polynomial converges fast:
    // sin(x) = sin(pi - x) for x > pi/2 (and the odd mirror for x < -pi/2).
    if (x > halfPi<T>)       x = pi<T> - x;
    else if (x < -halfPi<T>) x = -pi<T> - x;
 
    const T x2 = x * x;
    // Minimax coefficients for sin on [-pi/2, pi/2]. Stored to ~7 digits, so
    // the realised max abs error is ~3.2e-6 (coefficient-limited, identical in
    // float and double; measured M-002).
    return x * (T(0.9999999995)
         + x2 * (T(-0.1666666580)
         + x2 * (T(0.0083333075)
         + x2 * (T(-0.0001984090)
         + x2 *  T(0.0000027526)))));
}
 
template <FloatType T>
[[nodiscard]] inline T fastCos(T x) noexcept
{
    return fastSin(x + halfPi<T>);
}
 
template <FloatType T>
[[nodiscard]] inline T fastLog(T x) noexcept
{
    int e = 0;
    T m = std::frexp(x, &e);             // x = m * 2^e, m in [0.5, 1)
    // Normalise mantissa into [sqrt(0.5), sqrt(2)) so the series is centred.
    if (m < T(0.70710678118654752440)) { m *= T(2); --e; }
 
    // Classic atanh form (Cephes): ln(m) = 2*atanh(s), s = (m-1)/(m+1).
    // |s| <= 0.1716 so the odd series converges below 1e-7 with 4 terms.
    const T s  = (m - T(1)) / (m + T(1));
    const T s2 = s * s;
    const T lnm = T(2) * s * (T(1)
                + s2 * (T(1.0 / 3.0)
                + s2 * (T(1.0 / 5.0)
                + s2 *  T(1.0 / 7.0))));
 
    return lnm + static_cast<T>(e) * std::numbers::ln2_v<T>; // + e * ln2
}
 
// ============================================================================
// Utility
// ============================================================================
 
template <FloatType T>
[[nodiscard]] inline T wrapPhase(T phase) noexcept
{
    // floor() correctly handles negative values, preventing branching
    T wrapped = phase - twoPi<T> * std::floor(phase * invTwoPi<T>);
    // Rounding in the k * twoPi product can land wrapped just outside
    // [0, 2*pi) on either side: large phases where the product overshoots
    // the input leave a slightly negative result, and tiny negative phases
    // can round the sum up to exactly twoPi. Enforce the documented
    // half-open range. Order matters: the += branch can itself round up to
    // exactly twoPi, which the second branch then folds to 0.
    if (wrapped < T(0))      wrapped += twoPi<T>;
    if (wrapped >= twoPi<T>) wrapped -= twoPi<T>;
    return wrapped;
}
 
} // namespace dspark