special_functions.hpp

Go to the documentation of this file.

#pragma once
 
// Special mathematical functions for DSP and filter design.
//
// RT contract: the scalar free functions and nested `special::*` helpers are
// stateless and allocate no memory. They may be computationally expensive and
// are often better suited to prepare/design paths, but allocation behavior is
// fixed.
 
#include <algorithm>
#include <cmath>
#include <limits>
 
namespace pulp::signal {

inline float sinc(float x) {
    if (std::abs(x) < 1e-7f) return 1.0f;
    float px = 3.14159265358979f * x;
    return std::sin(px) / px;
}

inline float bessel_i0(float x) {
    float sum = 1.0f;
    float term = 1.0f;
    float x_half = x * 0.5f;
 
    for (int k = 1; k < 25; ++k) {
        term *= (x_half / static_cast<float>(k));
        term *= (x_half / static_cast<float>(k));
        sum += term;
        if (term < sum * 1e-10f) break;
    }
    return sum;
}

inline float gamma_fn(float x) {
    return std::tgamma(x);
}

inline float erf_fn(float x) {
    return std::erf(x);
}

inline float erfc_fn(float x) {
    return std::erfc(x);
}

inline float lanczos(float x, int a) {
    if (std::abs(x) < 1e-7f) return 1.0f;
    if (std::abs(x) >= static_cast<float>(a)) return 0.0f;
    return sinc(x) * sinc(x / static_cast<float>(a));
}

inline float db_to_linear(float db) {
    return std::pow(10.0f, db / 20.0f);
}

inline float linear_to_db(float linear) {
    return 20.0f * std::log10(std::max(linear, 1e-10f));
}

inline float freq_to_midi(float freq) {
    return 69.0f + 12.0f * std::log2(freq / 440.0f);
}

inline float midi_to_freq(float note) {
    return 440.0f * std::pow(2.0f, (note - 69.0f) / 12.0f);
}
 
// ---------------------------------------------------------------------------
// Elliptic / Jacobi special functions
//
// These primitives are needed by elliptic (Cauer) IIR filter design and by
// future analog-prototype filter routines. They live in a nested namespace
// so the surface is self-contained and easy to reference from filter design
// code (e.g. `pulp::signal::special::elliptic_K(m)`), without colliding with
// the existing free functions above.
//
// All routines are header-only and `double`-precision because filter design
// requires more headroom than the audio path: coefficient quantization is
// the bottleneck, so we keep design-time math in double and downcast at
// the biquad-cascade boundary.
//
// References:
//   * Abramowitz & Stegun, "Handbook of Mathematical Functions", §17.
//   * Numerical Recipes in C, §6.11 (elliptic integrals and functions).
//   * Vaidyanathan, "Multirate Systems and Filter Banks", §3.
//
// All functions take the *parameter* m = k^2 (Abramowitz convention), not
// the modulus k. Callers used to SciPy's `ellipk(m)` or MATLAB's
// `ellipke(m)` will find this familiar.
// ---------------------------------------------------------------------------
 
namespace special {

inline constexpr double kPi = 3.141592653589793238462643383279502884;

inline double elliptic_K(double m) {
    if (m < 0.0 || m >= 1.0) {
        if (m == 1.0) return std::numeric_limits<double>::infinity();
        // Caller passed something out of [0, 1); clamp into the valid range
        // rather than NaN-poisoning downstream filter design code. For m < 0
        // there are extension formulas, but we currently have no caller that
        // needs them.
        if (m < 0.0) m = 0.0;
        else m = std::nextafter(1.0, 0.0);
    }
    double a = 1.0;
    double b = std::sqrt(1.0 - m);
    for (int i = 0; i < 60; ++i) {
        const double a_next = 0.5 * (a + b);
        const double b_next = std::sqrt(a * b);
        if (std::abs(a - b) < 1e-16 * std::abs(a)) {
            a = a_next;
            break;
        }
        a = a_next;
        b = b_next;
    }
    return kPi / (2.0 * a);
}

inline double elliptic_E(double m) {
    if (m <= 0.0) return kPi / 2.0;
    if (m >= 1.0) return 1.0;
    double a = 1.0;
    double b = std::sqrt(1.0 - m);
    double c = std::sqrt(m);              // c_0
    double pow2 = 0.5;                    // 2^{n-1} for n=0
    double sum = pow2 * c * c;            // first term: m/2
    for (int i = 1; i < 60; ++i) {
        const double a_next = 0.5 * (a + b);
        const double b_next = std::sqrt(a * b);
        c = 0.5 * (a - b);                // c_i
        a = a_next;
        b = b_next;
        pow2 *= 2.0;                      // now 2^{i-1}
        sum += pow2 * c * c;
        if (std::abs(c) < 1e-17 * std::abs(a)) break;
    }
    const double K = kPi / (2.0 * a);
    return K * (1.0 - sum);
}

inline double jacobi_nome(double m) {
    if (m <= 0.0) return 0.0;
    if (m >= 1.0) return 1.0;
    const double K  = elliptic_K(m);
    const double Kp = elliptic_K(1.0 - m);
    return std::exp(-kPi * Kp / K);
}

inline void jacobi_sncndn(double u, double m,
                          double* sn_out, double* cn_out, double* dn_out) {
    if (m <= 0.0) {
        if (sn_out) *sn_out = std::sin(u);
        if (cn_out) *cn_out = std::cos(u);
        if (dn_out) *dn_out = 1.0;
        return;
    }
    if (m >= 1.0) {
        const double th = std::tanh(u);
        const double sech = 1.0 / std::cosh(u);
        if (sn_out) *sn_out = th;
        if (cn_out) *cn_out = sech;
        if (dn_out) *dn_out = sech;
        return;
    }
 
    // Numerical Recipes 3rd ed. §6.12 `sncndn`. The variable names stay
    // close to the published code so it's easy to cross-check, but Pulp
    // uses descriptive comments instead of one-letter mystery.
    constexpr int kMaxN = 14;
    constexpr double kEps = 1.0e-15;
    double em[kMaxN];                 // sequence of a values (AGM upper)
    double en[kMaxN];                 // sequence of complementary moduli
    double emc = 1.0 - m;             // complementary parameter (1 - m)
    double a   = 1.0;
    double dn  = 1.0;
    double c   = 0.0;
    int last_i = 0;
    for (int i = 0; i < kMaxN; ++i) {
        last_i = i;
        em[i] = a;
        emc   = std::sqrt(emc);
        en[i] = emc;
        c     = 0.5 * (a + emc);
        if (std::abs(a - emc) <= kEps * a) break;
        emc = emc * a;                // emc <- a * emc
        a   = c;
    }
    u *= c;                           // u <- final c * u (not a * u)
    double sn = std::sin(u);
    double cn = std::cos(u);
    if (sn != 0.0) {
        a = cn / sn;
        c *= a;
        for (int i = last_i; i >= 0; --i) {
            const double b = em[i];
            a *= c;
            c *= dn;
            dn = (en[i] + a) / (b + a);
            a  = c / b;
        }
        a  = 1.0 / std::sqrt(c * c + 1.0);
        sn = (sn >= 0.0) ? a : -a;
        cn = sn * c;
    }
    if (sn_out) *sn_out = sn;
    if (cn_out) *cn_out = cn;
    if (dn_out) *dn_out = dn;
}

inline double jacobi_sn(double u, double m) {
    double sn = 0.0;
    jacobi_sncndn(u, m, &sn, nullptr, nullptr);
    return sn;
}

inline double jacobi_cn(double u, double m) {
    double cn = 0.0;
    jacobi_sncndn(u, m, nullptr, &cn, nullptr);
    return cn;
}

inline double jacobi_dn(double u, double m) {
    double dn = 0.0;
    jacobi_sncndn(u, m, nullptr, nullptr, &dn);
    return dn;
}

inline double carlson_RF(double x, double y, double z) {
    if (x < 0.0) x = 0.0;
    if (y < 0.0) y = 0.0;
    if (z < 0.0) z = 0.0;
    constexpr double kErrTol = 0.0025;  // ~ machine eps^(1/6); NR recommendation.
    for (int i = 0; i < 60; ++i) {
        const double sqrt_x = std::sqrt(x);
        const double sqrt_y = std::sqrt(y);
        const double sqrt_z = std::sqrt(z);
        const double lambda = sqrt_x * sqrt_y + sqrt_y * sqrt_z + sqrt_z * sqrt_x;
        x = 0.25 * (x + lambda);
        y = 0.25 * (y + lambda);
        z = 0.25 * (z + lambda);
        const double avg = (x + y + z) / 3.0;
        const double dx = (avg - x) / avg;
        const double dy = (avg - y) / avg;
        const double dz = (avg - z) / avg;
        const double max_d = std::max({std::abs(dx), std::abs(dy), std::abs(dz)});
        if (max_d < kErrTol) {
            const double e2 = dx * dy - dz * dz;
            const double e3 = dx * dy * dz;
            // 7th-order series expansion (NR eq. 6.12.17):
            return (1.0 + e2 * (-0.1 + e2 / 24.0 - 3.0 * e3 / 44.0)
                        + e3 / 14.0) / std::sqrt(avg);
        }
    }
    return 1.0 / std::sqrt((x + y + z) / 3.0);  // fallback
}

inline double jacobi_asn(double x, double m) {
    if (x >  1.0) x =  1.0;
    if (x < -1.0) x = -1.0;
    if (m <= 0.0) return std::asin(x);
    if (m >= 1.0) return std::atanh(x);
    const double y = 1.0 - x * x;
    const double z = 1.0 - m * x * x;
    return x * carlson_RF(y, z, 1.0);
}
 
}  // namespace special
 
}  // namespace pulp::signal