string_loop.hpp

Go to the documentation of this file.

// sbl/dsp/pm/string_loop.hpp — A string as two waveguide segments (Physical modeling — composition)
//
// The excitation point is the junction between two segments: one running to
// the nut, one to the bridge. Each delay line folds the round trip of its own
// segment, so their lengths sum to the period (STK's convention):
//
//     d_nut = T * beta        d_bridge = T * (1 - beta)
//
// AP-032 Phase 1 settled this the expensive way: a single loop with a tap
// cannot be bowed. DC is a fixed point there — the string converges on the
// bow velocity and sticks. Two segments with inverting reflections at both
// ends are what gives Helmholtz motion: the injection enters at the junction
// and travels both ways; the nut segment returns it inverted; the bridge
// segment passes through the allpass and the loss filter, then the bridge's
// reflection, and returns it to the junction. The wiring as compiled is a
// describe() away — `sbl-diagram waveguide-loop` (AP-037).
//
// Per sample the exciter is handed the junction's state — the velocity
// arriving there, the point's displacement, and the loop's wave impedance —
// and returns the force it exerts; the loop turns the force into the wave
// it injects, F / 2Z. Energy passes both ways at the point (a bow does work,
// a finger damps, a hammer bounces). A callable of one float, the velocity,
// is accepted too, for tests and plucks that need nothing else.
//
// The displacement is the point's velocity integrated, with a slow leak so a
// bowed string's small DC does not walk it off (DISPLACEMENT_LEAK_HZ). It is
// in wave units × samples; a displacement law scales it as it likes.
//
// The bridge end has the same shape (AP-036): a second callable is handed the
// wave arriving at the bridge, after the loop's loss and allpass, and returns
// the reflection. Without one the bridge is a rigid wall, reflection −1. A
// LoadJunction with a lumped load makes it a body.
//
// ## Tuning
//
// The delay elements sum to the period: the integer parts go in the delay
// lines and the leftover fraction goes to the Thiran allpass, so the loop
// tunes between samples rather than to the nearest one. A segment is never
// shorter than one sample; when the split would make one shorter, the other
// gives up the difference so the total still holds.
//
// The loss filter's delay is not taken out of the period, so the loop plays
// flat by it — more so the darker the damping and the shorter the string.
// That is deliberate for now (AP-032 addendum 2026-09-12): the string is
// heard as it is before anything is added to correct it.
//
// ## Energy
//
// Both reflections are -1 and the loss stage has gain at or below unity, so
// the loop is passive with no excitation: it can only lose energy, never
// make it (RPT-029 Part 4 rule 1). The single source is the injection the
// caller returns. The loop gain is the passivity margin — energy leaving
// into the body at the bridge — and a value below 1 is what makes a plucked
// string decay.
//
// Usage:
//   loop.init(nut_buf, 1024, bridge_buf, 1024);
//   loop.set_period_samples(sample_rate / hz);
//   loop.set_junction(0.3f);               // bow position along the string
//   float out = loop.tick([&](float v_string) {
//       return bow.process(bow_velocity - v_string);
//   });
//   // with a body:
//   float out = loop.tick(bow_fn, [&](float arriving) { return bridge.reflect(arriving); });
 
#ifndef SBL_DSP_PM_STRING_LOOP_HPP_
#define SBL_DSP_PM_STRING_LOOP_HPP_
 
#include <cmath>
#include <cstdint>
#include <type_traits>
 
#include <sbl/dsp/diagram/graph.hpp>
#include <sbl/dsp/math/clamp.hpp>
#include <sbl/dsp/math/constants.hpp>
#include <sbl/dsp/pm/contracts.hpp>
#include <sbl/dsp/primitives/delay_line.hpp>
#include <sbl/dsp/primitives/one_pole.hpp>
#include <sbl/dsp/primitives/thiran_allpass.hpp>
#include <sbl/dsp/types/fixed.hpp>
 
namespace sbl::dsp::pm {
 
class StringLoop {
public:
    /// @note All public methods are ISR-safe — bounded computation, no I/O.
 
