zero_latency_convolver_support.hpp

Go to the documentation of this file.

#pragma once

 
// `convolver_non_uniform.hpp` is named above for the reconciliation, but is
// deliberately NOT included: it pulls the convolver message/slot machinery and
// its CHOC FIFO dependency, and nothing here uses its types.
#include <pulp/signal/denormal.hpp>
#include <pulp/signal/scoped_flush_denormals.hpp>
#include <pulp/signal/units.hpp>
#include <pulp/signal/windowed_sinc_design.hpp>
 
#include <algorithm>
#include <cmath>
#include <complex>
#include <cstddef>
#include <cstdint>
#include <limits>
#include <vector>
 
namespace pulp::signal {

enum class IrNormalizeMode {
    none,    
    peak,    
    energy,  
};
 
namespace detail {

template <typename SampleType>
class SlicedFftPlanT {
public:
    // Cost weights, in real flops per work item. A radix-2 butterfly is one
    // complex multiply (6 flops) plus two complex adds (4 flops); a permuted
    // gather or an in-place swap is one load/store pair; the inverse's
    // conjugate and conjugate-scale passes are a negate, and a negate plus two
    // multiplies, respectively.
    static constexpr long kGatherCost = 1;
    static constexpr long kButterflyCost = 10;
    static constexpr long kConjugateCost = 2;
    static constexpr long kScaleCost = 4;
    static constexpr long kMacCost = 8;
 
    void prepare(int size) {
        size_ = size;
        log2_size_ = 0;
        for (int n = size; n > 1; n >>= 1) ++log2_size_;
        reversed_.resize(static_cast<std::size_t>(size));
        for (int i = 0; i < size; ++i) {
            int r = 0;
            for (int b = 0; b < log2_size_; ++b)
                if (i & (1 << b)) r |= 1 << (log2_size_ - 1 - b);
            reversed_[static_cast<std::size_t>(i)] = r;
        }
        // Compute in double, then store the processing precision once. The
        // butterflies use these same rounded coefficients without converting
        // every table load; float plans also need half the table storage.
        twiddles_.resize(static_cast<std::size_t>(std::max(1, size / 2)));
        constexpr double two_pi = 6.283185307179586476925286766559;
        for (int i = 0; i < size / 2; ++i) {
            const double angle = -two_pi * static_cast<double>(i) / static_cast<double>(size);
            twiddles_[static_cast<std::size_t>(i)] = {static_cast<SampleType>(std::cos(angle)),
                                                      static_cast<SampleType>(std::sin(angle))};
        }
    }
 
    int size() const { return size_; }
    int log2_size() const { return log2_size_; }

    long forward_cost() const {
        return static_cast<long>(size_) * kGatherCost +
               static_cast<long>(log2_size_) * static_cast<long>(size_ / 2) * kButterflyCost;
    }
    long inverse_cost() const {
        return forward_cost() + static_cast<long>(size_) * (kConjugateCost + kScaleCost);
    }

    struct Cursor {
        enum class Phase { idle, gather, butterflies, conjugate, permute, scale };
 
        Phase phase = Phase::idle;
        int cursor = 0;
        int stage = 1;
        bool inverse = false;
        const SampleType* gather_src = nullptr;
        std::size_t gather_base = 0;
        std::size_t gather_mask = 0;
 
        bool done() const { return phase == Phase::idle; }

        void begin_forward_from_ring(const SampleType* ring, std::size_t base, std::size_t mask) {
            phase = Phase::gather;
            cursor = 0;
            stage = 1;
            inverse = false;
            gather_src = ring;
            gather_base = base;
            gather_mask = mask;
        }

        void begin_forward_in_place() {
            phase = Phase::permute;
            cursor = 0;
            stage = 1;
            inverse = false;
            gather_src = nullptr;
        }

        void begin_inverse_in_place() {
            phase = Phase::conjugate;
            cursor = 0;
            stage = 1;
            inverse = true;
            gather_src = nullptr;
        }
 
        void cancel() { phase = Phase::idle; }
 
