va.hpp

Go to the documentation of this file.

#pragma once

 
#include "blep.hpp"
#include "phase.hpp"
 
#include <cmath>
#include <cstddef>
#include <numbers>
#include <span>
 
namespace pulp::signal::osc {

enum class VaShape {
    sine,     
    saw,      
    square,   
    triangle, 
};

inline int threshold_crossings(double phase_before_advance, double increment,
                               double threshold, std::span<double> out) noexcept {
    if (!std::isfinite(increment) || increment == 0.0) return 0;
 
    // Shift so the threshold sits on an integer, then the crossings are the
    // integers strictly between the endpoints.
    const double q0 = phase_before_advance - threshold;
    const double first = std::floor(q0);
    const double last = std::floor(q0 + increment);
    if (!std::isfinite(first) || !std::isfinite(last)) return 0;
 
    const bool forward = increment > 0.0;
    const double magnitude = forward ? (last - first) : (first - last);
    if (!(magnitude >= 1.0)) return 0;
 
    const int capacity = static_cast<int>(out.size());
    const int count = magnitude > static_cast<double>(capacity)
                          ? capacity
                          : static_cast<int>(magnitude);
 
    for (int k = 0; k < count; ++k) {
        const double level =
            forward ? first + 1.0 + static_cast<double>(k) : first - static_cast<double>(k);
        out[static_cast<std::size_t>(k)] = (level - q0) / increment;
    }
    return count;
}

class VaOscillator {
public:
    void set_shape(VaShape shape) noexcept { shape_ = shape; }
    VaShape shape() const noexcept { return shape_; }

    void set_pulse_width(double width) noexcept { pulse_width_ = clamp_unit(width); }
    double pulse_width() const noexcept { return pulse_width_; }

    void reset(double phase = 0.0) noexcept {
        phase_.reset(phase);
        carry_ = 0.0;
    }

    double phase() const noexcept { return phase_.phase(); }

    double next(double increment) noexcept {
        const double entry_phase = phase_.phase();
        double out = value(entry_phase) + carry_;
        carry_ = 0.0;
 
        phase_.advance(increment);
        correct_events(increment, out);
        // No sync, so the phase ran the whole sample at one rate from one place:
        // the sample is a single segment.
        scan_threshold(entry_phase, increment, increment, 0.0, 1.0, out);
        return out;
    }

    double next_synced(double increment, double sync_frac, double sync_phase = 0.0) noexcept {
        const double entry_phase = phase_.phase();
        double out = value(entry_phase) + carry_;
        carry_ = 0.0;
 
        phase_.advance_synced(increment, sync_frac, sync_phase);
        correct_events(increment, out);
 
        // A sync SPLITS the sample in two. The phase runs from `entry_phase` over
        // sub-sample span [0, f], is then replaced outright, and runs from the
        // target over [f, 1]. An internal threshold must therefore be scanned
        // once per segment, from that segment's own starting phase: a single scan
        // of (entry_phase, increment) describes a trajectory the phase never
        // took, and would place the apex where it never was while missing the one
        // it actually crossed after the reset.
        //
        // This mirrors `advance_synced`'s own segmentation — the same split, the
        // same sub-sample spans — because that split IS the semantics, not an
        // implementation detail to approximate.
        const double f = clamp_unit(sync_frac);
        scan_threshold(entry_phase, increment * f, increment, 0.0, f, out);
        scan_threshold(wrap_unit(sync_phase), increment * (1.0 - f), increment, f, 1.0 - f, out);
        return out;
    }
 
private:
    void correct_events(double increment, double& out) noexcept {
        for (const PhaseEvent& e : phase_.events()) {
            const double value_step = value(e.phase_after) - value(e.phase_before);
            const double slope_step =
                (slope_per_cycle(e.phase_after) - slope_per_cycle(e.phase_before)) * increment;
            add_break(e.frac, value_step, slope_step, out);
        }
    }

    void scan_threshold(double start_phase, double delta, double rate, double t0, double span,
                        double& out) noexcept {
        Threshold t;
        if (!internal_threshold(t)) return;
 
        double fracs[PhaseAccumulator::max_events_per_sample];
        const int n = threshold_crossings(start_phase, delta, t.phase, std::span<double>(fracs));
        if (n == 0) return;
 
        // Time order decides which side is "before": going forward the phase
        // arrives from below the threshold, going backward from above.
        const bool forward = delta > 0.0;
        const double value_step =
            forward ? (t.value_above - t.value_below) : (t.value_below - t.value_above);
        const double slope_step =
            (forward ? (t.slope_above - t.slope_below) : (t.slope_below - t.slope_above)) * rate;
 
        for (int k = 0; k < n; ++k)
            add_break(t0 + fracs[k] * span, value_step, slope_step, out);
    }

