vco.hpp

Go to the documentation of this file.

#pragma once

 
#include "blep.hpp"
#include "phase.hpp"
#include "va.hpp"
 
#include <pulp/signal/dc_blocker.hpp>
#include <algorithm>
 
#include <array>
#include <cmath>
#include <cstddef>
#include <cstdint>
#include <limits>
#include <numbers>
 
namespace pulp::signal::osc {

struct VcoTuning {
    double reference_hz = 261.625565; // C4.
    double volts_per_octave = 1.0;
    double tune_offset_cents = 0.0;
    double scale_error = 1.0;
    double hf_compression = 0.0;
    double hf_knee_octaves = 3.0;

    double frequency_hz(double control_volts) const noexcept {
        double octaves = (control_volts / volts_per_octave) * scale_error
                       + tune_offset_cents / 1200.0;
        if (hf_compression > 0.0 && octaves > hf_knee_octaves)
            octaves -= hf_compression * (octaves - hf_knee_octaves);
        return reference_hz * std::exp2(octaves);
    }

    double phase_increment(double control_volts, double sample_rate) const noexcept {
        return frequency_hz(control_volts) / sample_rate;
    }
};

struct WaveshaperParams {
    double amount = 0.0;
    double drive = 1.0;
    double asymmetry = 0.0;
};

class VcoOscillator {
public:
    VcoOscillator() noexcept { update_drift_coeffs(); }
 
    void prepare(double sample_rate) noexcept {
        sample_rate_ = sample_rate > 0.0 ? sample_rate : sample_rate_;
        update_ac_pole();
        update_drift_coeffs();
    }
 
    void set_shape(VaShape shape) noexcept { core_.set_shape(shape); }
    VaShape shape() const noexcept { return core_.shape(); }
 
    void set_pulse_width(double width) noexcept { core_.set_pulse_width(width); }
    double pulse_width() const noexcept { return core_.pulse_width(); }
 
    // ── Integrator-leak bow ────────────────────────────────────────────────
    void set_bow(double curvature) noexcept { bow_ = curvature > 0.0 ? curvature : 0.0; }
    double bow() const noexcept { return bow_; }
 
    // ── Per-shape waveshaper ───────────────────────────────────────────────
    void set_waveshaper(VaShape shape, const WaveshaperParams& params) noexcept {
        // `amount` is a blend, so it is clamped to [0, 1] like every other
        // unit-range control here: outside that range it extrapolates the blend
        // (amount 2 reaches 3.0 on a full-scale shape) and breaks the [-1, 1]
        // output contract the boundedness gate relies on.
        WaveshaperParams p = params;
        p.amount = std::clamp(p.amount, 0.0, 1.0);
        shapers_[static_cast<std::size_t>(shape)] = p;
    }
    WaveshaperParams waveshaper(VaShape shape) const noexcept {
        return shapers_[static_cast<std::size_t>(shape)];
    }
 
    // ── Level vs pitch ─────────────────────────────────────────────────────
    void set_level_tilt(double db_per_octave) noexcept { level_tilt_db_per_octave_ = db_per_octave; }
    double level_tilt() const noexcept { return level_tilt_db_per_octave_; }

    void set_core_reset_time(double seconds) noexcept {
        core_reset_seconds_ = seconds > 0.0 ? seconds : 0.0;
    }
    double core_reset_time() const noexcept { return core_reset_seconds_; }
 
    // ── AC coupling (output DC blocker) ────────────────────────────────────
    void set_ac_coupling(double corner_hz) noexcept {
        ac_corner_hz_ = corner_hz > 0.0 ? corner_hz : 0.0;
        update_ac_pole();
    }
    double ac_coupling() const noexcept { return ac_corner_hz_; }
    double ac_pole() const noexcept { return ac_enabled_ ? ac_pole_ : 1.0; }
 
    // ── Drift / jitter (seeded pitch noise) ─────────────────────────────────
    void set_drift_depth(double depth_cents) noexcept { drift_depth_ = depth_cents; }
    double drift_depth() const noexcept { return drift_depth_; }

