rng.hpp

Go to the documentation of this file.

#pragma once

 
#include <algorithm>
#include <cmath>
#include <cstdint>
 
namespace pulp::signal {

class Xorshift32 {
public:
    static constexpr std::uint32_t kDefaultSeed = 0x9E3779B9u;
 
    constexpr Xorshift32() = default;
    constexpr explicit Xorshift32(std::uint32_t s)
        : state_(s == 0u ? kDefaultSeed : s), seed_(state_) {}
 
    void seed(std::uint32_t s) {
        seed_ = (s == 0u ? kDefaultSeed : s);
        state_ = seed_;
        // The Box-Muller cache holds a deviate from the OLD stream. Leaving it
        // would make the first gaussian() after a re-seed come from the
        // sequence the caller just replaced.
        cached_ = 0.0f;
        has_cached_ = false;
    }

    void set_seed(std::uint32_t s) { seed(s); }

    void reset() {
        state_ = seed_ == 0u ? kDefaultSeed : seed_;
        cached_ = 0.0f;
        has_cached_ = false;
    }
 
    std::uint32_t next_u32() {
        std::uint32_t x = state_;
        if (x == 0u) x = kDefaultSeed;
        x ^= x << 13;
        x ^= x >> 17;
        x ^= x << 5;
        state_ = x;
        return x;
    }

    std::uint32_t next_uint() { return next_u32(); }

    float next_unipolar() {
        return static_cast<float>(next_u32() >> 8) * (1.0f / 16777216.0f);
    }

    float next_bipolar() { return next_unipolar() * 2.0f - 1.0f; }
 
    template <typename T>
    T next_bipolar() {
        return static_cast<T>(static_cast<double>(next_u32()) *
                                  (1.0 / 2147483648.0) -
                              1.0);
    }

    float next_unit() { return next_unipolar(); }
 
    template <typename T>
    T next_unit() {
        return static_cast<T>(static_cast<double>(next_u32()) *
                              (1.0 / 4294967296.0));
    }

    float gaussian() {
        if (has_cached_) {
            has_cached_ = false;
            return cached_;
        }
        float u1 = next_unipolar();
        if (u1 <= 0.0f) u1 = 1.0f / 16777216.0f;
        const float u2 = next_unipolar();
        const float radius = std::sqrt(-2.0f * std::log(u1));
        const float angle = 6.2831853071795864769f * u2;
        cached_ = radius * std::sin(angle);
        has_cached_ = true;
        return radius * std::cos(angle);
    }
 
    std::uint32_t state() const { return state_ == 0u ? kDefaultSeed : state_; }
 
    std::uint32_t seed() const { return seed_ == 0u ? kDefaultSeed : seed_; }
 
private:
    std::uint32_t state_ = kDefaultSeed;
    std::uint32_t seed_ = kDefaultSeed;
    float cached_ = 0.0f;
    bool has_cached_ = false;
};

constexpr std::uint64_t mix64(std::uint64_t z) {
    z ^= z >> 30;
    z *= 0xBF58476D1CE4E5B9ull;
    z ^= z >> 27;
    z *= 0x94D049BB133111EBull;
    z ^= z >> 31;
    return z;
}

constexpr std::uint64_t mix64(std::uint64_t seed,
                              std::uint64_t index,
                              std::uint64_t field = 0u) {
    return mix64(seed + index * 0x9E3779B97F4A7C15ull +
                 field * 0xBF58476D1CE4E5B9ull);
}

inline float unit_from(std::uint64_t key) {
    return static_cast<float>(mix64(key) >> 40) * (1.0f / 16777216.0f);
}
 
template <typename T>
inline T unit_from(std::uint64_t key) {
    return static_cast<T>(static_cast<double>(mix64(key) >> 11) *
                          (1.0 / 9007199254740992.0));
}
 
template <typename T = float>
inline T bipolar_from(std::uint64_t key) {
    return static_cast<T>(unit_from<double>(key) * 2.0 - 1.0);
}

inline constexpr std::uint64_t kGoldenGamma = 0x9E3779B97F4A7C15ull;

constexpr std::uint64_t rng_key(std::uint64_t purpose, std::uint64_t index) {
    return mix64(purpose * kGoldenGamma + index + kGoldenGamma);
}

template <typename SampleType = float>
class OuWalkT {
public:
    static constexpr SampleType kClamp = SampleType{1.5};
    static constexpr double kMinTheta = 1e-6;
    static constexpr double kMaxTheta = 1.0;
    static constexpr double kDefaultTheta = 0.05;

