convolver_non_uniform.hpp

Go to the documentation of this file.

#pragma once

 
#include <pulp/signal/convolver.hpp>
#include <pulp/signal/convolver_messages.hpp>
 
#include <algorithm>
#include <array>
#include <cassert>
#include <complex>
#include <cstddef>
#include <cstdint>
#include <string_view>
#include <type_traits>
#include <vector>
 
#if defined(__APPLE__)
#include <TargetConditionals.h>
#include <sys/sysctl.h>
#endif
 
// The tail-multiplier table was measured on Apple Silicon Macs with the
// Accelerate backend in single precision; it applies only there.
#if defined(PULP_SIMD_DEFAULT_BACKEND_ACCELERATE) && defined(__aarch64__) &&                       \
    defined(TARGET_OS_OSX) && TARGET_OS_OSX
#define PULP_NON_UNIFORM_MEASURED_TABLE 1
#else
#define PULP_NON_UNIFORM_MEASURED_TABLE 0
#endif
 
namespace pulp::signal {

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

    static constexpr std::size_t kDefaultTailMultiplier = 4;

    static constexpr std::size_t choose_tail_multiplier(std::size_t ir_length,
                                                        std::size_t block_size) noexcept {
#if PULP_NON_UNIFORM_MEASURED_TABLE
        if constexpr (!std::is_same_v<SampleType, float>) {
            (void)ir_length;
            (void)block_size;
            return kDefaultTailMultiplier;
        } else {
            // Columns: IRs nearest 4096, 16384, 48000 and 96000 taps.
            constexpr std::array<std::uint64_t, 4> kMeasuredIr{4096, 16384, 48000, 96000};
            std::size_t column = 0;
            const auto n = static_cast<std::uint64_t>(ir_length);
            while (column + 1 < kMeasuredIr.size() &&
                   n * n > kMeasuredIr[column] * kMeasuredIr[column + 1])
                ++column;
            switch (detail::convolver_block_size(block_size)) {
            case 16:
                return std::array<std::size_t, 4>{16, 32, 64, 64}[column];
            case 32:
                return std::array<std::size_t, 4>{8, 16, 32, 64}[column];
            case 64:
                return std::array<std::size_t, 4>{8, 16, 16, 32}[column];
            case 128:
                return std::array<std::size_t, 4>{8, 8, 16, 32}[column];
            case 256:
                return std::array<std::size_t, 4>{4, 8, 16, 16}[column];
            case 512:
                return std::array<std::size_t, 4>{8, 8, 8, 16}[column];
            default:
                return kDefaultTailMultiplier;
            }
        }
#else
        (void)ir_length;
        (void)block_size;
        return kDefaultTailMultiplier;
#endif
    }

    static constexpr bool table_measured_on(std::string_view cpu_brand) noexcept {
        for (const std::string_view family :
             {std::string_view{"Apple M1"}, std::string_view{"Apple M5"}}) {
            if (cpu_brand.substr(0, family.size()) != family)
                continue;
            if (cpu_brand.size() == family.size() || cpu_brand[family.size()] == ' ')
                return true;
        }
        return false;
    }

    static constexpr std::size_t default_tail_multiplier_for(std::string_view cpu_brand,
                                                             std::size_t ir_length,
                                                             std::size_t block_size) noexcept {
        return table_measured_on(cpu_brand) ? choose_tail_multiplier(ir_length, block_size)
                                            : kDefaultTailMultiplier;
    }

    static std::size_t default_tail_multiplier(std::size_t ir_length, std::size_t block_size) {
#if PULP_NON_UNIFORM_MEASURED_TABLE
        static const bool measured = [] {
            char brand[128] = {};
            std::size_t size = sizeof(brand) - 1;
            if (sysctlbyname("machdep.cpu.brand_string", brand, &size, nullptr, 0) != 0)
                return false;
            return table_measured_on(std::string_view{brand});
        }();
        if (measured)
            return choose_tail_multiplier(ir_length, block_size);
#endif
        (void)ir_length;
        (void)block_size;
        return kDefaultTailMultiplier;
    }

    void load_ir(const SampleType* ir, std::size_t ir_length, std::size_t block_size) {
        load_ir(ir, ir_length, block_size, default_tail_multiplier(ir_length, block_size));
    }

