dco.hpp

Go to the documentation of this file.

#pragma once

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

enum class DcoDivider {
    integer_n,    
    fractional_n, 
};

struct DcoProfile {
    double master_clock_hz = 8'000'000.0;
    DcoDivider divider_scheme = DcoDivider::integer_n;
    int accumulator_bits = 24;
    int octave_shift = 0;
    double fine_tune_cents = 0.0;
};

class DcoOscillator {
public:
    void prepare(double sample_rate) noexcept {
        sample_rate_ = sample_rate > 0.0 ? sample_rate : sample_rate_;
        recompute();
    }
 
    void set_profile(const DcoProfile& profile) noexcept {
        profile_ = profile;
        recompute();
    }
    const DcoProfile& profile() const noexcept { return profile_; }

    void set_note_hz(double note_hz) noexcept {
        note_hz_ = note_hz;
        recompute();
    }
    double note_hz() const noexcept { return note_hz_; }
 
    // ── Shared shape stage (delegated; OSC-DCO shapes nothing) ───────────────
    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(); }
 
    // ── Derived timing facts (the exact equality-test targets) ───────────────
    long long divider_n() const noexcept {
        return profile_.divider_scheme == DcoDivider::integer_n ? divider_ : 0;
    }

    long long tuning_word() const noexcept {
        return profile_.divider_scheme == DcoDivider::fractional_n
                   ? static_cast<long long>(delta_)
                   : 0;
    }

    double realized_hz() const noexcept { return realized_hz_; }
    double effective_note_hz() const noexcept { return effective_note_hz_; }
    double detune_cents() const noexcept {
        return (realized_hz_ > 0.0 && effective_note_hz_ > 0.0)
                   ? 1200.0 * std::log2(realized_hz_ / effective_note_hz_)
                   : 0.0;
    }

    double quantization_bound_cents() const noexcept {
        if (!(profile_.master_clock_hz > 0.0) || !(effective_note_hz_ > 0.0)) return 0.0;
        const double ideal_divider =
            profile_.divider_scheme == DcoDivider::integer_n
                ? profile_.master_clock_hz / effective_note_hz_
                : static_cast<double>(two_b_) * effective_note_hz_ / profile_.master_clock_hz;
        // A note at/above 2·f_clk drives the ideal divider below ½ LSB, where the
        // divider clamp — not rounding — governs and the ½-LSB envelope no longer
        // applies. Report no bound there rather than a false one (unreachable for
        // any audio note below the master clock).
        if (!(ideal_divider > 0.5)) return 0.0;
        // Fractional-N also clamps the tuning word at the TOP (`Δ ≤ 2^B − 1`): as
        // the note approaches the clock the ideal `Δ` exceeds that and the clamp,
        // not ½-LSB rounding, governs the realized pitch — so the envelope does not
        // apply there either. Report no bound rather than a false tiny one.
        if (profile_.divider_scheme == DcoDivider::fractional_n &&
            !(ideal_divider < static_cast<double>(two_b_) - 0.5))
            return 0.0;
        return 1200.0 * std::log2(ideal_divider / (ideal_divider - 0.5));
    }

    int reset_intervals(std::span<long long> out) const noexcept {
        const int n = static_cast<int>(out.size());
        if (profile_.divider_scheme == DcoDivider::integer_n) {
            for (int i = 0; i < n; ++i) out[static_cast<std::size_t>(i)] = divider_;
            return n;
        }
        std::uint64_t residue = 0;
        for (int i = 0; i < n; ++i)
            out[static_cast<std::size_t>(i)] = advance_divider(residue);
        return n;
    }

    void reset(double phase = 0.0) noexcept {
        core_.reset(phase);
        residue_ = 0;
        if (profile_.divider_scheme == DcoDivider::fractional_n) {
            current_interval_ = advance_divider(residue_);
            clocks_to_reset_ = static_cast<double>(current_interval_);
        } else {
            current_interval_ = divider_;
            clocks_to_reset_ = 0.0;
        }
    }
 
    double phase() const noexcept { return core_.phase(); }

    double next() noexcept {
        if (profile_.divider_scheme == DcoDivider::integer_n) {
            // The natural wrap IS the reset; a constant increment keeps the period
            // exactly N clocks, so integer-N carries no jitter.
            return core_.next(increment_integer_);
        }
 
        // Fractional-N: the reset lands on the master-clock grid. The per-cycle
        // increment carries the phase to exactly 1.0 at the grid edge, so a single
        // clean step is emitted there; re-anchoring to 0 each cycle keeps the
        // ±1-clock jitter exact rather than accumulating slope error.
        const double increment = clocks_per_sample_ / static_cast<double>(current_interval_);
        double out;
        if (clocks_to_reset_ < clocks_per_sample_) {
            const double frac = clocks_to_reset_ > 0.0 ? clocks_to_reset_ / clocks_per_sample_ : 0.0;
            out = core_.next_synced(increment, frac, 0.0);
            const long long next_interval = advance_divider(residue_);
            clocks_to_reset_ += static_cast<double>(next_interval);
            current_interval_ = next_interval;
        } else {
            out = core_.next(increment);
        }
        clocks_to_reset_ -= clocks_per_sample_;
        return out;
    }
 