    void prepare(double update_rate_hz) {
        rate_ = update_rate_hz > 0.0 ? update_rate_hz : 1.0;
        dt_ = 1.0 / rate_;
        update_coefficients_();
    }
 
    void set_theta(double theta_per_second) {
        if (!std::isfinite(theta_per_second)) return;
        theta_ = std::clamp(theta_per_second, kMinTheta, kMaxTheta);
        update_coefficients_();
    }
    void set_sigma(double sigma) {
        if (!std::isfinite(sigma)) return;
        sigma_ = std::max(0.0, sigma);
        update_coefficients_();
    }
    void set_mu(double mu) {
        if (std::isfinite(mu)) mu_ = mu;
    }
    void seed(std::uint32_t s) { rng_.seed(s); }
    void set_seed(std::uint32_t s) { seed(s); }
 
    void set_correlation_time(double seconds, double sample_rate) {
        if (!(std::isfinite(seconds) && seconds > 0.0 &&
              std::isfinite(sample_rate) && sample_rate > 0.0))
            return;
        prepare(sample_rate);
        set_theta(1.0 / seconds);
    }
 
    void reset() {
        value_ = mu_;
        rng_.reset();
    }

    void reset_value() { value_ = mu_; }
 
    SampleType next() {
        // Exact transition: the decay is always in (0, 1), and the innovation
        // scale preserves the requested stationary deviation at every rate.
        const double diffusion = innovation_std_ * static_cast<double>(rng_.gaussian());
        value_ = std::clamp(mu_ + decay_ * (value_ - mu_) + diffusion,
                            -static_cast<double>(kClamp),
                            static_cast<double>(kClamp));
        return static_cast<SampleType>(value_);
    }
 
    SampleType current() const { return static_cast<SampleType>(value_); }
    SampleType value() const { return current(); }
    double theta() const { return theta_; }
    double sigma() const { return sigma_; }

    double stationary_std() const {
        return sigma_;
    }
 
private:
    void update_coefficients_() {
        decay_ = std::exp(-theta_ * dt_);
        // Exact discrete-time OU transition. `sigma_` is the requested
        // stationary standard deviation, so changing update rate or theta does
        // not silently change the modulation depth.
        innovation_std_ = sigma_ * std::sqrt(std::max(0.0, 1.0 - decay_ * decay_));
    }
 
    Xorshift32 rng_{};
    double value_ = 0.0;
    double mu_ = 0.0;
    double theta_ = kDefaultTheta;
    double sigma_ = 0.1;
    // An unprepared walk advances in unit time. Audio/control-rate users call
    // `prepare(rate)`, while the direct per-step Round 2 surface can state
    // theta without inheriting an arbitrary hidden 48 kHz assumption.
    double rate_ = 1.0;
    double dt_ = 1.0;
    double decay_ = std::exp(-kDefaultTheta);
    double innovation_std_ =
        0.1 * std::sqrt(1.0 - std::exp(-2.0 * kDefaultTheta));
};
 
using OuWalk = OuWalkT<float>;
using OuWalk64 = OuWalkT<double>;

template <typename SampleType = float>
class DriftT {
public:
    static constexpr int kDecimation = 32;
    static constexpr double kDefaultDepthPercent = 0.1;
    static constexpr double kDefaultCorrelationSeconds = 2.0;
 
    void prepare(double sample_rate) {
        const double sr = sample_rate > 0.0 ? sample_rate : 1.0;
        const double control_rate = sr / static_cast<double>(kDecimation);
        pitch_walk_.prepare(control_rate);
        fraction_walk_.prepare(control_rate);
        multiplier_walk_.prepare(control_rate);
        multiplier_walk_.set_theta(1.0 / correlation_s_);
        multiplier_walk_.set_sigma(depth_percent_ * 0.01);
        reset();
    }