        long advance(const SlicedFftPlanT& plan, std::complex<SampleType>* data, long budget) {
            long spent = 0;
            const int n = plan.size_;
            while (phase != Phase::idle && spent < budget) {
                switch (phase) {
                case Phase::gather: {
                    const int take = affordable(budget - spent, kGatherCost, n - cursor);
                    if (take <= 0) return spent;
                    for (int i = cursor; i < cursor + take; ++i) {
                        const std::size_t src =
                            (gather_base +
                             static_cast<std::size_t>(plan.reversed_[static_cast<std::size_t>(i)])) &
                            gather_mask;
                        data[i] = {gather_src[src], SampleType{0}};
                    }
                    cursor += take;
                    spent += static_cast<long>(take) * kGatherCost;
                    if (cursor >= n) { phase = Phase::butterflies; cursor = 0; stage = 1; }
                    break;
                }
                case Phase::butterflies: {
                    const int total = n / 2;
                    const int take = affordable(budget - spent, kButterflyCost, total - cursor);
                    if (take <= 0) return spent;
                    const int len = 1 << stage;
                    const int half = len / 2;
                    const int step = n / len;
                    for (int b = cursor; b < cursor + take; ++b) {
                        // half is 2^(stage - 1); keep the butterfly index
                        // arithmetic as shifts/masks in this hot loop.
                        const int group = b >> (stage - 1);
                        const int j = b & (half - 1);
                        const int i = group * len;
                        const auto& tw = plan.twiddles_[static_cast<std::size_t>(j * step)];
                        const auto w = tw;
                        const auto u = data[i + j];
                        const auto v = data[i + j + half] * w;
                        data[i + j] = u + v;
                        data[i + j + half] = u - v;
                    }
                    cursor += take;
                    spent += static_cast<long>(take) * kButterflyCost;
                    if (cursor >= total) {
                        cursor = 0;
                        ++stage;
                        if (stage > plan.log2_size_)
                            phase = inverse ? Phase::scale : Phase::idle;
                    }
                    break;
                }
                case Phase::conjugate: {
                    const int take = affordable(budget - spent, kConjugateCost, n - cursor);
                    if (take <= 0) return spent;
                    for (int i = cursor; i < cursor + take; ++i)
                        data[i] = std::conj(data[i]);
                    cursor += take;
                    spent += static_cast<long>(take) * kConjugateCost;
                    if (cursor >= n) { phase = Phase::permute; cursor = 0; }
                    break;
                }
                case Phase::permute: {
                    const int take = affordable(budget - spent, kGatherCost, n - cursor);
                    if (take <= 0) return spent;
                    for (int i = cursor; i < cursor + take; ++i) {
                        const int j = plan.reversed_[static_cast<std::size_t>(i)];
                        if (i < j) std::swap(data[i], data[j]);
                    }
                    cursor += take;
                    spent += static_cast<long>(take) * kGatherCost;
                    if (cursor >= n) { phase = Phase::butterflies; cursor = 0; stage = 1; }
                    break;
                }
                case Phase::scale: {
                    const int take = affordable(budget - spent, kScaleCost, n - cursor);
                    if (take <= 0) return spent;
                    const SampleType inv_n = SampleType{1} / static_cast<SampleType>(n);
                    for (int i = cursor; i < cursor + take; ++i)
                        data[i] = std::conj(data[i]) * inv_n;
                    cursor += take;
                    spent += static_cast<long>(take) * kScaleCost;
                    if (cursor >= n) phase = Phase::idle;
                    break;
                }
                case Phase::idle:
                    return spent;
                }
            }
            return spent;
        }
 
    private:
        static int affordable(long budget, long unit_cost, int remaining) {
            if (budget <= 0 || remaining <= 0) return 0;
            return static_cast<int>(std::min<long>(budget / unit_cost, remaining));
        }
    };

    void transform_real_padded(const double* src, int src_len,
                               std::complex<SampleType>* dst) const {
        for (int i = 0; i < size_; ++i) {
            const int r = reversed_[static_cast<std::size_t>(i)];
            const double v = r < src_len ? src[r] : 0.0;
            dst[i] = {static_cast<SampleType>(v), SampleType{0}};
        }
        for (int stage = 1; stage <= log2_size_; ++stage) {
            const int len = 1 << stage;
            const int half = len / 2;
            const int step = size_ / len;
            for (int i = 0; i < size_; i += len) {
                for (int j = 0; j < half; ++j) {
                    const auto& tw = twiddles_[static_cast<std::size_t>(j * step)];
                    const auto w = tw;
                    const auto u = dst[i + j];
                    const auto v = dst[i + j + half] * w;
                    dst[i + j] = u + v;
                    dst[i + j + half] = u - v;
                }
            }
        }
    }
 
private:
    int size_ = 0;
    int log2_size_ = 0;
    std::vector<int> reversed_;
    std::vector<std::complex<SampleType>> twiddles_;
};
 
}  // namespace detail

struct ZeroLatencyConvolverDesign {
    // ── Design parameters ────────────────────────────────────────────────────
    // Every value below is ours to tune and is declared with default + range.

    static constexpr int kHeadLenMin = 64;
    static constexpr int kHeadLenMax = 512;

    static constexpr double kIrGainDbDefault = 0.0;
    static constexpr double kIrGainDbMin = -24.0;
    static constexpr double kIrGainDbMax = 24.0;

    static constexpr double kPredelayMsDefault = 0.0;
    static constexpr double kPredelayMsMin = 0.0;
    static constexpr double kPredelayMsMax = 250.0;

    static constexpr double kWetPercentDefault = 100.0;
    static constexpr double kDryPercentDefault = 0.0;

    static constexpr double kLowcutHzDefault = 20.0;
    static constexpr double kLowcutHzMin = 20.0;
    static constexpr double kLowcutHzMax = 500.0;
    static constexpr double kHighcutHzDefault = 20000.0;
    static constexpr double kHighcutHzMin = 1000.0;
    static constexpr double kHighcutHzMax = 20000.0;

    static constexpr double kWidthPercentDefault = 100.0;
    static constexpr double kWidthPercentMin = 0.0;
    static constexpr double kWidthPercentMax = 200.0;

