table_friction.hpp

Go to the documentation of this file.

// sbl/dsp/pm/table_friction.hpp — A bow: friction solved against the string (Physical modeling — primitive)
//
// Friction at a bow depends on the relative velocity between bow and string,
// but the string's velocity depends on the force applied. The two are one
// equation:
//
//     dv = dv_free - F / (2Z)      dv_free = v_bow - v_incoming
//     F  = f(dv)
//
// f(dv) is the friction characteristic; the string constrains the operating
// point to a line of slope -2Z through (dv_free, 0) — the load line. The
// answer is where they cross. McIntyre, Schumacher and Woodhouse (1983)
// showed the crossing is a function of dv_free alone, so lute solves it
// offline and this reads the answer: one lookup, no iteration. STK's bow
// table (Smith, Cook) and Elements' Resonator::BowTable (Emilie Gillet, MIT)
// are the same object, with the impedance folded into constants.
//
// Units are the waveguide's wave variables, normalised so 2Z = 1. Holding
// the string stuck then costs exactly F = dv_free, which makes the stick
// test one comparison.
//
// This is the "friction curve model": force a function of sliding speed
// alone. Woodhouse & Galluzzo (Acta Acustica 90, 2004, §3.5) call it the
// one part of the bowed-string model that measurements contradict — real
// rosin traces a hysteresis loop in the force–velocity plane, and the
// thermal plastic model is the candidate replacement. It is the right first
// layer (every classic result rests on it) and a known simplification, not a
// finished bow. RPT-029: provide both, architect around neither.
//
// ## Stick and slip
//
// Stuck, the string travels with the bow and the contact carries whatever
// force that takes — until it exceeds what static friction can supply, and
// the string releases. Sliding, the force follows the solved branch until
// the relative velocity collapses and the bow recaptures the string. The
// release and capture points differ: that hysteresis is stick-slip, and it
// is what gives a bowed attack its bite. One bool of state carries it.
//
// ## Energy
//
// A bow is a source: it does work on the string, so this has no passivity
// proof and is one of the components RPT-029 Part 4 rule 6 names — the
// SoftLimiter on the injection and the NaN guard downstream are
// load-bearing. Force never exceeds dv_free (the table is generated with
// that property asserted), so the bow cannot drive the string faster than
// itself.
//
// A law, not an exciter: a Bow<TableFriction> owns the speed, the weight and
// the junction; this owns the friction. In wave units: the force returned is
// F / 2Z, what a string of impedance Z takes from it.
//
// Usage:
//   sbl::dsp::pm::TableFriction law;
//   law.set_bow(lut::bow_stribeck_128, lut::bow_stribeck_limits, 8, 128, 129);
//   law.set_force_su(pressure);
//   float injection = law.process(v_bow - v_incoming);
 
#ifndef SBL_DSP_PM_TABLE_FRICTION_HPP_
#define SBL_DSP_PM_TABLE_FRICTION_HPP_
 
#include <cmath>
#include <cstdint>
 
#include <sbl/dsp/diagram/graph.hpp>
#include <sbl/dsp/pm/contracts.hpp>
 
namespace sbl::dsp::pm {
 
/**
 * @brief One bow: a lute-solved family and the shape of its arrays
 *
 * The generated headers give a family and its limits; this carries them
 * together with the geometry needed to read them, so a bow can be passed
 * around as one value and a set of bows is an array.
 */
struct BowTable {
    const float* family = nullptr;   ///< force_slices rows, light bow first
    const float* limits = nullptr;   ///< static-friction limit per slice
    uint8_t slices = 0;
    uint16_t size = 0;               ///< samples per row
    uint16_t stride = 0;             ///< row pitch (size + guard points)
    float dv_max = 1.0f;             ///< relative velocity at the last sample
};
 
class TableFriction {
public:
    /// @note All public methods are ISR-safe — bounded computation, no I/O.
 
    /**
     * @brief Install a bow: a lute-solved family and its per-slice limits
     *
     * @param family  force_slices rows of `stride` floats, light bow first
     * @param limits  one static-friction limit per slice
     * @param slices  number of force slices
     * @param size    samples per row spanning dv_free in [0, dv_max]
     * @param stride  row pitch (size + guard points)
     * @param dv_max  relative velocity the last sample corresponds to
     */
    void set_bow(const float* family, const float* limits, uint8_t slices,
                 uint16_t size, uint16_t stride, float dv_max = 1.0f) {
        family_ = family;
        limits_ = limits;
        slices_ = slices;
        size_ = size;
        stride_ = stride;
        dv_max_ = dv_max > 0.0f ? dv_max : 1.0f;
        inv_dv_max_ = 1.0f / dv_max_;   // the per-sample lookup multiplies, never divides
        set_force_su(force_su_);
    }
 
