convolver.hpp

Go to the documentation of this file.

#pragma once
 
#include <pulp/signal/convolver_messages.hpp>
#include <pulp/signal/fft.hpp>
#include <pulp/signal/transition_mixer.hpp>
 
#include <algorithm>
#include <cassert>
#include <complex>
#include <cstddef>
#include <cstdint>
#include <memory>
#include <vector>
 
namespace pulp::signal {

template <typename SampleType = float>
class PartitionedConvolverT {
public:
    PartitionedConvolverT() = default;

    void load_ir(const SampleType* ir, std::size_t ir_length, std::size_t block_size) {
        state_ = detail::build_convolver_ir_state(ir, ir_length, block_size);
        // Adopt the freshly-built state's spare as THE input history. Loading an
        // IR outright (as opposed to swapping one in) deliberately starts the
        // stream cold — this is a prepare()-time call, not a live change.
        history_ = state_ ? std::move(state_->history) : nullptr;
    }

    bool try_swap_ir(ConvolverIrSwapperT<SampleType>& swapper) {
        // Gate the swap on retire-ring capacity FIRST so the audio
        // thread never has to free the displaced IR inline if the
        // ring is saturated. A single retired slot would let the audio
        // thread deallocate when 2+ swaps happened between drain_old()
        // calls. Refusing the swap when the ring is full is RT-safe:
        // pending stays in the swapper for the next try_swap_ir attempt;
        // in-flight state_ continues uninterrupted; worker thread catches
        // up on its next drain tick.
        if (state_ && !swapper.has_retire_capacity()) {
            return false;
        }
 
        // A crossfade still running: don't start another swap and
        // don't touch the in-flight fade-out — the pending IR waits for the next
        // call. Human-paced IR swaps rarely overlap the ~ms fade window.
        if (fading_ && !fade_mixer_.done())
            return false;
        // A completed fade-out can now be parked for the worker to free (never
        // on the audio thread). Refuse the swap if the ring is full so we never
        // free inline.
        if (fading_) {
            // Reserve TWO slots: one for this completed fade-out, one for the
            // just-displaced IR should the incoming swap fall back to an instant
            // retire below. Without this a near-full ring could accept the
            // fade-out then fail the second retire and strand (leak) the
            // displaced IR on the audio thread.
            if (!swapper.has_retire_capacity(2)) return false;
            (void)swapper.retire(fading_);
            (void)fading_.release();
        }
 
        auto next = swapper.try_consume();
        if (!next)
            return false;
 
        auto previous = std::move(state_);
        state_ = std::move(next);
 
        // The input delay line does NOT move with the swap. It belongs to the
        // input stream, whose identity is unchanged by an IR change, so the
        // convolver keeps the one it has and renders the incoming IR — and, under
        // a crossfade, the outgoing one too — against that same history. Only the
        // ring's SHAPE can force a change: a first load, a different block/FFT
        // size, or an IR with more partitions than the ring has slots.
        auto spare = std::move(state_->history);
        const bool ring_serves_new =
            history_ != nullptr &&
            history_->block_size == state_->block_size &&
            history_->fft_size == state_->fft_size &&
            history_->spectra.size() >= state_->num_partitions;
        if (!ring_serves_new) {
            // Adopt the incoming state's spare ring — allocated off the audio
            // thread precisely so this step never has to — carrying across
            // whatever the old ring held. If the incoming state has no spare (a
            // hand-built state that skipped it) the old ring is kept; process()
            // detects the too-short ring and fails safe rather than aliasing it.
            if (spare) {
                if (history_) detail::migrate_input_history(*history_, *spare);
                history_.swap(spare);
            }
        }
 
        if (previous) {
            // Crossfade path: keep the displaced IR rendering in parallel from the
            // shared history and blend old->new. Only when a fade is configured AND
            // the scratch fits both IRs' block size; otherwise fall back to the
            // instant retire below. Require MATCHING block sizes: the fade-out
            // render reads a full `previous->block_size` of `input`, but process()
            // only guarantees num_samples == state_->block_size — a larger old
            // block would read past the caller's input buffer. A mismatch falls
            // back to instant.
            const bool can_fade =
                fade_samples_ > 0 &&
                previous->block_size == state_->block_size &&
                fade_scratch_.size() >= static_cast<std::size_t>(state_->block_size);
            // Whichever ring lost the swap above is now rendered by nobody. Park
            // it on the displaced state so it is reclaimed off the audio thread
            // along with it — the audio thread must never free. The live state
            // never holds a spare (it is taken above on every path that installs
            // one), so this overwrites nothing and frees nothing.
            assert(!previous->history);
            previous->history = std::move(spare);
            if (can_fade) {
                fading_ = std::move(previous);
                fade_mixer_.configure(fade_samples_, fade_curve_);
            } else {
                // Pre-checked capacity above; the false branch is defense-in-
                // depth — leaks rather than freeing on RT.
                const bool ok = swapper.retire(previous);
                (void)ok;
                (void)previous.release();
            }
        } else {
            // `previous` is null only on a first load into an unloaded convolver,
            // where `history_` was null too — so the adopt branch above left
            // `spare` empty and nothing is freed on the audio thread here.
            assert(!spare);
        }
        return true;
    }

