cyclic_stretch.hpp

Go to the documentation of this file.

#pragma once

 
#include <pulp/signal/crossfade.hpp>
#include <pulp/signal/denormal.hpp>
#include <pulp/signal/units.hpp>
 
#include <algorithm>
#include <cmath>
#include <cstddef>
#include <cstdint>
#include <vector>
 
namespace pulp::signal {

inline long long round_half_even(double x) noexcept {
    const double floored = std::floor(x);
    const double fraction = x - floored;
    const auto below = static_cast<long long>(floored);
    if (fraction > 0.5) return below + 1;
    if (fraction < 0.5) return below;
    return (below & 1LL) == 0 ? below : below + 1;  // tie → the even neighbour
}

struct CyclicStretchRegime {
    double cycle_hz = 200.0;
    int grain_periods = 1;
    double crossfade_pct = 25.0;
    double crossfade_shape = 1.0;
};

inline constexpr CyclicStretchRegime kCyclicStretchShortFrame{200.0, 1, 45.0, 1.0};

inline constexpr CyclicStretchRegime kCyclicStretchLongFrame{200.0, 8, 8.0, 0.0};

template <typename SampleType = float>
class CyclicStretchT {
public:
    // ── Design parameters (the complete roster) ───────────────────────────

    static constexpr double kCycleHzMin = 20.0;
    static constexpr double kCycleHzMax = 2000.0;
    static constexpr double kCycleHzDefault = 200.0;

    static constexpr int kGrainPeriodsMin = 1;
    static constexpr int kGrainPeriodsMax = 16;
    static constexpr int kGrainPeriodsDefault = 1;

    static constexpr double kCrossfadePctMin = 1.0;
    static constexpr double kCrossfadePctMax = 50.0;
    static constexpr double kCrossfadePctDefault = 25.0;

    static constexpr double kCrossfadeShapeDefault = 1.0;

    static constexpr double kStretchRatioMin = 0.25;
    static constexpr double kStretchRatioMax = 4.0;
    static constexpr double kStretchRatioDefault = 1.0;

    static constexpr double kCaptureMsMin = 250.0;
    static constexpr double kCaptureMsMax = 8000.0;
    static constexpr double kCaptureMsDefault = 2000.0;

    static constexpr double kMixDefault = 100.0;

    static constexpr double kOutputDbMin = -24.0;
    static constexpr double kOutputDbMax = 12.0;
    static constexpr double kOutputDbDefault = 0.0;

    static constexpr int kMinCycleSamples = 4;

    static constexpr int kMinHopSamples = 2;

    static constexpr double kWorstCaseGain = 1.41421356237309504880;
 
    // ── Lifecycle ─────────────────────────────────────────────────────────

    void prepare(double sample_rate) {
        if (std::isfinite(sample_rate) && sample_rate > 0.0) sample_rate_ = sample_rate;
 
        // Sized for the extremes of the declared ranges, not for the current
        // settings: `set_capture_ms` and `set_cycle_hz` must never reallocate,
        // and the only way to promise that is to buy the worst case once.
        const auto ring_needed = static_cast<std::size_t>(
            std::ceil(units::ms_to_samples(kCaptureMsMax, sample_rate_))) + 2u;
        ring_mask_ = power_of_two_mask(ring_needed);
        ring_.assign(ring_mask_ + 1u, SampleType{0});
 
        const auto longest_cycle = static_cast<long long>(
            std::llround(sample_rate_ / kCycleHzMin));
        max_grain_ = std::max<long long>(
            1, longest_cycle * kGrainPeriodsMax);
        // A grain writes `N` samples starting at the current write cursor, and
        // that cursor never runs more than `N` ahead of the read cursor (grains
        // are produced only as far as the current output sample needs). So the
        // live span is `N`, and twice that leaves every slot a full lap between
        // being zeroed and being written again.
        accum_mask_ = power_of_two_mask(static_cast<std::size_t>(2 * max_grain_) + 2u);
        accum_.assign(accum_mask_ + 1u, 0.0);
        accum_epoch_.assign(accum_mask_ + 1u, 0u);
 
        update_lengths();
        reset();
    }

