mass_spring_damper.hpp

Go to the documentation of this file.

// sbl/dsp/pm/mass_spring_damper.hpp — A lumped resonator that can be a load (Physical modeling — primitive)
//
//     m x'' + c x' + k x = F
//
// One mass on a spring with a damper: the smallest thing that resonates and
// the smallest thing a string can terminate on. Integrated with the
// trapezoidal rule, which is passive — the discrete system cannot make
// energy, only store or lose it (RPT-029 Part 4 rule 1) — and which gives
// the load a form a junction can solve against without a delay-free loop:
// at each sample the velocity is an instantaneous part plus a history part,
//
//     v[n] = g · F[n] + h[n]
//
// where g is the instantaneous admittance (set by the parameters) and h is
// what the past state contributes.
//
// Two ways to say what the load is. The physical one — frequency, Q, mass —
// is exact and awkward: mass in string-impedance units means nothing at a
// panel, and a light mass is a free end at every musical frequency (AP-036,
// 2026-09-13). The musical one — frequency, Q, match — names what a string
// feels: where the resonance is, how wide it is, and how much of the string
// leaves into it at resonance. Match is the damper relative to the string's
// impedance, c/Z; at resonance the load is that damper alone and the
// reflection is (1 − match)/(1 + match) — 1 absorbs everything, 0.1 keeps
// 82 % of the wave, 10 is nearly a wall again. Mass follows: m = match·Q/ω0.
// Away from resonance the load then looks rigid where match·Q·(ω/ω0) ≫ 1. A caller that knows F outright uses
// commit(); a junction that has to solve for F first reads admittance()
// and history(), solves, then commit()s the force it found.
//
// Units are the waveguide's: time in samples, velocity in wave units, force
// in wave units times the string's wave impedance, which the loop
// normalises to 1. A mass of 1 is the mass whose impedance at one radian per
// sample (7.6 kHz at 48 kHz) equals the string's — which makes it a free end
// at any musical frequency. Rigid needs m·ω ≫ 1 at the lowest note played:
// m ≫ 30 for a string at 250 Hz, m ≫ 230 down at C1. Between resonance and
// rigid the load's impedance is stiffness-like below ω0 (k/ω) and mass-like
// above (m·ω); at ω0 it is the damper alone, c = m·ω0/Q, and that is where
// the string leaks — most when c matches the string's 1. A light, high-Q
// load is a wolf, not a body (AP-036, 2026-09-13).
//
// Usage:
//   sbl::dsp::pm::MassSpringDamper body;
//   body.set_frequency(220.0f);
//   body.set_q(12.0f);
//   body.set_mass(4.0f);
//   float v = body.commit(force);             // as a resonator
//   // or, from a junction:
//   const float g = body.admittance(), h = body.history();
//   ... solve for F ...
//   float v = body.commit(F);
 
#ifndef SBL_DSP_PM_MASS_SPRING_DAMPER_HPP_
#define SBL_DSP_PM_MASS_SPRING_DAMPER_HPP_
 
#include <sbl/dsp/diagram/graph.hpp>
#include <sbl/dsp/math/clamp.hpp>
#include <sbl/dsp/pm/contracts.hpp>
#include <sbl/dsp/math/constants.hpp>
#include <sbl/dsp/types/fixed.hpp>
 
namespace sbl::dsp::pm {
 
class MassSpringDamper {
public:
    /// @note All public methods are ISR-safe — bounded computation, no I/O.
 
    static constexpr float MIN_FREQUENCY_HZ = 1.0f;
    static constexpr float MAX_FREQUENCY_FRACTION = 0.45f;   ///< of the sample rate
    static constexpr float MIN_Q = 0.05f;
    static constexpr float MAX_Q = 1000.0f;
    static constexpr float MIN_MASS = 1e-4f;
    static constexpr float MAX_MASS = 1e6f;
    static constexpr float MIN_MATCH = 1e-3f;
    static constexpr float MAX_MATCH = 100.0f;
 
    MassSpringDamper() { update(); }
 
    // ─── Parameters (natural units) ──────────────────────────────────
 
    /// Undamped resonance, Hz.
    void set_frequency(float hz) {
        const float max_hz = types::SAMPLE_RATE_F * MAX_FREQUENCY_FRACTION;
        frequency_hz_ = math::clamp(hz, MIN_FREQUENCY_HZ, max_hz);   // NaN → min
        update();
    }
 