    static constexpr double kTailTrimDbDefault = -72.0;
    static constexpr double kTailTrimDbMin = -90.0;
    static constexpr double kTailTrimDbMax = -40.0;

    static constexpr double kTailFadeMsDefault = 5.0;
    static constexpr double kTailFadeMsMin = 2.0;
    static constexpr double kTailFadeMsMax = 20.0;

    static constexpr int kResampTapsPerPhaseDefault = 32;
    static constexpr int kResampTapsPerPhaseMin = 16;
    static constexpr int kResampTapsPerPhaseMax = 64;

    static bool valid_resample_geometry(int source_length, double source_rate,
                                        double destination_rate, int taps_per_phase) noexcept {
        if (source_length <= 0 || taps_per_phase <= 0 ||
            !std::isfinite(source_rate) || source_rate <= 0.0 ||
            !std::isfinite(destination_rate) || destination_rate <= 0.0)
            return false;
        constexpr double int_max = static_cast<double>(std::numeric_limits<int>::max());
        const double ratio = destination_rate / source_rate;
        if (!std::isfinite(ratio) || ratio <= 0.0) return false;
 
        const double destination_length =
            std::ceil(static_cast<double>(source_length) * ratio);
        if (!std::isfinite(destination_length) || destination_length < 1.0 ||
            destination_length > int_max)
            return false;
 
        const double cutoff = std::min(1.0, ratio);
        const double half_width = 0.5 * static_cast<double>(taps_per_phase) / cutoff;
        const double reach = std::ceil(half_width) + 1.0;
        if (!std::isfinite(half_width) || !std::isfinite(reach) || reach < 1.0 ||
            reach > int_max)
            return false;
 
        // The kernel loop forms `i0 +/- reach` in signed-int arithmetic. Its
        // largest possible centre is the final source sample; reserve one more
        // value so the loop increment after the last iteration cannot overflow.
        return reach <= int_max - static_cast<double>(source_length);
    }

    static constexpr double kResampKaiserBeta = 8.0;

    static constexpr double kBlockCostSlackDefault = 1.5;
    static constexpr double kBlockCostSlackMin = 1.2;
    static constexpr double kBlockCostSlackMax = 2.0;

    static constexpr double kSilenceEpsilon = 1e-24;
 
 
};
 
namespace detail {
 
// ── Partition group: two length-L partitions covering IR [2L, 4L) ────────
template <typename SampleType>
struct ZeroLatencyConvolverGroup {
    using FftCursor = typename SlicedFftPlanT<SampleType>::Cursor;
    int block_len = 0;       // L
    int fft_size = 0;        // K = 2L
    int ir_start = 0;        // d of the first partition = 2L
    int num_partitions = 0;  // 1 or 2
    int slot = 0;            // FDL slot the in-flight frame writes to
 
    detail::SlicedFftPlanT<SampleType> plan;
 
    // Per audio channel: FDL of forward spectra, depth num_partitions.
    std::vector<std::vector<std::complex<SampleType>>> spectra;     // [ch*np + p][K]
    // Per cell: IR partition spectra, accumulator, ping-pong output.
    std::vector<std::vector<std::complex<SampleType>>> ir_spectra;  // [cell*np + p][K]
    std::vector<std::vector<std::complex<SampleType>>> accum;       // [cell][K]
    std::vector<std::vector<SampleType>> emit;                      // [cell][L]
    std::vector<std::vector<SampleType>> ready;                     // [cell][L]
 
    // Sliced-execution state for the period currently in flight.
    long total_cost = 0;
    long granted = 0;
    std::size_t frame_base = 0;  // ring index of the closed frame's first sample
    int job = 0;
    int job_step = 0;
    int job_cursor = 0;
    FftCursor fft;
    bool armed = false;
    bool has_result = false;
 
    void reset() {
        for (auto& s : spectra) std::fill(s.begin(), s.end(), std::complex<SampleType>{});
        for (auto& a : accum) std::fill(a.begin(), a.end(), std::complex<SampleType>{});
        for (auto& e : emit) std::fill(e.begin(), e.end(), SampleType{0});
        for (auto& r : ready) std::fill(r.begin(), r.end(), SampleType{0});
        slot = 0;
        granted = 0;
        frame_base = 0;
        job = 0;
        job_step = 0;
        job_cursor = 0;
        fft.cancel();
        armed = false;
        has_result = false;
    }
};
 
template <typename SampleType>
struct ZeroLatencyConvolverCell {
    int src = 0;
    int dst = 0;
    int ir_channel = 0;
    double l1 = 0.0;
    std::vector<SampleType> head; // N0 direct-FIR taps; head[t] weighs x[n - t]
    // head reversed (oldest-first) for pulp::simd::correlate over a linear
    // window of the input history.
    std::vector<SampleType> head_reversed;
};
 
struct ZeroLatencyConvolverOnePole {
    double z = 0.0;
    double prev_x = 0.0;
};
 
}  // namespace detail
 
}  // namespace pulp::signal