    /// Install a bow described in one value.
    void set_bow(const BowTable& bow) {
        set_bow(bow.family, bow.limits, bow.slices, bow.size, bow.stride, bow.dv_max);
    }
 
    /// Bow force across the family [0, 1]; crossfades adjacent slices.
    void set_force_su(float su) {
        force_su_ = su < 0.0f ? 0.0f : (su > 1.0f ? 1.0f : su);
        if (family_ == nullptr || slices_ == 0) return;
 
        const float pos = force_su_ * static_cast<float>(slices_ - 1);
        slice_ = static_cast<uint8_t>(pos);
        if (slice_ >= slices_ - 1 && slices_ > 1) slice_ = static_cast<uint8_t>(slices_ - 2);
        blend_ = slices_ > 1 ? pos - static_cast<float>(slice_) : 0.0f;
 
        limit_ = slices_ > 1
            ? limits_[slice_] + blend_ * (limits_[slice_ + 1] - limits_[slice_])
            : limits_[0];
    }
 
    /**
     * @brief One sample: relative velocity in, injection force out
     *
     * @param dv_free v_bow - v_incoming, in wave variables
     */
    float process(float dv_free) {
        if (family_ == nullptr) return 0.0f;
 
        const float a = std::fabs(dv_free);
 
        if (stuck_) {
            if (a <= limit_) return dv_free;   // travelling with the bow
            stuck_ = false;                    // static friction gave way
        }
 
        const float f = lookup(a);
        if (a - f <= CAPTURE_VELOCITY) stuck_ = true;   // the bow caught it again
        return dv_free < 0.0f ? -f : f;
    }
 
    /// True while the string travels with the bow.
    bool sticking() const { return stuck_; }
 
    /// Static-friction limit at the current bow force — where the string releases.
    float limit() const { return limit_; }
 
    void reset() { stuck_ = true; }
 
    /// One Junction node: the friction nonlinearity. Never called from audio code (AP-037).
    diagram::Ports describe(diagram::Graph& g, uint8_t group, const char* name = "bow friction") const {
        diagram::Ports p;
        p.in = p.out = g.add_node(diagram::Kind::Junction, name, group);
        p.group = group;
        return p;
    }
 
private:
    /// Relative velocity below which the bow is considered to have recaptured.
    static constexpr float CAPTURE_VELOCITY = 1e-4f;
 
    /// Solved force at |dv_free|, interpolated within and across slices.
    float lookup(float a) const {
        const float last = static_cast<float>(size_ - 1);
        float x = a * inv_dv_max_ * last;
        // NaN fails every comparison, so test the positive case: a NaN
        // velocity reads the table's origin rather than casting NaN, which
        // is undefined.
        x = (x > 0.0f) ? (x < last ? x : last) : 0.0f;
 
        const uint16_t i = static_cast<uint16_t>(x);
        const float frac = x - static_cast<float>(i);
        const uint16_t j = (i + 1 < size_) ? static_cast<uint16_t>(i + 1) : i;
 
        const float* row0 = family_ + static_cast<uint32_t>(slice_) * stride_;
        const float v0 = row0[i] + frac * (row0[j] - row0[i]);
        if (slices_ < 2) return v0;
 
        const float* row1 = row0 + stride_;
        const float v1 = row1[i] + frac * (row1[j] - row1[i]);
        return v0 + blend_ * (v1 - v0);
    }
 
    const float* family_ = nullptr;
    const float* limits_ = nullptr;
    uint8_t slices_ = 0;
    uint16_t size_ = 0;
    uint16_t stride_ = 0;
    float dv_max_ = 1.0f;
    float inv_dv_max_ = 1.0f;
 
    float force_su_ = 1.0f;
    uint8_t slice_ = 0;
    float blend_ = 0.0f;
    float limit_ = 0.0f;
    bool stuck_ = true;
};
 
static_assert(is_friction_law_v<TableFriction>, "TableFriction is a friction law (physical-modeling.md §4.5)");
 
} // namespace sbl::dsp::pm
 
#endif // SBL_DSP_PM_TABLE_FRICTION_HPP_