MinBlepTable.h

// DSPark -- Professional Audio DSP Framework
// Copyright (c) 2026 Cristian Moresi -- MIT License
 
#pragma once
 
#include "FFT.h"
#include "WindowFunctions.h"
 
#include <array>
#include <cmath>
#include <cstddef>
#include <numbers>
#include <type_traits>
#include <vector>
 
namespace dspark {
 
template <typename T>
class MinBlepTable
{
    static_assert(std::is_floating_point_v<T>, "MinBlepTable requires float or double");
 
public:
    static constexpr int kTaps = 64;
    static constexpr int kOversample = 64;
    static constexpr int kTableSize = kTaps * kOversample + 1;
 
    static const MinBlepTable& instance() noexcept
    {
        static const MinBlepTable table;
        return table;
    }
 
    [[nodiscard]] T residual(T t) const noexcept
    {
        const T pos = t * static_cast<T>(kOversample);
        const auto idx = static_cast<int>(pos);
        if (idx < 0 || idx >= kTableSize - 1)
            return T(0);
        const T frac = pos - static_cast<T>(idx);
        const T a = residual_[static_cast<size_t>(idx)];
        const T b = residual_[static_cast<size_t>(idx) + 1];
        return a + frac * (b - a);
    }
 
    [[nodiscard]] T dcDelay() const noexcept { return dcDelay_; }
 
private:
    MinBlepTable() noexcept { build(); }
 
    void build() noexcept
    {
        constexpr int    n       = kTableSize;     // 4097 oversampled points
        constexpr size_t fftSize = 65536;          // ~16x n keeps cepstral aliasing low
        constexpr double pi      = std::numbers::pi_v<double>;
 
        // -- 1. Linear-phase BLIT: sinc cutting at the base-rate Nyquist
        //       (zeros every kOversample points), Blackman-Harris windowed
        //       (-92 dB stopband), normalised to unit DC gain so the
        //       integrated step settles at 1.
        std::vector<double> x(fftSize, 0.0);
        std::vector<double> win(n);
        WindowFunctions<double>::blackmanHarris(win.data(), n, false);
        const double center = 0.5 * (n - 1);
        double dc = 0.0;
        for (int i = 0; i < n; ++i)
        {
            const double t = (static_cast<double>(i) - center) / kOversample;
            const double s = (std::abs(t) < 1e-9) ? 1.0 : std::sin(pi * t) / (pi * t);
            x[static_cast<size_t>(i)] = s * win[static_cast<size_t>(i)];
            dc += s * win[static_cast<size_t>(i)];
        }
        for (int i = 0; i < n; ++i)
            x[static_cast<size_t>(i)] /= dc;
 
        // -- 2. Real cepstrum. log|X| is a real, even function of frequency,
        //       so running it through the real inverse FFT (imaginary parts
        //       zero) yields the even cepstrum directly.
        FFTReal<double> fft(fftSize);
        std::vector<double> spec(fftSize + 2);
        std::vector<double> cep(fftSize);
        fft.forward(x.data(), spec.data());
 
        const size_t bins = fftSize / 2 + 1;
        for (size_t k = 0; k < bins; ++k)
        {
            const double re  = spec[2 * k];
            const double im  = spec[2 * k + 1];
            // Soft floor at -100 dB, blended in power. The floor must sit just
            // below the window's -92 dB stopband, not far below it: hard, deep
            // corners in log|X| make the cepstrum decay like 1/q^2 and alias
            // through the causal fold, which comes back as percent-level error
            // in the reconstructed step (measured, not hypothetical).
            constexpr double floorMag = 1e-5;
            const double mag = std::sqrt(re * re + im * im + floorMag * floorMag);
            spec[2 * k]     = std::log(mag);
            spec[2 * k + 1] = 0.0;
        }
        fft.inverse(spec.data(), cep.data());
 
        // -- 3. Fold onto the causal part (Hilbert relation between log
        //       magnitude and minimum phase): keep q=0 and q=N/2, double
        //       1..N/2-1, zero the anticausal half.
        for (size_t q = 1; q < fftSize / 2; ++q)
        {
            cep[q]           *= 2.0;
            cep[fftSize - q]  = 0.0;
        }
 
        // -- 4. Back to the spectrum, exponentiate, back to time: the
        //       minimum-phase BLIT, energy packed at the front.
        fft.forward(cep.data(), spec.data());
        for (size_t k = 0; k < bins; ++k)
        {
            const double m  = std::exp(spec[2 * k]);
            const double ph = spec[2 * k + 1];
            spec[2 * k]     = m * std::cos(ph);
            spec[2 * k + 1] = m * std::sin(ph);
        }
        fft.inverse(spec.data(), x.data());
 
        // -- 5. Integrate into a step and store the residual. Pinning the
        //       table end exactly at step==1 closes the truncated residual
        //       at zero (no leftover micro-step when the kernel expires).
        double acc = 0.0;
        std::vector<double> step(n);
        for (int i = 0; i < n; ++i)
        {
            acc += x[static_cast<size_t>(i)];
            step[static_cast<size_t>(i)] = acc;
        }
        const double settle = step[static_cast<size_t>(n - 1)];
        double area = 0.0;
        for (int i = 0; i < n; ++i)
        {
            const double r = step[static_cast<size_t>(i)] / settle - 1.0;
            residual_[static_cast<size_t>(i)] = static_cast<T>(r);
            // Trapezoidal integral of the stored residual (what residual()
            // interpolates), in samples.
            area += (i == 0 || i == n - 1) ? 0.5 * r : r;
        }
        dcDelay_ = static_cast<T>(-area / kOversample);
    }
 
    std::array<T, kTableSize> residual_{};
    T dcDelay_ = T(0);
};
 
} // namespace dspark