modal_bank.hpp

Go to the documentation of this file.

#pragma once
 
#include <pulp/signal/denormal.hpp>
#include <pulp/signal/simd_buffer.hpp>
 
#include <algorithm>
#include <cmath>
#include <span>
 
namespace pulp::signal {

struct ModalMode {
    float freq_hz = 440.0f;
    float t60_s = 1.0f;   
    float gain = 1.0f;    
};

template <typename SampleType = float>
inline void fill_strike_pulse(std::span<SampleType> dst, SampleType velocity) {
    const auto n = dst.size();
    if (n == 0) return;
    if (n == 1) { dst[0] = velocity; return; }
    constexpr SampleType two_pi = SampleType(6.28318530717958647692);
    const SampleType phase_step = two_pi / static_cast<SampleType>(n);
    for (std::size_t i = 0; i < n; ++i) {
        dst[i] = velocity * SampleType(0.5)
               * (SampleType(1) - std::cos(phase_step * static_cast<SampleType>(i)));
    }
}

struct ModalShape {
    float half_waves = 1.0f;  
    float end_phase = 0.0f;   
};

inline float mode_shape_at(const ModalShape& shape, float position) {
    constexpr double pi = 3.14159265358979323846;
    return static_cast<float>(std::sin(pi * static_cast<double>(shape.half_waves)
                                            * static_cast<double>(position)
                                       + static_cast<double>(shape.end_phase)));
}

inline void fill_ideal_string_shapes(std::span<ModalShape> dst) {
    for (std::size_t i = 0; i < dst.size(); ++i)
        dst[i] = ModalShape{static_cast<float>(i + 1), 0.0f};
}

template <typename SampleType = float>
inline void fill_mode_weights(std::span<SampleType> dst,
                              std::span<const ModalShape> shapes, float position) {
    const std::size_t n = std::min(dst.size(), shapes.size());
    for (std::size_t i = 0; i < n; ++i)
        dst[i] = static_cast<SampleType>(mode_shape_at(shapes[i], position));
    for (std::size_t i = n; i < dst.size(); ++i) dst[i] = SampleType{0};
}

template <typename SampleType = float>
class ModalBankT {
public:
    static constexpr int kLanes = 16;

    static constexpr int kMaxPickups = 8;

    void prepare(double sample_rate, int max_modes, int max_pickups = 1) {
        sample_rate_ = sample_rate > 0.0 ? sample_rate : 48000.0;
        max_modes_ = max_modes > 0 ? max_modes : 0;
        max_pickups_ = std::clamp(max_pickups, 1, kMaxPickups);
        padded_ = static_cast<std::size_t>((max_modes_ + kLanes - 1) / kLanes) * kLanes;
        cu_.resize(padded_);
        cv_.resize(padded_);
        rc_.resize(padded_);
        rs_.resize(padded_);
        u_.resize(padded_);
        v_.resize(padded_);
        pickup_.resize(padded_ * static_cast<std::size_t>(max_pickups_));
        std::fill(pickup_.begin(), pickup_.end(), SampleType{1});
        mode_count_ = 0;
    }

    void set_modes(std::span<const ModalMode> modes) {
        const int count = static_cast<int>(
            std::min<std::size_t>(modes.size(), static_cast<std::size_t>(max_modes_)));
        const double fs = sample_rate_;
        for (int m = 0; m < count; ++m) {
            const auto i = static_cast<std::size_t>(m);
            const double freq = std::clamp(static_cast<double>(modes[i].freq_hz),
                                           1.0, 0.49 * fs);
            const double t60 = std::max(static_cast<double>(modes[i].t60_s), 1.0e-3);
            const double r = std::pow(10.0, -3.0 / (t60 * fs));
            const double w = 2.0 * 3.14159265358979323846 * freq / fs;
            const double cos_w = std::cos(w);
            const double sin_w = std::sin(w);
            const double gain = static_cast<double>(modes[i].gain);
            rc_[i] = static_cast<SampleType>(r * cos_w);
            rs_[i] = static_cast<SampleType>(r * sin_w);
            cu_[i] = static_cast<SampleType>(gain * cos_w);
            cv_[i] = static_cast<SampleType>(gain * sin_w);
        }
        // Zero the coefficients of unused padding lanes so they contribute
        // exact silence regardless of stale state.
        for (std::size_t m = static_cast<std::size_t>(count); m < padded_; ++m) {
            cu_[m] = SampleType{0};
            cv_[m] = SampleType{0};
            rc_[m] = SampleType{0};
            rs_[m] = SampleType{0};
        }
        mode_count_ = count;
    }