    void reset() {
        std::fill(ring_.begin(), ring_.end(), SampleType{0});
        std::fill(accum_.begin(), accum_.end(), 0.0);
        std::fill(accum_epoch_.begin(), accum_epoch_.end(), 0u);
        accum_generation_ = 1u;
        reset_runtime_state();
    }
 
private:
    void audio_fault_reset() noexcept {
        ++accum_generation_;
        if (accum_generation_ == 0u) accum_generation_ = 1u;
        reset_runtime_state();
    }
 
    void reset_runtime_state() noexcept {
        write_index_ = 0;
        total_written_ = 0;
        out_read_ = 0;
        grain_ = 0;
        grain_out_pos_ = 0;
        read_pos_ = 0;
        prev_out_pos_ = 0;
        previous_fade_out_ = crossfade_length_;
        end_previous_grain_ = 0;
        end_two_grains_ago_ = 0;
    }
 
public:
 
    // ── Parameters ────────────────────────────────────────────────────────
    //
    // None of these resets the schedule or the ring. That is deliberate: a live
    // automation move on `cycle_hz` should slide the buzz, not restart the
    // machine with a click. `reset()` is the only full clear.
 
    void set_cycle_hz(double hz) {
        if (!std::isfinite(hz)) return;
        cycle_hz_ = std::clamp(hz, kCycleHzMin, kCycleHzMax);
        update_lengths();
    }
 
    void set_grain_periods(int periods) {
        grain_periods_ = std::clamp(periods, kGrainPeriodsMin, kGrainPeriodsMax);
        update_lengths();
    }
 
    void set_crossfade_pct(double pct) {
        if (!std::isfinite(pct)) return;
        crossfade_pct_ = std::clamp(pct, kCrossfadePctMin, kCrossfadePctMax);
        update_lengths();
    }
 
    void set_crossfade_shape(double shape) noexcept {
        if (!std::isfinite(shape)) return;
        crossfade_shape_ = std::clamp(shape, 0.0, 1.0);
    }
 
    void set_stretch_ratio(double ratio) noexcept {
        if (!std::isfinite(ratio)) return;
        stretch_ratio_ = std::clamp(ratio, kStretchRatioMin, kStretchRatioMax);
    }
 
    void set_capture_ms(double ms) {
        if (!std::isfinite(ms)) return;
        capture_ms_ = std::clamp(ms, kCaptureMsMin, kCaptureMsMax);
        update_lengths();
    }
 
    void set_mix(double percent) noexcept {
        if (!std::isfinite(percent)) return;
        mix_ = std::clamp(percent, 0.0, 100.0) * 0.01;
    }
 
    void set_output_db(double db) noexcept {
        if (!std::isfinite(db)) return;
        output_gain_ = units::db_to_linear(std::clamp(db, kOutputDbMin, kOutputDbMax));
    }

    void set_regime(const CyclicStretchRegime& regime) {
        set_cycle_hz(regime.cycle_hz);
        set_grain_periods(regime.grain_periods);
        set_crossfade_pct(regime.crossfade_pct);
        set_crossfade_shape(regime.crossfade_shape);
    }
 
    // ── Resolved lengths ──────────────────────────────────────────────────

    int cycle_samples() const noexcept { return static_cast<int>(cycle_length_); }
    int grain_samples() const noexcept { return static_cast<int>(grain_length_); }
    int crossfade_samples() const noexcept { return static_cast<int>(crossfade_length_); }
    int hop_samples() const noexcept { return static_cast<int>(hop_length_); }
    double flutter_hz() const noexcept {
        return hop_length_ > 0 ? sample_rate_ / static_cast<double>(hop_length_) : 0.0;
    }

    long long capture_window_samples() const noexcept { return capture_window_; }
 
    double sample_rate() const noexcept { return sample_rate_; }
    double stretch_ratio() const noexcept { return stretch_ratio_; }
    double crossfade_shape() const noexcept { return crossfade_shape_; }

    int latency_samples() const noexcept { return static_cast<int>(grain_length_); }