    /// Shortest loop: eight samples, 6 kHz at 48 kHz. Pitch requests above it play this.
    static constexpr float MIN_PERIOD = 8.0f;
    /// The loop's wave impedance unless told otherwise: the unit the laws and loads are written in.
    static constexpr float DEFAULT_IMPEDANCE = 1.0f;
    /// The junction displacement forgets DC below this: a leaky integrator.
    static constexpr float DISPLACEMENT_LEAK_HZ = 5.0f;
 
    /**
     * @brief Install the two segments' buffers
     *
     * Each buffer bounds its own segment, so the longest period depends on
     * the junction: the nut segment needs `period * beta` samples and the
     * bridge segment the rest. A segment asked for more than its buffer is
     * clamped to the buffer, and the loop then plays sharp of the request;
     * `delay_a_samples()` / `delay_b_samples()` report the lengths
     * actually read, so a caller can tell.
     */
    void init(float* a_buffer, uint32_t a_max,
              float* b_buffer, uint32_t b_max) {
        end_a_.init(a_buffer, a_max);
        end_b_.init(b_buffer, b_max);
        // The delay line reads up to max - 1; the last slot is the write.
        limit_a_ = static_cast<float>(a_max > 1 ? a_max - 1 : 1);
        limit_b_ = static_cast<float>(b_max > 1 ? b_max - 1 : 1);
        loss_.set_range(DAMPING_MIN_HZ, DAMPING_MAX_HZ);
        loss_.set_cutoff_su(DEFAULT_DAMPING_SU);
        retune();
    }
 
    /// Loop period in samples — the pitch. Sample rate over frequency.
    void set_period_samples(float period) {
        period_ = period > MIN_PERIOD ? period : MIN_PERIOD;
        retune();
    }
 
    /// Excitation point along the string, 0 at end a, 1 at end b (a string: nut and bridge).
    void set_junction(float beta) {
        // A junction at an end is not a junction; NaN lands on the edge.
        beta_ = math::clamp(beta, JUNCTION_EDGE, 1.0f - JUNCTION_EDGE);
        retune();
    }
 
    /// Loss filter cutoff: how fast the high partials give up. Signal units.
    void set_damping_su(float su) {
        loss_.set_cutoff_su(math::clamp01(su));   // NaN → darkest
    }
 
    /// The wave impedance Z the junction reports and converts forces with.
    void set_impedance(float z) { impedance_ = z > 0.0f ? z : DEFAULT_IMPEDANCE; }
    float impedance() const { return impedance_; }
 
    /// The junction's displacement right now (wave units × samples).
    float displacement() const { return x_; }
 
    /**
     * @brief Round-trip gain, the passivity margin
     *
     * Below 1 the loop decays; at 1 it is lossless and only the damping
     * filter removes energy. Clamped at 1 so the loop can never generate.
     */
    void set_loop_gain(float gain) {
        gain_ = math::clamp01(gain);   // NaN → silent
    }
 
    /**
     * @brief Advance one sample
     *
     * @param excite  an Exciter (excite(const JunctionState&) → force), or a
     *                callable of the arriving velocity returning the wave to
     *                inject
     * @param end_b   called with the wave arriving at end b (after the loop's
     *                loss and allpass); returns the reflection
     * @return the wave arriving at end b — the pickup
     */
    template<typename Excite, typename EndB>
    float tick(Excite&& excite, EndB&& end_b) {
        static_assert(is_wave_callable_v<EndB>, "the end takes a wave and returns one: float(float)");
        const float to_b = end_b_.read(d_b_);
        const float to_a = end_a_.read(d_a_);
 
        // End a is a hard reflection. Leg b carries the loss and the
        // fractional tuning, then whatever end b sends back.
        const float from_a = -to_a;
        const float arriving = loss_.process(allpass_.process(to_b));
        const float from_b = gain_ * end_b(arriving);
 
        const float v_free = from_a + from_b;
        const float injection = inject(static_cast<Excite&&>(excite), v_free);
 
        end_a_.write(from_b + injection);
        end_b_.write(from_a + injection);
        x_ = leak_ * x_ + (v_free + injection);   // the point moved by its velocity this sample
        return to_b;
    }
 
    /// Advance one sample against a rigid end b: reflection −1.
    template<typename Excite>
    float tick(Excite&& excite) {
        return tick(static_cast<Excite&&>(excite), rigid_end);
    }
 
    /// Advance one sample with nothing injected: the string rings on alone.
    float tick() {
        return tick([](float) { return 0.0f; }, rigid_end);
    }
 