    /// Quality factor: how many cycles the ringing lasts.
    void set_q(float q) {
        q_ = math::clamp(q, MIN_Q, MAX_Q);   // NaN → min
        update();
    }
 
    /// Mass in string-impedance units; heavier is closer to rigid.
    void set_mass(float m) {
        mass_ = math::clamp(m, MIN_MASS, MAX_MASS);   // NaN → min
        update();
    }
 
    /**
     * @brief The musical parametrisation: where, how wide, how much it takes
     *
     * @param hz     resonance
     * @param q      how long it rings — also how rigid it looks away from
     *               resonance, since the mass follows from q and match
     * @param match  damper relative to the string's impedance at resonance;
     *               1 absorbs the resonant wave completely
     */
    void set_resonance(float hz, float q, float match) {
        const float m = math::clamp(match, MIN_MATCH, MAX_MATCH);   // NaN → min
        set_frequency(hz);
        set_q(q);
        const float w0 = math::TWO_PI * frequency_hz_ / types::SAMPLE_RATE_F;
        set_mass(m * q_ / w0);
    }
 
    float frequency() const { return frequency_hz_; }
    float q() const { return q_; }
    float mass() const { return mass_; }
 
    /// The damper relative to the string's impedance: c / Z with Z = 1.
    float match() const { return c_; }
 
    // ─── The junction interface ──────────────────────────────────────
 
    /// Instantaneous admittance g: how much of this sample's force shows up in this sample's velocity.
    float admittance() const { return g_; }
 
    /// What the past state contributes to this sample's velocity, before any
    /// force. Fixed once the previous sample is committed, so it is computed
    /// there and read here — the junction reads it, then commits.
    float history() const { return h_; }
 
    /// Apply this sample's force and advance: returns the velocity.
    float commit(float force) {
        const float v = g_ * force + h_;
        x_ += HALF * (v + v_);
        v_ = v;
        a_ = (force - c_ * v - k_ * x_) * inv_mass_;
        update_history();
        return v;
    }
 
    float velocity() const { return v_; }
    float position() const { return x_; }
 
    /// Stored energy ½mv² + ½kx², for passivity tests.
    float energy() const { return HALF * (mass_ * v_ * v_ + k_ * x_ * x_); }
 
    void reset() {
        x_ = 0.0f;
        v_ = 0.0f;
        a_ = 0.0f;
        update();   // rate-derived coefficients belong here (FDP-055 Addendum A); h_ follows
    }
 
    /// One Load node: frequency and Q. Never called from audio code (AP-037).
    diagram::Ports describe(diagram::Graph& g, uint8_t group, const char* name = "mode") const {
        diagram::Ports p;
        p.in = p.out = g.add_node(diagram::Kind::Load, name, group, frequency_hz_, q_);
        p.group = group;
        return p;
    }
 
private:
    static constexpr float HALF = 0.5f;
    static constexpr float QUARTER = 0.25f;
 
    /// Coefficients from the parameters; time step is one sample.
    void update() {
        const float w0 = math::TWO_PI * frequency_hz_ / types::SAMPLE_RATE_F;   // rad/sample
        k_ = mass_ * w0 * w0;
        c_ = mass_ * w0 / q_;
        inv_mass_ = 1.0f / mass_;
        // Trapezoidal: v[n] (m + c/2 + k/4) = F/2 + history terms
        const float d = mass_ + HALF * c_ + QUARTER * k_;
        inv_d_ = 1.0f / d;
        g_ = HALF * inv_d_;
        update_history();
    }
 
    void update_history() {
        h_ = (mass_ * v_ - QUARTER * k_ * v_ - HALF * k_ * x_ + HALF * mass_ * a_) * inv_d_;
    }
 
    float frequency_hz_ = 220.0f;
    float q_ = 10.0f;
    float mass_ = 1.0f;
 
    float k_ = 0.0f;
    float c_ = 0.0f;
    float inv_mass_ = 1.0f;
    float inv_d_ = 1.0f;
    float g_ = 0.0f;
 
    float x_ = 0.0f;   ///< position
    float v_ = 0.0f;   ///< velocity
    float a_ = 0.0f;   ///< acceleration at the last sample
    float h_ = 0.0f;   ///< history for the next sample, from the state above
};
 
static_assert(is_load_v<MassSpringDamper> && is_described_v<MassSpringDamper>, "MassSpringDamper is a Load and describes itself (physical-modeling.md §4.1, §7)");
 
} // namespace sbl::dsp::pm
 
#endif // SBL_DSP_PM_MASS_SPRING_DAMPER_HPP_