Interpolation.h

Go to the documentation of this file.

// DSPark - Professional Audio DSP Framework
// Copyright (c) 2026 Cristian Moresi - MIT License
 
#pragma once
 
#include "DspMath.h"
 
#include <cassert>
#include <cmath>
 
namespace dspark {
 
template <FloatType T>
[[nodiscard]] inline T interpolateLinear(T y0, T y1, T frac) noexcept {
    return y0 + frac * (y1 - y0);
}
 
template <FloatType T>
[[nodiscard]] inline T interpolateLinear(const T* buffer, int length, T position) noexcept {
    assert(length > 0 && position >= T(0) && position < static_cast<T>(length));
    if (length <= 0) return T(0);
 
    // An out-of-range position (including NaN) would make the int cast below
    // undefined and the buffer reads out of bounds: clamp before casting.
    int idx0;
    T frac;
    if (!(position >= T(0))) { idx0 = 0; frac = T(0); }
    else if (position >= static_cast<T>(length)) { idx0 = length - 1; frac = T(0); }
    else {
        idx0 = static_cast<int>(position);
        frac = position - static_cast<T>(idx0);
    }
 
    int idx1 = idx0 + 1;
    if (idx1 >= length) idx1 = 0;
 
    return interpolateLinear(buffer[idx0], buffer[idx1], frac);
}
 
template <FloatType T>
[[nodiscard]] inline T interpolateHermite(T y0, T y1, T y2, T y3, T frac) noexcept {
    T c = (y2 - y0) * T(0.5);
    T v = y1 - y2;
    T w = c + v;
    T a = w + v + (y3 - y1) * T(0.5);
    T b = w + a;
    return (((a * frac) - b) * frac + c) * frac + y1;
}
 
template <FloatType T>
[[nodiscard]] inline T interpolateHermite(const T* buffer, int length, T position) noexcept {
    assert(length > 0 && position >= T(0) && position < static_cast<T>(length));
    if (length <= 0) return T(0);
 
    // Same release-safe clamp as the linear overload (undefined int cast /
    // out-of-bounds reads otherwise).
    int idx1;
    T frac;
    if (!(position >= T(0))) { idx1 = 0; frac = T(0); }
    else if (position >= static_cast<T>(length)) { idx1 = length - 1; frac = T(0); }
    else {
        idx1 = static_cast<int>(position);
        frac = position - static_cast<T>(idx1);
    }
 
    int idx0 = (idx1 > 0) ? idx1 - 1 : length - 1;
    int idx2 = idx1 + 1; if (idx2 >= length) idx2 -= length;
    int idx3 = idx2 + 1; if (idx3 >= length) idx3 -= length;
 
    return interpolateHermite(buffer[idx0], buffer[idx1], buffer[idx2], buffer[idx3], frac);
}
 
template <FloatType T>
[[nodiscard]] inline T interpolateCubic(const T* buffer, int length, T position) noexcept {
    return interpolateHermite(buffer, length, position);
}
 
template <FloatType T>
[[nodiscard]] inline T interpolateLagrange(T y0, T y1, T y2, T y3, T frac) noexcept {
    T d = frac;
    T dm1 = d - T(1);
    T dm2 = d - T(2);
    T dp1 = d + T(1);
 
    // Precomputed division multipliers: 1/6 ~ 0.16666667, 1/2 = 0.5
    static constexpr T inv6 = T(1.0 / 6.0);
    static constexpr T inv2 = T(0.5);
 
    T l0 = -(d * dm1 * dm2) * inv6;
    T l1 =  (dp1 * dm1 * dm2) * inv2;
    T l2 = -(dp1 * d * dm2) * inv2;
    T l3 =  (dp1 * d * dm1) * inv6;
 
    return l0 * y0 + l1 * y1 + l2 * y2 + l3 * y3;
}
 
template <FloatType T>
[[nodiscard]] inline T interpolateLagrange(const T* buffer, int length, T position) noexcept {
    assert(length > 0 && position >= T(0) && position < static_cast<T>(length));
    if (length <= 0) return T(0);
 
    // Same release-safe clamp as the linear overload (undefined int cast /
    // out-of-bounds reads otherwise).
    int idx1;
    T frac;
    if (!(position >= T(0))) { idx1 = 0; frac = T(0); }
    else if (position >= static_cast<T>(length)) { idx1 = length - 1; frac = T(0); }
    else {
        idx1 = static_cast<int>(position);
        frac = position - static_cast<T>(idx1);
    }
 
    int idx0 = (idx1 > 0) ? idx1 - 1 : length - 1;
    int idx2 = idx1 + 1; if (idx2 >= length) idx2 -= length;
    int idx3 = idx2 + 1; if (idx3 >= length) idx3 -= length;
 
    return interpolateLagrange(buffer[idx0], buffer[idx1], buffer[idx2], buffer[idx3], frac);
}
 
template <FloatType T>
[[nodiscard]] inline T interpolateAllpass(T currentSample, T previousSample,
                                          T frac, T& state) noexcept
{
    // frac == 0 would place the pole exactly at z = -1 (a marginally stable
    // resonator at Nyquist) and frac == -1 would divide by zero: clamp away
    // from both. The clamp changes the delivered delay only for out-of-range
    // requests.
    T safeFrac = (frac < T(0.001)) ? T(0.001) : frac;
    T coeff = (T(1) - safeFrac) / (T(1) + safeFrac);
 
    T output = coeff * (currentSample - state) + previousSample;
 
    // Cut the recursive tail below ~-300 dBFS so a decaying state can never
    // reach the denormal range (matters on targets without FTZ/DAZ).
    if (std::abs(output) < T(1e-15)) output = T(0);
 
    state = output;
    return output;
}
 
} // namespace dspark