    void set_cents(double cents) { cents_ = cents; }
 
    void set_theta(double theta_per_second) {
        pitch_walk_.set_theta(theta_per_second);
        fraction_walk_.set_theta(theta_per_second);
    }
 
    void set_sigma(double sigma) {
        pitch_walk_.set_sigma(sigma);
        fraction_walk_.set_sigma(sigma);
    }

    void seed(std::uint32_t s) {
        pitch_walk_.seed(static_cast<std::uint32_t>(mix64(s) & 0xFFFFFFFFu));
        fraction_walk_.seed(static_cast<std::uint32_t>(mix64(s ^ 0xA5A5A5A5u) >> 32));
        multiplier_walk_.seed(static_cast<std::uint32_t>(mix64(s ^ 0xC6BC2796u) & 0xFFFFFFFFu));
    }
 
    void set_seed(std::uint32_t s) { seed(s); }
 
    void set_depth_percent(double percent) {
        if (!std::isfinite(percent)) return;
        depth_percent_ = std::max(0.0, percent);
        multiplier_walk_.set_sigma(depth_percent_ * 0.01);
    }
 
    void set_correlation_time(double seconds) {
        if (!(std::isfinite(seconds) && seconds > 0.0)) return;
        correlation_s_ = seconds;
        multiplier_walk_.set_theta(1.0 / correlation_s_);
    }
 
    void reset() {
        pitch_walk_.reset();
        fraction_walk_.reset();
        multiplier_walk_.reset();
        counter_ = 0;
        pitch_prev_ = pitch_next_ = pitch_walk_.current();
        fraction_prev_ = fraction_next_ = fraction_walk_.current();
        multiplier_prev_ = multiplier_next_ = multiplier_walk_.current();
        pitch_now_ = pitch_prev_;
        fraction_now_ = fraction_prev_;
        multiplier_now_ = multiplier_prev_;
    }

    SampleType next() {
        if (counter_ == 0) {
            pitch_prev_ = pitch_next_;
            fraction_prev_ = fraction_next_;
            pitch_next_ = pitch_walk_.next();
            fraction_next_ = fraction_walk_.next();
            multiplier_prev_ = multiplier_next_;
            multiplier_next_ = multiplier_walk_.next();
        }
        const SampleType t = static_cast<SampleType>(counter_)
                             / static_cast<SampleType>(kDecimation);
        pitch_now_ = pitch_prev_ + (pitch_next_ - pitch_prev_) * t;
        fraction_now_ = fraction_prev_ + (fraction_next_ - fraction_prev_) * t;
        multiplier_now_ = multiplier_prev_ + (multiplier_next_ - multiplier_prev_) * t;
        if (++counter_ >= kDecimation) counter_ = 0;
        return static_cast<SampleType>(
            std::max(kMinMultiplier, 1.0 + static_cast<double>(multiplier_now_)));
    }

    SampleType pitch_factor() const {
        return static_cast<SampleType>(
            std::exp2(cents_ * static_cast<double>(pitch_now_) / 1200.0));
    }

    SampleType fraction() const {
        return std::clamp(fraction_now_, SampleType{-1}, SampleType{1});
    }
 
private:
    static constexpr double kMinMultiplier = 1e-3;
    OuWalkT<SampleType> pitch_walk_{};
    OuWalkT<SampleType> fraction_walk_{};
    OuWalkT<SampleType> multiplier_walk_{};
    double cents_ = 5.0;
    double depth_percent_ = kDefaultDepthPercent;
    double correlation_s_ = kDefaultCorrelationSeconds;
    int counter_ = 0;
    SampleType pitch_prev_{};
    SampleType pitch_next_{};
    SampleType fraction_prev_{};
    SampleType fraction_next_{};
    SampleType pitch_now_{};
    SampleType fraction_now_{};
    SampleType multiplier_prev_{};
    SampleType multiplier_next_{};
    SampleType multiplier_now_{};
};
 
using Drift = DriftT<float>;
using Drift64 = DriftT<double>;
 
} // namespace pulp::signal