Hilbert.h

// DSPark - Professional Audio DSP Framework
// Copyright (c) 2026 Cristian Moresi - MIT License
 
#pragma once
 
#include "DspMath.h"
#include "DenormalGuard.h"
#include "SimdOps.h" // SIMD dot product for the FIR convolution
 
#include <algorithm>
#include <array>
#include <cassert>
#include <cmath>
#include <cstdint>
#include <span>
 
namespace dspark {
 
namespace detail {
// Unwindowed discrete Hilbert impulse, shared by finite-source operators and
// the real-time windowed FIR. Keep the parity test in the integer source clock.
[[nodiscard]] inline double hilbertIdealImpulse(std::int64_t index) noexcept
{
    constexpr double kPi = 3.14159265358979323846;
    return index % 2 == 0 ? 0.0 : 2.0 / (kPi * static_cast<double>(index));
}
} // namespace detail
 
template <FloatType T>
class Hilbert
{
public:
    struct Result
    {
        T real; 
        T imag; 
    };
 
    void prepare(double sampleRate) noexcept
    {
        assert(sampleRate > 0.0);
        if (!(sampleRate > 0.0)) return;
 
        sampleRate_ = sampleRate;
        reset();
        isPrepared_ = true;
    }
 
    [[nodiscard]] inline Result process(T input) noexcept
    {
        delay_[static_cast<size_t>(writePos_)]         = input;
        delay_[static_cast<size_t>(writePos_ + kTaps)] = input;
 
        // Contiguous window of the last kTaps samples, oldest first:
        // window[j] = x[n - (kTaps - 1 - j)].
        const T* window = delay_.data() + writePos_ + 1;
 
        // imag = sum_k h[k] * x[n-k]: with an oldest-first window this is the
        // dot product against the REVERSED kernel. The pointer is cached in a
        // member so the hot path never touches the magic-static init guard.
        const T imag = simd::dotProduct(kernelData_, window, kTaps);
 
        // real = x[n - kCenter]: index kTaps-1-kCenter == kCenter (odd length).
        const T re = window[kCenter];
 
        writePos_ = (writePos_ + 1 == kTaps) ? 0 : (writePos_ + 1);
        return { re, imag };
    }
 
    void processBlock(std::span<const T> input,
                      std::span<T> outReal,
                      std::span<T> outImag) noexcept
    {
        assert(isPrepared_);
        assert(outReal.size() >= input.size() && outImag.size() >= input.size());
        DenormalGuard dg;
 
        const size_t numSamples =
            std::min(input.size(), std::min(outReal.size(), outImag.size()));
        for (size_t i = 0; i < numSamples; ++i)
        {
            const auto res = process(input[i]);
            outReal[i] = res.real;
            outImag[i] = res.imag;
        }
    }
 
    void reset() noexcept
    {
        delay_.fill(T(0));
        writePos_ = 0;
    }
 
    [[nodiscard]] static constexpr int getLatencySamples() noexcept { return kCenter; }
 
private:
    // 191-tap FIR: < 0.1% analytic ripple above ~1 kHz at 48 kHz (< 2.5% at
    // 500 Hz), 95-sample group delay. Odd length keeps a true integer center tap.
    static constexpr int kTaps   = 191;
    static constexpr int kCenter = kTaps / 2;
 
    [[nodiscard]] static const std::array<T, kTaps>& reversedKernel() noexcept
    {
        static const std::array<T, kTaps> kernel = []
        {
            std::array<T, kTaps> k {};
            constexpr double kPi = 3.14159265358979323846;
            for (int i = 0; i < kTaps; ++i)
            {
                // Ideal Hilbert kernel h[n] = 2/(pi*n) for odd n, 0 for even n,
                // windowed with Blackman to control ripple / sideband leakage.
                const int n = i - kCenter;
                double ideal = detail::hilbertIdealImpulse(n);
                double w = 0.42
                         - 0.5  * std::cos(2.0 * kPi * static_cast<double>(i) / static_cast<double>(kTaps - 1))
                         + 0.08 * std::cos(4.0 * kPi * static_cast<double>(i) / static_cast<double>(kTaps - 1));
                k[static_cast<size_t>(kTaps - 1 - i)] = static_cast<T>(ideal * w);
            }
            return k;
        }();
        return kernel;
    }
 
    double sampleRate_ = 0.0;
    bool isPrepared_ = false;
    int writePos_ = 0;
 
    // Cached at construction so process() skips the magic-static guard check
    // (the static outlives every instance; copies stay valid).
    const T* kernelData_ = reversedKernel().data();
 