    void set_crossfade(std::size_t samples, TransitionCurve curve = TransitionCurve::Smoothstep) {
        fade_samples_ = samples;
        fade_curve_ = curve;
        if (state_) {
            fade_scratch_.assign(static_cast<std::size_t>(state_->block_size),
                                 SampleType{0.0f});
        }
    }

    void process(const SampleType* input,
                 SampleType* output,
                 std::size_t num_samples) {
        // No IR to convolve with — or no input history long enough to serve it,
        // which only a hand-built state that skipped its `history` spare can
        // produce. Either way the convolver is a wire rather than something that
        // would read ring slots that do not exist.
        if (!state_ || state_->ir_spectra.empty() || !history_ ||
            history_->block_size != state_->block_size ||
            history_->fft_size != state_->fft_size ||
            history_->spectra.size() < state_->num_partitions) {
            std::copy_n(input, num_samples, output);
            return;
        }
        // Loaded, but the caller's block size does not match the IR partitioning.
        if (static_cast<int>(num_samples) != state_->block_size) {
            std::fill_n(output, num_samples, SampleType{0.0f});
            ++block_size_violations_;
            history_torn_ = true;
            return;
        }
        // Recover from a previous violation before folding a valid block in.
        if (history_torn_) reset();
 
        // Fold this block into the shared input history ONCE — one forward FFT
        // per block however many IRs render from it.
        capture_input(*history_, *state_->fft, input);
 
        render_ir(*state_, *history_, output);
 
        // Crossfade the retiring IR: render it in parallel from the SAME history
        // into scratch and blend old→new. RT-safe: scratch is pre-sized, no
        // alloc; the fade-out IR is retired on the ring after the fade (or
        // opportunistically once done), never freed on the audio thread.
        if (fading_ && !fade_mixer_.done()) {
            // Everything the parallel render needs. All three hold for any fade
            // try_swap_ir() actually starts — it only begins one when the scratch
            // fits and the block sizes match, and the shared ring is always at
            // least as long as the outgoing IR needs.
            const bool can_render =
                fade_scratch_.size() >= num_samples &&
                static_cast<std::size_t>(fading_->block_size) == num_samples &&
                history_->spectra.size() >= fading_->num_partitions;
            if (can_render) {
                render_ir(*fading_, *history_, fade_scratch_.data());
                const std::size_t base = fade_mixer_.position();
                for (std::size_t i = 0; i < num_samples; ++i) {
                    SampleType old_gain, new_gain;
                    fade_mixer_.gains_at(base + i, old_gain, new_gain);
                    output[i] = output[i] * new_gain + fade_scratch_[i] * old_gain;
                }
            }
            // Advance even when the render was skipped. A fade that cannot render
            // must still COMPLETE: `fading_` is only retired once the mixer is
            // done, and try_swap_ir() refuses every swap while a fade is running,
            // so a fade that never finishes would wedge the convolver on the IR
            // it happens to be holding.
            fade_mixer_.advance(num_samples);
        }
        // When the fade completes the fade-out IR is simply held (not rendered
        // anymore); it is parked on the swapper's retire ring by the next
        // try_swap_ir() — never freed on the audio thread here.
 
        // Step the shared ring once, after every IR has read this block.
        advance_input(*history_);
    }