    void set_drift_rate_hz(double corner_hz) noexcept {
        drift_rate_hz_ = corner_hz > 0.0 ? corner_hz : 0.0;
        update_drift_coeffs();
    }
    double drift_rate_hz() const noexcept { return drift_rate_hz_; }

    void set_jitter_depth(double depth_cents) noexcept { jitter_depth_ = depth_cents; }
    double jitter_depth() const noexcept { return jitter_depth_; }

    void set_seed(std::uint64_t seed) noexcept {
        seed_ = seed;
        seed_streams();
    }
    std::uint64_t seed() const noexcept { return seed_; }

    static constexpr double kDefaultDriftRateHz = 0.4;

    void reset(double phase = 0.0) noexcept {
        core_.reset(phase);
        dc_blocker_.reset();
        drift_state_ = 0.0;
        seed_streams();
    }
 
    double phase() const noexcept { return core_.phase(); }

    double next(double increment) noexcept {
        const double modulated = increment * pitch_noise_factor();
        double v = core_.next(modulated);
        v = apply_bow(v);
        v = apply_waveshaper(v);
        v *= level_for(modulated);
        return ac_enabled_ ? dc_blocker_.process(v) : v;
    }

    static constexpr double kLevelReferenceHz = 1000.0;
 
private:
    double apply_bow(double v) const noexcept {
        if (bow_ == 0.0 || core_.shape() != VaShape::saw) return v;
        // Saw value in [-1, 1] ↔ phase in [0, 1); clamp the tiny polyBLEP
        // overshoot so the charge curve stays monotone and bounded.
        double p = 0.5 * (v + 1.0);
        p = p < 0.0 ? 0.0 : (p > 1.0 ? 1.0 : p);
        // Normalized capacitor charge (1 - e^{-k p}) / (1 - e^{-k}); k → 0 is the
        // linear ramp, guarded above by the exact bypass.
        const double charged = (1.0 - std::exp(-bow_ * p)) / (1.0 - std::exp(-bow_));
        return 2.0 * charged - 1.0;
    }

    double apply_waveshaper(double v) const noexcept {
        const WaveshaperParams& w = shapers_[static_cast<std::size_t>(core_.shape())];
        if (w.amount == 0.0) return v;
        // A tanh saturator biased by `asymmetry`, mapped so the full-scale ends
        // ±1 land exactly on ±1 and the curve is DC-free (its midpoint removed).
        // Anchoring on the two ends rather than a single-point norm is what keeps
        // the shaped signal inside [-1, 1] for any bias — a monotone curve pinned
        // at both endpoints cannot overshoot between them.
        const double lo = std::tanh(-w.drive + w.asymmetry);
        const double hi = std::tanh(w.drive + w.asymmetry);
        const double mid = 0.5 * (hi + lo);   // the DC the bias introduces.
        const double half = 0.5 * (hi - lo);  // > 0 for drive > 0.
        const double shaped =
            half != 0.0 ? (std::tanh(w.drive * v + w.asymmetry) - mid) / half : v;
        return (1.0 - w.amount) * v + w.amount * shaped;
    }

    double level_for(double increment) const noexcept {
        if (level_tilt_db_per_octave_ == 0.0 && core_reset_seconds_ == 0.0) return 1.0;
        const double frequency = std::fabs(increment) * sample_rate_;
        if (!(frequency > 0.0)) return 1.0;
 
        double gain = 1.0;
        if (level_tilt_db_per_octave_ != 0.0) {
            const double octaves = std::log2(frequency / kLevelReferenceHz);
            gain *= std::pow(10.0, -level_tilt_db_per_octave_ * octaves / 20.0);
        }
        if (core_reset_seconds_ != 0.0) {
            const double eaten = core_reset_seconds_ * frequency; // fraction of the period.
            gain *= eaten < 1.0 ? (1.0 - eaten) : 0.0;
        }
        return gain;
    }