    // Mirrored delay line (double-write) for contiguous SIMD reads. No
    // over-alignment: the window starts at a variable offset every sample, so
    // the SIMD dot product uses unaligned loads regardless.
    std::array<T, kTaps * 2> delay_{};
};
 
template <FloatType T>
class HilbertIIR
{
public:
    struct Result
    {
        T real; 
        T imag; 
    };
 
    void reset() noexcept
    {
        for (auto& s : sections_) s = {};
        delayedImag_ = 0.0;
    }
 
    [[nodiscard]] inline Result process(T input) noexcept
    {
        double re = 0.0, im = 0.0;
        step(static_cast<double>(input), re, im);
        return { static_cast<T>(re), static_cast<T>(im) };
    }
 
    [[nodiscard]] inline T magnitude(T input) noexcept
    {
        double re = 0.0, im = 0.0;
        step(static_cast<double>(input), re, im);
        return static_cast<T>(std::sqrt(re * re + im * im));
    }
 
    void magnitudeBlock(const T* in, T* out, int numSamples) noexcept
    {
        constexpr int kChunk = 64;
        double re[kChunk];
        double im[kChunk];
        for (int start = 0; start < numSamples; start += kChunk)
        {
            const int n = std::min(kChunk, numSamples - start);
            for (int i = 0; i < n; ++i)
                re[i] = im[i] = static_cast<double>(in[start + i]);
            for (int k = 0; k < kSections; ++k)
                Section::runPair(sections_[static_cast<size_t>(k)], re, kReal[k],
                                 sections_[static_cast<size_t>(kSections + k)], im, kImag[k], n);
            for (int i = 0; i < n; ++i)
            {
                const double imDelayed = delayedImag_;
                delayedImag_ = im[i];
                out[start + i] = static_cast<T>(std::sqrt(re[i] * re[i] + imDelayed * imDelayed));
            }
        }
    }
 
private:
    static constexpr int kSections = 7;
 
    inline void step(double input, double& re, double& im) noexcept
    {
        re = input;
        double quad = input;
        for (int k = 0; k < kSections; ++k)
        {
            re = sections_[static_cast<size_t>(k)].process(re, kReal[k]);
            quad = sections_[static_cast<size_t>(kSections + k)].process(quad, kImag[k]);
        }
        im = delayedImag_;
        delayedImag_ = quad;
    }
 
    struct Section
    {
        double x1 = 0.0, x2 = 0.0, y1 = 0.0, y2 = 0.0;
 
        // y = (c x[n] - x[n-2]) + c y[n-2]: the recursive term enters last,
        // so the loop-carried path is one multiply and one add.
        inline double process(double x, double c) noexcept
        {
            const double y = (c * x - x2) + c * y2;
            x2 = x1; x1 = x;
            y2 = y1; y1 = y;
            return y;
        }
 
        static inline void runPair(Section& p, double* dp, double cp,
                                   Section& q, double* dq, double cq, int n) noexcept
        {
            double px1 = p.x1, px2 = p.x2, py1 = p.y1, py2 = p.y2;
            double qx1 = q.x1, qx2 = q.x2, qy1 = q.y1, qy2 = q.y2;
            for (int i = 0; i < n; ++i)
            {
                const double xp = dp[i];
                const double xq = dq[i];
                const double yp = (cp * xp - px2) + cp * py2;
                const double yq = (cq * xq - qx2) + cq * qy2;
                px2 = px1; px1 = xp; py2 = py1; py1 = yp;
                qx2 = qx1; qx1 = xq; qy2 = qy1; qy1 = yq;
                dp[i] = yp;
                dq[i] = yq;
            }
            p.x1 = px1; p.x2 = px2; p.y1 = py1; p.y2 = py2;
            q.x1 = qx1; q.x2 = qx2; q.y1 = qy1; q.y2 = qy2;
        }
    };
 
    static constexpr double kReal[kSections] = {
        0.0849750635881296, 0.5135957829646070, 0.8197816632773615, 0.9420612863757530,
        0.9822486171863649, 0.9947076541867014, 0.9986500580725346 };
    static constexpr double kImag[kSections] = {  // followed by the one-sample delay
        0.2887645620083060, 0.6955314021190576, 0.8968265992273647, 0.9678214299445055,
        0.9902622348692460, 0.9971984927047717, 0.9995996402566010 };
 
    std::array<Section, 2 * kSections> sections_ {};
    double delayedImag_ = 0.0;
};
 
} // namespace dspark