    /// An end with nothing on it: a wall.
    static float rigid_end(float arriving) { return -arriving; }
 
    /// What the exciter sees this sample, for a scope.
    JunctionState junction_state(float v_free) const { return JunctionState{v_free, x_, impedance_}; }
 
    float period_samples() const { return period_; }
    float junction() const { return beta_; }
 
    /// Segment lengths actually read — the split the junction makes, clamped to the buffers.
    float delay_a_samples() const { return d_a_; }
    float delay_b_samples() const { return d_b_; }
 
    /// The allpass's share of the period: its nominal sample plus the rounding residual.
    float allpass_delay_samples() const { return allpass_.group_delay_samples(); }
 
    /**
     * @brief Energy in the two segments, in wave units: the sum of squares of
     *        every sample in flight (physical-modeling.md §5.1, §8.3)
     *
     * The loss filter's and allpass's one-sample states are left out; they
     * are a sample each against hundreds in the lines. Never called from
     * audio code.
     */
    float energy() const {
        float e = 0.0f;
        const uint32_t n_a = static_cast<uint32_t>(d_a_);
        const uint32_t n_b = static_cast<uint32_t>(d_b_);
        for (uint32_t i = 1; i <= n_a; ++i) { const float v = end_a_.read(static_cast<float>(i)); e += v * v; }
        for (uint32_t i = 1; i <= n_b; ++i) { const float v = end_b_.read(static_cast<float>(i)); e += v * v; }
        return e;
    }
 
    void reset() {
        end_a_.reset();
        end_b_.reset();
        loss_.reset();
        allpass_.reset();
        x_ = 0.0f;
        leak_ = 1.0f - math::TWO_PI * DISPLACEMENT_LEAK_HZ / types::SAMPLE_RATE_F;   // rate-derived, at reset
    }
 
    /**
     * @brief The loop's wiring as data (AP-037). Never called from audio code.
     *
     * @param external_end_b  true when an owner supplies the bridge callable:
     *        the wall is left out and the named ports "to end b" and "from
     *        end b" are the two ends the owner's termination connects to. `in` is the excitation
     *        junction; `out` is the pickup (the wave arriving at the bridge).
     */
    diagram::Ports describe(diagram::Graph& g, uint8_t parent = diagram::NO_GROUP,
                            bool external_end_b = false) const {
        using diagram::Kind;
        using diagram::Wave;
        diagram::Ports p;
        p.group = g.add_group("StringLoop", parent);
 
        const uint8_t junction = g.add_node(Kind::Junction, "excitation", p.group);
        const uint8_t seg_a = g.add_node(Kind::DelayLine, "segment a", p.group, d_a_);
        const uint8_t wall_a = g.add_node(Kind::Reflection, "end a", p.group, -1.0f);
        const uint8_t seg_b = g.add_node(Kind::DelayLine, "segment b", p.group, d_b_);
        const uint8_t allpass = g.add_node(Kind::Allpass, "tuning", p.group,
                                           allpass_.group_delay_samples());
        const uint8_t loss = g.add_node(Kind::Loss, "loss", p.group, loss_cutoff_hz());
 
        g.add_edge(junction, seg_a, Wave::Velocity);
        g.add_edge(seg_a, wall_a, Wave::Velocity);
        g.add_edge(wall_a, junction, Wave::Velocity);
        g.add_edge(junction, seg_b, Wave::Velocity);
        g.add_edge(seg_b, allpass, Wave::Velocity);
        g.add_edge(allpass, loss, Wave::Velocity);
 
        if (external_end_b) {
            p.add("to end b", loss);
            p.add("from end b", junction);
        } else {
            const uint8_t wall = g.add_node(Kind::Reflection, "end b", p.group, -1.0f);
            g.add_edge(loss, wall, Wave::Velocity);
            g.add_edge(wall, junction, Wave::Velocity);
        }
        p.in = junction;
        p.out = seg_b;   // the pickup reads the wave arriving at the bridge
        return p;
    }
 
private:
    using DelayLine = sbl::dsp::primitives::DelayLine;
    using OnePole = sbl::dsp::primitives::OnePole;
    using Allpass = sbl::dsp::primitives::ThiranAllpass<1>;
 