    double pitch_noise_factor() noexcept {
        if (drift_depth_ == 0.0 && jitter_depth_ == 0.0) return 1.0;
 
        double cents = 0.0;
        if (drift_depth_ != 0.0) {
            drift_state_ = drift_pole_ * drift_state_ + drift_norm_ * drift_rng_.next_gaussian();
            cents += drift_depth_ * drift_state_;
        }
        if (jitter_depth_ != 0.0)
            cents += jitter_depth_ * jitter_rng_.next_gaussian();
 
        return std::exp2(cents / 1200.0);
    }

    struct NoiseSource {
        std::uint64_t state = 0;
 
        void seed(std::uint64_t s) noexcept { state = s; }
 
        std::uint64_t next_u64() noexcept {
            std::uint64_t z = (state += 0x9E3779B97F4A7C15ull);
            z = (z ^ (z >> 30)) * 0xBF58476D1CE4E5B9ull;
            z = (z ^ (z >> 27)) * 0x94D049BB133111EBull;
            return z ^ (z >> 31);
        }

        double next_uniform() noexcept {
            const double u =
                static_cast<double>(next_u64() >> 11) * (1.0 / 9007199254740992.0);
            return u > 0.0 ? u : std::numeric_limits<double>::min();
        }
 
        double next_gaussian() noexcept {
            const double u1 = next_uniform();
            const double u2 = next_uniform();
            return std::sqrt(-2.0 * std::log(u1)) *
                   std::cos(2.0 * std::numbers::pi * u2);
        }
    };

    void seed_streams() noexcept {
        drift_rng_.seed(seed_);
        jitter_rng_.seed(seed_ ^ 0xD1B54A32D192ED03ull);
    }

    void update_drift_coeffs() noexcept {
        // Lower rate is a slower wander, so the pole rises toward 1 as the rate
        // falls; the limit as rate -> 0 is a frozen walk (pole 1), which the
        // unit-variance gain below turns into no wander at all. A non-positive
        // rate takes that limit rather than the opposite extreme: mapping it to
        // pole 0 would make drift full-depth per-sample white noise -- the
        // FASTEST wander -- so a caller passing 0 to mean "off" or "very slow"
        // would get audible hiss-like FM instead.
        drift_pole_ = drift_rate_hz_ > 0.0
                          ? std::exp(-2.0 * std::numbers::pi * drift_rate_hz_ / sample_rate_)
                          : 1.0;
        drift_norm_ = std::sqrt(1.0 - drift_pole_ * drift_pole_);
    }
 
    void update_ac_pole() noexcept {
        ac_enabled_ = ac_corner_hz_ > 0.0;
        if (ac_enabled_) {
            ac_pole_ = std::exp(-2.0 * std::numbers::pi * ac_corner_hz_ / sample_rate_);
            dc_blocker_.set_pole(ac_pole_);
        } else {
            ac_pole_ = 1.0;
        }
    }
 
    VaOscillator core_;
    double sample_rate_ = 48000.0;
 
    double bow_ = 0.0;
    std::array<WaveshaperParams, 4> shapers_{};
 
    double level_tilt_db_per_octave_ = 0.0;
    double core_reset_seconds_ = 0.0;
 
    double ac_corner_hz_ = 0.0;
    double ac_pole_ = 1.0;
    bool ac_enabled_ = false;
    pulp::signal::DcBlocker<double> dc_blocker_;
 
    double drift_depth_ = 0.0;
    double jitter_depth_ = 0.0;
    double drift_rate_hz_ = kDefaultDriftRateHz;
    double drift_pole_ = 0.0;
    double drift_norm_ = 1.0;
    double drift_state_ = 0.0;
    std::uint64_t seed_ = 0x0123456789ABCDEFull;
    NoiseSource drift_rng_{seed_};
    NoiseSource jitter_rng_{seed_ ^ 0xD1B54A32D192ED03ull};
};
 
} // namespace pulp::signal::osc