    void reset() {
        clear_history(history_.get());
        history_torn_ = false;
        block_size_violations_ = 0;
    }

    std::size_t latency() const { return 0; }

    std::size_t block_size() const {
        return state_ ? static_cast<std::size_t>(state_->block_size) : 0;
    }

    std::uint64_t block_size_violations() const { return block_size_violations_; }
 
    std::size_t num_partitions() const {
        return state_ ? state_->num_partitions : 0;
    }
 
    bool is_loaded() const {
        return state_ && !state_->ir_spectra.empty();
    }
 
  private:
    // Fold one input block into the shared history: copy it into the upper
    // half of the real overlap window and take its real forward transform into
    // the current ring slot (DC..Nyquist only; every spectrum here is the
    // transform of real data, so the upper half is never needed). Called ONCE
    // per processed block, before any IR renders.
    static void capture_input(ConvolverInputHistoryT<SampleType>& h,
                              FftT<SampleType>& fft,
                              const SampleType* input) {
        SampleType* window = overlap_samples(h);
        std::copy_n(input, h.block_size, window + h.block_size);
        // The transform is a property of the FFT size, which every IR rendering
        // against this history shares, so any loaded state's `fft` will do.
        fft.forward_real_half(window, h.spectra[h.write_pos].data());
    }
 
    // Slide the overlap window down and step the ring cursor. Called ONCE per
    // processed block, after every IR has read this block — a crossfade renders
    // twice from one history, so stepping per render would double-advance it.
    static void advance_input(ConvolverInputHistoryT<SampleType>& h) {
        SampleType* window = overlap_samples(h);
        std::copy_n(window + h.block_size, h.block_size, window);
        h.write_pos = (h.write_pos + 1) % h.spectra.size();
    }
 
    // The overlap window as `fft_size` real samples, packed into the leading
    // scalars of the history's complex `overlap` storage.
    // (std::complex<T> is layout-compatible with T[2], [complex.numbers].)
    static SampleType* overlap_samples(ConvolverInputHistoryT<SampleType>& h) {
        return reinterpret_cast<SampleType*>(h.overlap.data());
    }
 