    static constexpr float DAMPING_MIN_HZ = 500.0f;
    static constexpr float DAMPING_MAX_HZ = 16000.0f;
    static constexpr float DEFAULT_DAMPING_SU = 0.8f;
    static constexpr float DEFAULT_PERIOD = 100.0f;
    static constexpr float DEFAULT_JUNCTION = 0.3f;
    static constexpr float DEFAULT_LOOP_GAIN = 0.995f;
    static constexpr float JUNCTION_EDGE = 0.02f;
    static constexpr float ALLPASS_NOMINAL = 1.0f;   ///< Thiran<1> sits at 1 sample
    static constexpr float MIN_SEGMENT = 1.0f;       ///< a delay line reads at least one sample back
    static constexpr float HALF_SAMPLE = 0.5f;       ///< the allpass tunes within this of nominal
    static constexpr float TWO = 2.0f;
 
    /// The loss filter's cutoff back in Hz from its coefficient: c = 1 − e^(−2πf/fs).
    float loss_cutoff_hz() const {
        const float c = loss_.coefficient();
        if (c <= 0.0f || c >= 1.0f) return 0.0f;
        return -std::log(1.0f - c) * sbl::dsp::types::SAMPLE_RATE_F / math::TWO_PI;
    }
 
    /// The exciter's force becomes the wave the loop injects, F / 2Z; a bare
    /// callable of the velocity returns the wave itself.
    template<typename Excite>
    float inject(Excite&& excite, float v_free) {
        if constexpr (is_exciter_v<std::decay_t<Excite>>) {
            return excite.excite(JunctionState{v_free, x_, impedance_}) / (TWO * impedance_);
        } else {
            static_assert(is_wave_callable_v<Excite>, "excite takes the junction state or the velocity");
            return excite(v_free);
        }
    }
 
    static float round_to_sample(float x) {
        return static_cast<float>(static_cast<int32_t>(x + HALF_SAMPLE));
    }
 
    static float clamp_segment(float d, float limit) {
        if (d < MIN_SEGMENT) return MIN_SEGMENT;
        if (d > limit) return limit;
        return d;
    }
 
    /**
     * @brief Split the period between the segments; the remainder tunes the allpass
     *
     * The nut segment is rounded and clamped first, the bridge segment takes
     * what is left and is rounded and clamped in turn, and whatever fraction
     * remains goes to the allpass — so a clamp on either side comes out of
     * the other, not out of the total. The allpass can only absorb half a
     * sample either way; past that the bridge segment moves by one.
     */
    void retune() {
        const float total = period_ - ALLPASS_NOMINAL;
 
        float a = clamp_segment(round_to_sample(total * beta_), limit_a_);
        float b = clamp_segment(round_to_sample(total - a), limit_b_);
 
        // Whichever segment can still move takes the whole sample; segment
        // b first, segment a when b is pinned at an end.
        float residual = total - a - b;
        if (residual > HALF_SAMPLE) {
            if (b < limit_b_)         { b += 1.0f; residual -= 1.0f; }
            else if (a < limit_a_)    { a += 1.0f; residual -= 1.0f; }
        } else if (residual < -HALF_SAMPLE) {
            if (b > MIN_SEGMENT)      { b -= 1.0f; residual += 1.0f; }
            else if (a > MIN_SEGMENT) { a -= 1.0f; residual += 1.0f; }
        }
 
        d_a_ = a;
        d_b_ = b;
        allpass_.set_delay(ALLPASS_NOMINAL + residual);   // clamps to ±half a sample
    }
 
    DelayLine end_a_;
    DelayLine end_b_;
    OnePole loss_;
    Allpass allpass_;
 
    float period_ = DEFAULT_PERIOD;
    float beta_ = DEFAULT_JUNCTION;
    float gain_ = DEFAULT_LOOP_GAIN;
    float limit_a_ = 1.0f;
    float limit_b_ = 1.0f;
    float d_a_ = MIN_SEGMENT;
    float d_b_ = MIN_SEGMENT;
    float impedance_ = DEFAULT_IMPEDANCE;
    float x_ = 0.0f;       ///< junction displacement
    float leak_ = 1.0f;    ///< set at reset from the sample rate
};
 
} // namespace sbl::dsp::pm
 
#endif // SBL_DSP_PM_STRING_LOOP_HPP_