    void set_pickup_gains(int channel, std::span<const SampleType> gains) {
        if (channel < 0 || channel >= max_pickups_) return;
        const std::size_t base = static_cast<std::size_t>(channel) * padded_;
        const std::size_t n = std::min(gains.size(), padded_);
        for (std::size_t m = 0; m < n; ++m) pickup_[base + m] = gains[m];
    }

    void process_add(const SampleType* input, SampleType* output, int num_samples,
                     int channel = 0) {
        SampleType* outs[1] = {output};
        const int c = std::clamp(channel, 0, max_pickups_ - 1);
        run(input, outs, &c, 1, num_samples);
    }

    void process_add(const SampleType* input, SampleType* const* outputs,
                     int num_channels, int num_samples) {
        int channels[kMaxPickups];
        const int n = std::clamp(num_channels, 0, max_pickups_);
        for (int c = 0; c < n; ++c) channels[c] = c;
        run(input, outputs, channels, n, num_samples);
    }

    void reset() {
        u_.clear();
        v_.clear();
    }
 
    int mode_count() const { return mode_count_; }
    int max_modes() const { return max_modes_; }
    int max_pickups() const { return max_pickups_; }
    double sample_rate() const { return sample_rate_; }
 
private:
    void run(const SampleType* input, SampleType* const* outputs,
             const int* channels, int num_channels, int num_samples) {
        if (num_channels <= 0) return;
        const int groups = (mode_count_ + kLanes - 1) / kLanes;
        for (int g = 0; g < groups; ++g) {
            const std::size_t base = static_cast<std::size_t>(g) * kLanes;
            SampleType cu[kLanes], cv[kLanes], rc[kLanes], rs[kLanes];
            SampleType u[kLanes], v[kLanes];
            SampleType pick[kMaxPickups][kLanes];
            for (int l = 0; l < kLanes; ++l) {
                const std::size_t i = base + static_cast<std::size_t>(l);
                cu[l] = cu_[i];
                cv[l] = cv_[i];
                rc[l] = rc_[i];
                rs[l] = rs_[i];
                u[l] = u_[i];
                v[l] = v_[i];
            }
            for (int c = 0; c < num_channels; ++c) {
                const std::size_t pbase =
                    static_cast<std::size_t>(channels[c]) * padded_ + base;
                for (int l = 0; l < kLanes; ++l) pick[c][l] = pickup_[pbase + l];
            }
            for (int i = 0; i < num_samples; ++i) {
                const SampleType x = input[i];
                SampleType y[kLanes];
                for (int l = 0; l < kLanes; ++l) {
                    const SampleType nu = rc[l] * u[l] - rs[l] * v[l] + cu[l] * x;
                    const SampleType nv = rs[l] * u[l] + rc[l] * v[l] + cv[l] * x;
                    u[l] = nu;
                    v[l] = nv;
                    y[l] = nv;
                }
                for (int c = 0; c < num_channels; ++c) {
                    SampleType t[kLanes];
                    for (int l = 0; l < kLanes; ++l) t[l] = pick[c][l] * y[l];
                    // Fixed pairwise reduction order: deterministic output and
                    // vector-width friendly (elementwise adds, then one scalar).
                    for (int step = kLanes / 2; step > 0; step /= 2)
                        for (int l = 0; l < step; ++l) t[l] += t[l + step];
                    outputs[c][i] += t[0];
                }
            }
            for (int l = 0; l < kLanes; ++l) {
                const std::size_t i = base + static_cast<std::size_t>(l);
                u_[i] = snap_to_zero(u[l]);
                v_[i] = snap_to_zero(v[l]);
            }
        }
    }
 
    AlignedBufferT<SampleType> cu_, cv_, rc_, rs_, u_, v_, pickup_;
    double sample_rate_ = 48000.0;
    std::size_t padded_ = 0;
    int max_modes_ = 0;
    int max_pickups_ = 1;
    int mode_count_ = 0;
};
 
using ModalBank = ModalBankT<float>;
using ModalBank64 = ModalBankT<double>;
 
} // namespace pulp::signal