private:
    long long advance_divider(std::uint64_t& residue) const noexcept {
        // Smallest m with residue + m·Δ >= 2^B (ceil division), which alternates
        // between floor(2^B/Δ) and ceil(2^B/Δ) as the residue walks.
        const std::uint64_t m = (two_b_ - residue + delta_ - 1) / delta_;
        residue = residue + m * delta_ - two_b_;
        return static_cast<long long>(m);
    }
 
    void recompute() noexcept {
        clocks_per_sample_ = sample_rate_ > 0.0 ? profile_.master_clock_hz / sample_rate_ : 0.0;
 
        const int bits = profile_.accumulator_bits < 1    ? 1
                         : profile_.accumulator_bits > 32 ? 32
                                                          : profile_.accumulator_bits;
        two_b_ = std::uint64_t{1} << bits;
 
        effective_note_hz_ =
            note_hz_ * std::exp2(static_cast<double>(profile_.octave_shift) +
                                 profile_.fine_tune_cents / 1200.0);
 
        const double f_clk = profile_.master_clock_hz;
        if (profile_.divider_scheme == DcoDivider::integer_n) {
            const double ideal = effective_note_hz_ > 0.0 ? f_clk / effective_note_hz_ : 1.0;
            divider_ = static_cast<long long>(std::llround(ideal));
            if (divider_ < 1) divider_ = 1;
            realized_hz_ = divider_ > 0 ? f_clk / static_cast<double>(divider_) : 0.0;
            increment_integer_ = sample_rate_ > 0.0 ? realized_hz_ / sample_rate_ : 0.0;
            delta_ = 1;
        } else {
            const double ideal =
                f_clk > 0.0 ? effective_note_hz_ * static_cast<double>(two_b_) / f_clk : 1.0;
            long long d = static_cast<long long>(std::llround(ideal));
            if (d < 1) d = 1;
            if (static_cast<std::uint64_t>(d) >= two_b_) d = static_cast<long long>(two_b_ - 1);
            delta_ = static_cast<std::uint64_t>(d);
            realized_hz_ = f_clk * static_cast<double>(delta_) / static_cast<double>(two_b_);
            divider_ = 0;
            increment_integer_ = 0.0;
        }
 
        reset(phase());
    }
 
    VaOscillator core_;
 
    DcoProfile profile_{};
    double sample_rate_ = 48000.0;
    double note_hz_ = 261.625565; // C4.
 
    // Derived.
    double clocks_per_sample_ = 0.0;
    double effective_note_hz_ = 261.625565;
    double realized_hz_ = 0.0;
    long long divider_ = 1;         // integer-N.
    std::uint64_t delta_ = 1;       // fractional-N tuning word Δ.
    std::uint64_t two_b_ = 1u << 24;
    double increment_integer_ = 0.0;
 
    // Fractional-N running schedule.
    std::uint64_t residue_ = 0;
    long long current_interval_ = 1;
    double clocks_to_reset_ = 0.0;
};
 
} // namespace pulp::signal::osc