    double worst_case_gain() const noexcept {
        return (1.0 - crossfade_shape_) + crossfade_shape_ * kWorstCaseGain;
    }

    long long schedule_input_position(long long out_pos) const noexcept {
        const double ideal = static_cast<double>(out_pos) / stretch_ratio_;
        return cycle_length_ * round_half_even(ideal / static_cast<double>(cycle_length_));
    }

    long long grain_input_position(long long g) const noexcept {
        return schedule_input_position(g * hop_length_);
    }

    long long grain_count() const noexcept { return grain_; }

    long long next_grain_out_pos() const noexcept { return grain_out_pos_; }

    long long read_position() const noexcept { return read_pos_; }

    long long total_captured() const noexcept { return total_written_; }
 
    // ── Audio ─────────────────────────────────────────────────────────────

    void process(const SampleType* in, SampleType* out, int n) {
        if (n <= 0) return;
        if (ring_.empty() || accum_.empty() || in == nullptr) {
            std::fill_n(out, n, SampleType{0});
            return;
        }
        for (int i = 0; i < n; ++i) {
            if (!std::isfinite(static_cast<double>(in[i]))) {
                audio_fault_reset();
                out[i] = SampleType{0};
                continue;
            }
            ring_[write_index_] = in[i];
            write_index_ = (write_index_ + 1u) & ring_mask_;
            ++total_written_;
 
            // The initial fill. Grain 0 cannot exist until its N source samples
            // do, and that wait IS the reported latency.
            if (total_written_ <= grain_length_) {
                out[i] = SampleType{0};
                continue;
            }
 
            // Produce exactly the grains the current output sample depends on.
            // Output sample p is final once every grain starting at or before p
            // has been written, because a later grain starts later than p by
            // construction.
            while (grain_out_pos_ <= out_read_) emit_grain();
 
            const auto slot = static_cast<std::size_t>(out_read_) & accum_mask_;
            const double wet = accum_epoch_[slot] == accum_generation_
                                   ? snap_to_zero(accum_[slot])
                                   : 0.0;
            accum_[slot] = 0.0;
            accum_epoch_[slot] = accum_generation_;
            ++out_read_;
 
            // The dry path is delayed by exactly N, so `mix` blends two things
            // that are talking about the same moment. The `+ 1` is load-bearing
            // and was measured, not reasoned: `write_index_` has ALREADY been
            // advanced past the sample just captured, so the newest sample sits
            // at `write_index_ - 1`, and reading `write_index_ - N` delays by
            // N−1. That off-by-one is invisible on a tone and shows up as a
            // one-sample pre-echo on an impulse at `mix` between 0 and 1.
            const auto dry_index = static_cast<std::size_t>(
                (write_index_ + ring_mask_ + 1u -
                 static_cast<std::size_t>(grain_length_ + 1)) & ring_mask_);
            const double dry = static_cast<double>(ring_[dry_index]);
 
            const double result = output_gain_ * (mix_ * wet + (1.0 - mix_) * dry);
            if (!std::isfinite(result)) {
                audio_fault_reset();
                out[i] = SampleType{0};
            } else {
                out[i] = static_cast<SampleType>(result);
            }
        }
    }
 
private:
    static std::size_t power_of_two_mask(std::size_t at_least) {
        std::size_t size = 2;
        while (size < at_least) size <<= 1;
        return size - 1u;
    }

    void advance_read_cursor(long long out_pos) noexcept {
        read_pos_ += schedule_input_position(out_pos) -
                     schedule_input_position(prev_out_pos_);
        prev_out_pos_ = out_pos;
 
        const long long newest = total_written_ - grain_length_;
        const long long oldest = std::max<long long>(0, total_written_ - capture_window_);
        const long long span = std::max<long long>(cycle_length_, newest - oldest);
 
        if (read_pos_ > newest) {
            const long long laps = (read_pos_ - newest + span - 1) / span;
            read_pos_ -= laps * span;
        } else if (read_pos_ < oldest) {
            const long long laps = (oldest - read_pos_ + span - 1) / span;
            read_pos_ += laps * span;
        }
 
        // Grid-snap, then a final hard clamp: the snap can push a position that
        // was exactly on an edge one cycle past it.
        read_pos_ = (read_pos_ / cycle_length_) * cycle_length_;
        read_pos_ = std::clamp(read_pos_, std::max<long long>(0, oldest),
                               std::max<long long>(0, newest));
    }