    void load_ir(const SampleType* ir, std::size_t ir_length, std::size_t block_size,
                 std::size_t tail_multiplier) {
        reset_state();
        if (ir == nullptr || ir_length == 0 || block_size == 0) return;
 
        // Pow-of-2 block size.
        block_size_ = round_up_pow2(block_size);
        if (tail_multiplier == 0) tail_multiplier = 1;
        // Tail multiplier must be a power of 2 so the larger FFT also
        // works on a pow-of-2 block.
        tail_multiplier_ = round_up_pow2(tail_multiplier);
 
        std::size_t head_samples = tail_multiplier_ * block_size_;
        head_samples = std::min(head_samples, ir_length);
 
        // ── head stage ──
        head_.load_ir(ir, head_samples, block_size_);
        head_samples_ = head_samples;
 
        // ── tail stage (optional) ──
        if (head_samples < ir_length) {
            std::size_t tail_block = block_size_ * tail_multiplier_;
            std::size_t tail_len = ir_length - head_samples;
            // The tail IR is the slice ir[head_samples .. ir_length)
            // with NO leading zero-padding. The K-small-block input
            // buffering delay equals the K·B = head_samples IR-tap
            // offset, so the natural sum lines up.
            tail_ir_slice_.assign(ir + head_samples, ir + ir_length);
            tail_.load_ir(tail_ir_slice_.data(), tail_len, tail_block);
            tail_block_ = tail_block;
            // Input/output ring buffers for the tail stage.
            tail_input_buf_.assign(tail_block_, SampleType{0.0f});
            tail_output_buf_.assign(tail_block_, SampleType{0.0f});
            tail_fill_ = 0;
            tail_stream_pos_ = tail_block_;  // nothing to stream yet
        }
        loaded_ = true;
    }

    void process(const SampleType* input,
                 SampleType* output,
                 std::size_t num_samples) {
        // No IR to convolve with: pass through.
        if (!loaded_) {
            if (input != output)
                std::copy_n(input, num_samples, output);
            return;
        }
        // Loaded, but the caller's block size does not match the partitioning.
        if (num_samples != 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();
 
        // ── feed this block into the tail input buffer ──
        // Captured before the head writes `output`, so `input` may be the
        // same buffer as `output`.
        if (tail_block_ > 0) {
            std::copy_n(input, block_size_, tail_input_buf_.data() + tail_fill_);
            tail_fill_ += block_size_;
        }
 
        // ── head stage: live ──
        head_.process(input, output, num_samples);
 
        if (tail_block_ == 0) return;
 
        // ── stream out previously-computed tail samples ──
        if (tail_stream_pos_ < tail_block_) {
            for (std::size_t i = 0; i < block_size_; ++i)
                output[i] += tail_output_buf_[tail_stream_pos_ + i];
            tail_stream_pos_ += block_size_;
        }
 
        if (tail_fill_ >= tail_block_) {
            // Process one tail block worth of input. Result will be
            // streamed out over the NEXT K small blocks (starting on
            // the next process() call).
            tail_.process(tail_input_buf_.data(),
                          tail_output_buf_.data(), tail_block_);
            tail_fill_ = 0;
            tail_stream_pos_ = 0;
        }
    }
 
    void reset() {
        head_.reset();
        if (tail_block_ > 0) tail_.reset();
        std::fill(tail_input_buf_.begin(), tail_input_buf_.end(), SampleType{0.0f});
        std::fill(tail_output_buf_.begin(), tail_output_buf_.end(), SampleType{0.0f});
        tail_fill_ = 0;
        tail_stream_pos_ = tail_block_;
        history_torn_ = false;
        block_size_violations_ = 0;
    }

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

    std::uint64_t block_size_violations() const { return block_size_violations_; }
 
    std::size_t head_samples() const { return head_samples_; }
    std::size_t tail_block() const { return tail_block_; }
    std::size_t tail_multiplier() const { return tail_multiplier_; }
 
    bool is_loaded() const { return loaded_; }
 
private:
    static std::size_t round_up_pow2(std::size_t v) {
        if (v <= 1) return 1;
        if ((v & (v - 1)) == 0) return v;
        std::size_t p = 1;
        while (p < v) p <<= 1;
        return p;
    }
 
    void reset_state() {
        head_ = PartitionedConvolverT<SampleType>{};
        tail_ = PartitionedConvolverT<SampleType>{};
        tail_ir_slice_.clear();
        tail_input_buf_.clear();
        tail_output_buf_.clear();
        block_size_ = 0;
        tail_block_ = 0;
        tail_multiplier_ = 0;
        head_samples_ = 0;
        tail_fill_ = 0;
        tail_stream_pos_ = 0;
        loaded_ = false;
    }
 
    PartitionedConvolverT<SampleType> head_;
    PartitionedConvolverT<SampleType> tail_;
    std::vector<SampleType> tail_ir_slice_;
    std::vector<SampleType> tail_input_buf_;
    std::vector<SampleType> tail_output_buf_;
    std::size_t block_size_ = 0;
    std::size_t tail_block_ = 0;
    std::size_t tail_multiplier_ = 0;
    std::size_t head_samples_ = 0;
    std::size_t tail_fill_ = 0;
    std::size_t tail_stream_pos_ = 0;  // tail_block_ ⇒ nothing to stream
    bool loaded_ = false;
    // Set when process() is handed a block whose size does not match the loaded
    // IR partitioning. The next valid block resets both stages before folding
    // itself in.
    bool history_torn_ = false;
    std::uint64_t block_size_violations_ = 0;
};
 
using NonUniformPartitionedConvolver = NonUniformPartitionedConvolverT<float>;
using NonUniformPartitionedConvolver64 = NonUniformPartitionedConvolverT<double>;
 
} // namespace pulp::signal