    // Multiply-accumulate one IR's partitions against the shared input history
    // and write the block's time-domain result. Shared by the live path (state_)
    // and the parallel fade-out render (fading_). Reads the history; never writes
    // or advances it. Caller guarantees `s` has a loaded IR and that the ring is
    // at least `s.num_partitions` long.
    static void render_ir(ConvolverIrStateT<SampleType>& s,
                          const ConvolverInputHistoryT<SampleType>& h,
                          SampleType* output) {
        // Input and IR are transforms of real data, so every spectrum is
        // Hermitian and only DC..Nyquist (N/2+1 bins) is independent. Only those
        // bins are multiply-accumulated, and the real inverse transform below
        // reads only those bins.
        const int half = s.fft_size / 2;   // Nyquist bin; MAC bins [0, half]
        const std::size_t capacity = h.spectra.size();
        std::complex<SampleType>* accum = s.accum.data();
        const auto spectrum = [&](std::size_t p) {
            // Partition p convolves with the input block of age p. `write_pos`
            // holds the block being processed, so age p is p slots back.
            return h.spectra[(h.write_pos + capacity - p) % capacity].data();
        };
        // Complex MACs as explicit float ops: the finite-value result of
        // std::complex::operator* without its libm NaN/Inf-recovery branch, so
        // the loops vectorize. Matches operator* for the finite spectra a
        // convolver sees (differences are sub-ulp). Partition 0 writes rather
        // than accumulates, so the accumulator never needs clearing; the rest
        // go four at a time so each pass loads and stores the accumulator once
        // for four products.
        {
            const std::complex<SampleType>* in = spectrum(0);
            const std::complex<SampleType>* ir = s.ir_spectra[0].data();
            for (int i = 0; i <= half; ++i)
                accum[i] = cmul(in[i], ir[i]);
        }
        std::size_t p = 1;
        for (; p + 4 <= s.num_partitions; p += 4) {
            const std::complex<SampleType>* in0 = spectrum(p);
            const std::complex<SampleType>* in1 = spectrum(p + 1);
            const std::complex<SampleType>* in2 = spectrum(p + 2);
            const std::complex<SampleType>* in3 = spectrum(p + 3);
            const std::complex<SampleType>* ir0 = s.ir_spectra[p].data();
            const std::complex<SampleType>* ir1 = s.ir_spectra[p + 1].data();
            const std::complex<SampleType>* ir2 = s.ir_spectra[p + 2].data();
            const std::complex<SampleType>* ir3 = s.ir_spectra[p + 3].data();
            for (int i = 0; i <= half; ++i)
                accum[i] += (cmul(in0[i], ir0[i]) + cmul(in1[i], ir1[i])) +
                            (cmul(in2[i], ir2[i]) + cmul(in3[i], ir3[i]));
        }
        for (; p < s.num_partitions; ++p) {
            const std::complex<SampleType>* in = spectrum(p);
            const std::complex<SampleType>* ir = s.ir_spectra[p].data();
            for (int i = 0; i <= half; ++i)
                accum[i] += cmul(in[i], ir[i]);
        }
 
        // Real inverse in the accumulator's own storage, 1/N folded in; overlap-
        // save keeps the second half of the window.
        // std::complex<T> is layout-compatible with T[2] ([complex.numbers]),
        // so the accumulator's N/2+1 complex slots hold the N real samples.
        auto* time = reinterpret_cast<SampleType*>(accum);
        s.fft->inverse_real_half(accum, time);
        std::copy_n(time + s.block_size, s.block_size, output);
    }
 
    static std::complex<SampleType> cmul(std::complex<SampleType> a,
                                         std::complex<SampleType> b) noexcept {
        return {a.real() * b.real() - a.imag() * b.imag(),
                a.real() * b.imag() + a.imag() * b.real()};
    }
 
    // Zero the shared input history so the next block starts a new stream.
    // Null-tolerant so reset() can call it on an unloaded convolver.
    static void clear_history(ConvolverInputHistoryT<SampleType>* h) {
        if (!h) return;
        constexpr std::complex<SampleType> kZero{SampleType{0.0f}, SampleType{0.0f}};
        h->spectra.fill_zero();
        std::fill(h->overlap.begin(), h->overlap.end(), kZero);
        h->write_pos = 0;
    }
 
    std::unique_ptr<ConvolverIrStateT<SampleType>> state_;
    // The input delay line, shared by every IR this convolver renders. Outlives
    // any one IR: a swap replaces `state_`, never this. Its ring capacity is the
    // high-water mark of partition counts seen since the last `load_ir()`, since
    // growing it on the audio thread means adopting an incoming state's spare and
    // shrinking it would mean throwing history away for nothing.
    std::unique_ptr<ConvolverInputHistoryT<SampleType>> history_;
    // Set when process() is handed a block whose size does not match the loaded
    // IR partitioning. The next valid block resets before folding itself in.
    bool history_torn_ = false;
    std::uint64_t block_size_violations_ = 0;
 
    // ── Crossfade state (opt-in via set_crossfade) ────────────────────────────
    std::unique_ptr<ConvolverIrStateT<SampleType>> fading_; // IR fading out (parallel render)
    signal::TransitionMixerT<SampleType> fade_mixer_; // shared click-free blend
    std::size_t fade_samples_ = 0;                // configured fade length (0 = instant)
    TransitionCurve fade_curve_ = TransitionCurve::Smoothstep;
    std::vector<SampleType> fade_scratch_;        // parallel-render output, pre-sized off-RT
};
 
using PartitionedConvolver = PartitionedConvolverT<float>;
using PartitionedConvolver64 = PartitionedConvolverT<double>;
 
} // namespace pulp::signal