    void splice_weights(double theta, double& fade_in, double& fade_out) const noexcept {
        double linear_out = 0.0, linear_in = 0.0;
        double power_out = 0.0, power_in = 0.0;
        crossfade_gains(theta, CrossfadeGainLaw::EqualGain, linear_out, linear_in);
        crossfade_gains(theta, CrossfadeGainLaw::EqualPower, power_out, power_in);
        fade_in = linear_in + crossfade_shape_ * (power_in - linear_in);
        fade_out = linear_out + crossfade_shape_ * (power_out - linear_out);
    }

    void emit_grain() noexcept {
        // A grain may not start before the one BEFORE its predecessor has
        // finished. That is the two-grain overlap invariant stated directly,
        // and it is what every gain bound in this file rests on.
        //
        // At settled lengths it provably never binds: the grain two back ended
        // at `out_pos − 2S + N = out_pos − N + 2X`, and `2X ≤ N` is clamped, so
        // it ended at or before `out_pos` already. It binds only while the
        // length controls are being swept, which is exactly when a long grain
        // can still be sounding while short grains march past it at a hop that
        // has nothing to do with the length of the one in flight. Measured
        // without this guard, a full-range per-block sweep of every length
        // control reached 2.35× — well past the √2 the design proves — because
        // four and five grains were summing at once.
        const long long out_pos = std::max(grain_out_pos_, end_two_grains_ago_);
        grain_out_pos_ = out_pos;
        advance_read_cursor(out_pos);
        const long long in_pos = read_pos_;
 
        // This grain fades IN over exactly as many samples as its predecessor
        // faded OUT, not over its own `X`.
        //
        // At settled parameters those are the same number and this is a no-op.
        // They differ only while `crossfade_pct`, `cycle_hz` or `grain_periods`
        // are being automated, and that is precisely when it matters: the two
        // halves of a splice have to be COMPLEMENTARY for their weights to sum
        // to the bound the design proves. Give a splice a short fade-out and a
        // long fade-in and the outgoing grain is still at unity while the
        // incoming one is already most of the way up — the pair sums toward 2,
        // not √2, and the level jumps audibly on an automation move. Matching
        // the pair to the outgoing side keeps every splice complementary no
        // matter how violently the lengths are swept.
        // The first grain has no predecessor to match, so it uses the current
        // overlap. Deriving that from `grain_` rather than from a value seeded
        // at `reset()` is deliberate: a setter called AFTER `prepare()` — which
        // is the normal way to configure this thing — would otherwise leave the
        // seeded value describing the lengths that were in force during
        // `prepare`, and a later `reset()` would re-seed it from the CURRENT
        // lengths instead. Two renders that should be identical then differ in
        // their very first grain.
        const long long previous_fade =
            grain_ == 0 ? crossfade_length_ : previous_fade_out_;
        const long long fade_in = std::clamp<long long>(
            previous_fade, 1, std::max<long long>(1, grain_length_ - crossfade_length_));
        const long long fade_out = crossfade_length_;
        const long long interior_start = grain_length_ - fade_out;
 
        for (long long k = 0; k < grain_length_; ++k) {
            const auto source = static_cast<std::size_t>(in_pos + k) & ring_mask_;
            double weight = 1.0;
            if (k < fade_in) {
                double rising = 0.0, falling = 0.0;
                splice_weights(static_cast<double>(k) / static_cast<double>(fade_in),
                               rising, falling);
                weight = rising;
            } else if (k >= interior_start) {
                double rising = 0.0, falling = 0.0;
                splice_weights(static_cast<double>(k - interior_start) /
                                   static_cast<double>(fade_out),
                               rising, falling);
                weight = falling;
            }
            const auto slot = static_cast<std::size_t>(out_pos + k) & accum_mask_;
            if (accum_epoch_[slot] != accum_generation_) {
                accum_[slot] = 0.0;
                accum_epoch_[slot] = accum_generation_;
            }
            accum_[slot] += static_cast<double>(ring_[source]) * weight;
        }
 
        // Advance so the NEXT grain's fade-in lands exactly on this grain's
        // fade-out. The `fade_in` floor is the anti-pile-up guard: without it a
        // hop that has just been automated far below the previous grain's fade
        // length would let a third grain start before the first had finished,
        // and three overlapping grains break the two-term sum the bound rests
        // on. At settled parameters `fade_in ≤ X ≤ N/2 ≤ S`, so the floor never
        // binds and the hop is exactly `S`.
        grain_out_pos_ += std::max(hop_length_, fade_in);
        previous_fade_out_ = fade_out;
        end_two_grains_ago_ = end_previous_grain_;
        end_previous_grain_ = out_pos + grain_length_;
        ++grain_;
    }