    struct Threshold {
        double phase = 0.0;
        double value_below = 0.0;
        double value_above = 0.0;
        double slope_below = 0.0; 
        double slope_above = 0.0;
    };
 
    bool internal_threshold(Threshold& t) const noexcept {
        switch (shape_) {
            case VaShape::square:
                // No edge at all when the pulse fills the period or none of it.
                if (!(pulse_width_ > 0.0 && pulse_width_ < 1.0)) return false;
                t = {pulse_width_, 1.0, -1.0, 0.0, 0.0};
                return true;
            case VaShape::triangle:
                // The apex: value is continuous, the slope reverses.
                t = {0.5, 1.0, 1.0, 4.0, -4.0};
                return true;
            case VaShape::sine:
            case VaShape::saw:
                return false;
        }
        return false;
    }

    double value(double p) const noexcept {
        switch (shape_) {
            case VaShape::sine:
                // sin is periodic, so its limit at the top of the period is
                // sin(0) = 0 EXACTLY. `std::sin(2*pi)` returns -2.4e-16 instead,
                // because 2*pi is not exactly representable — and handing that to
                // the corrector applies a BLEP at every wrap of a shape that does
                // not step there at all. Mapping the sentinel to its exact limit
                // keeps a wrap's step exactly zero, which the corrector then
                // skips outright, leaving a plain sine bit-identical to
                // `std::sin`. That bit-exactness is the null proving the
                // corrector is inert on a continuous shape, and it is worth more
                // than the 2.4e-16 it costs to protect.
                //
                // A SYNC on a sine is a different matter and IS corrected: the
                // reset jumps to an unrelated phase, so the value genuinely
                // steps. It is the wrap that is continuous here, not the shape.
                return p >= 1.0 ? 0.0 : std::sin(2.0 * std::numbers::pi * p);
            case VaShape::saw:
                // Limit at 1 is +1, which is what `2p - 1` gives at p = 1.
                return 2.0 * p - 1.0;
            case VaShape::square:
                // A width of exactly 1 fills the period: the shape is constantly
                // high, so its limit at the top is high too. Testing `p < width`
                // alone would report LOW at the p = 1 sentinel and manufacture a
                // +2 step at a wrap that has no edge on either side of it.
                return (pulse_width_ >= 1.0 || p < pulse_width_) ? 1.0 : -1.0;
            case VaShape::triangle:
                // Limit at 1 is -1, which is what the descending leg gives at 1.
                return p < 0.5 ? (4.0 * p - 1.0) : (3.0 - 4.0 * p);
        }
        return 0.0;
    }

    double slope_per_cycle(double p) const noexcept {
        switch (shape_) {
            case VaShape::sine:
                return 2.0 * std::numbers::pi * std::cos(2.0 * std::numbers::pi * p);
            case VaShape::saw:
                return 2.0;
            case VaShape::triangle:
                return p < 0.5 ? 4.0 : -4.0;
            case VaShape::square:
                return 0.0;
        }
        return 0.0;
    }

    void add_break(double frac, double value_step, double slope_step, double& out) noexcept {
        if (value_step != 0.0) {
            const Correction c = poly_blep(frac, value_step);
            out += c.before;
            carry_ += c.after;
        }
        if (slope_step != 0.0) {
            const Correction c = poly_blamp(frac, slope_step);
            out += c.before;
            carry_ += c.after;
        }
    }
 
    static double clamp_unit(double x) noexcept {
        if (!(x > 0.0)) return 0.0;
        if (x > 1.0) return 1.0;
        return x;
    }

    static double wrap_unit(double x) noexcept {
        const double w = x - std::floor(x);
        return (w >= 0.0 && w < 1.0) ? w : 0.0;
    }
 
    PhaseAccumulator phase_;
    double carry_ = 0.0;
    double pulse_width_ = 0.5;
    VaShape shape_ = VaShape::saw;
};
 
} // namespace pulp::signal::osc