Hilbert.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 "DenormalGuard.h"
#include "SimdOps.h" // SIMD dot product for the FIR convolution
 
#include <algorithm>
#include <array>
#include <cassert>
#include <cmath>
#include <span>
 
namespace dspark {
 
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 = ((n == 0) || (n % 2 == 0)) ? 0.0 : (2.0 / (kPi * static_cast<double>(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_{};
};
 
} // namespace dspark