    void update_lengths() noexcept {
        const double max_cycle_hz = sample_rate_ / static_cast<double>(kMinCycleSamples);
        const double hz = std::min(cycle_hz_, max_cycle_hz);
        cycle_length_ = std::max<long long>(kMinCycleSamples,
                                            std::llround(sample_rate_ / hz));
 
        grain_length_ = cycle_length_ * grain_periods_;
        if (max_grain_ > 0) grain_length_ = std::min(grain_length_, max_grain_);
 
        crossfade_length_ = std::llround(crossfade_pct_ * 0.01 *
                                         static_cast<double>(grain_length_));
        // The 50 % ceiling keeps the overlap-add to two grains, which is what
        // both the interior-unity assumption and the √2 bound rest on.
        crossfade_length_ = std::clamp<long long>(crossfade_length_, 1, grain_length_ / 2);
 
        hop_length_ = std::max<long long>(kMinHopSamples,
                                          grain_length_ - crossfade_length_);
 
        capture_window_ = std::max<long long>(
            grain_length_ + 1,
            static_cast<long long>(std::llround(
                units::ms_to_samples(capture_ms_, sample_rate_))));
        capture_window_ = std::min<long long>(
            capture_window_, static_cast<long long>(ring_mask_));
    }
 
    // Controls.
    double sample_rate_ = 48000.0;
    double cycle_hz_ = kCycleHzDefault;
    int grain_periods_ = kGrainPeriodsDefault;
    double crossfade_pct_ = kCrossfadePctDefault;
    double crossfade_shape_ = kCrossfadeShapeDefault;
    double stretch_ratio_ = kStretchRatioDefault;
    double capture_ms_ = kCaptureMsDefault;
    double mix_ = kMixDefault * 0.01;
    double output_gain_ = 1.0;
 
    // Resolved lengths.
    long long cycle_length_ = 240;
    long long grain_length_ = 240;
    long long crossfade_length_ = 60;
    long long hop_length_ = 180;
    long long capture_window_ = 96000;
    long long max_grain_ = 0;
 
    // Buffers and cursors. `total_written_` and `out_read_` are ABSOLUTE, which
    // is what makes the render independent of block size.
    std::vector<SampleType> ring_{};
    std::vector<double> accum_{};
    std::vector<std::uint64_t> accum_epoch_{};
    std::uint64_t accum_generation_ = 1u;
    std::size_t ring_mask_ = 0;
    std::size_t accum_mask_ = 0;
    std::size_t write_index_ = 0;
    long long total_written_ = 0;
    long long out_read_ = 0;
    long long grain_ = 0;

    long long grain_out_pos_ = 0;

    long long read_pos_ = 0;
    long long prev_out_pos_ = 0;

    long long previous_fade_out_ = 60;

    long long end_previous_grain_ = 0;
    long long end_two_grains_ago_ = 0;
};
 
using CyclicStretch = CyclicStretchT<float>;
using CyclicStretch64 = CyclicStretchT<double>;
 
}  // namespace pulp::signal