harmony_engine.hpp

Go to the documentation of this file.

#pragma once
 
#include <pulp/music/pitch.hpp>

 
#include <pulp/signal/dc_blocker.hpp>
#include <pulp/signal/delay_line.hpp>
#include <pulp/signal/denormal.hpp>
#include <pulp/signal/pitch_shifter.hpp>
#include <pulp/signal/rng.hpp>
#include <pulp/signal/slew_limiter.hpp>
#include <pulp/signal/units.hpp>
#include <pulp/signal/yin_tracker.hpp>
 
#include <algorithm>
#include <cmath>
#include <cstddef>
#include <cstdint>
#include <vector>
 
namespace pulp::signal {

enum class ScaleType : std::uint8_t {
    major = 0,         
    natural_minor = 1, 
    dorian = 2,        
    phrygian = 3,      
    lydian = 4,        
    mixolydian = 5,    
    harmonic_minor = 6, 
    melodic_minor = 7,  
    major_pentatonic = 8, 
    minor_pentatonic = 9, 
};

inline constexpr int kScaleCount = 10;

inline constexpr std::uint16_t kScaleTable[kScaleCount] = {
    music::scale_intervals(music::kPulpSignalScales[0].scale)->mask(),
    music::scale_intervals(music::kPulpSignalScales[1].scale)->mask(),
    music::scale_intervals(music::kPulpSignalScales[2].scale)->mask(),
    music::scale_intervals(music::kPulpSignalScales[3].scale)->mask(),
    music::scale_intervals(music::kPulpSignalScales[4].scale)->mask(),
    music::scale_intervals(music::kPulpSignalScales[5].scale)->mask(),
    music::scale_intervals(music::kPulpSignalScales[6].scale)->mask(),
    music::scale_intervals(music::kPulpSignalScales[7].scale)->mask(),
    music::scale_intervals(music::kPulpSignalScales[8].scale)->mask(),
    music::scale_intervals(music::kPulpSignalScales[9].scale)->mask(),
};

enum class OffScalePolicy : std::uint8_t {
    nearest_lower,  
    nearest,        
    mute_wet,       
};

struct DiatonicMapping {
    int shift_semitones = 0;  
    int input_midi = 0;       
    int snapped_midi = 0;     
    int target_midi = 0;      
    int degree = 0;           
    bool chromatic = false;   
    bool clamped = false;     
};
 
// ── Stage 2 — diatonic interval mapping ───────────────────────────────────

template <typename SampleType = float>
class DiatonicMapT {
public:
    static constexpr int kMaxDegrees = 12;

    static constexpr int kIntervalStepsMax = 14;

    static constexpr int kShiftSemitonesMax = 12;
 
    DiatonicMapT() { rebuild(); }

    void set_key(int root_pitch_class) {
        root_ = ((root_pitch_class % 12) + 12) % 12;
    }
    int key() const { return root_; }
 
    void set_scale(ScaleType scale) {
        scale_ = scale;
        rebuild();
    }
    ScaleType scale() const { return scale_; }
 
    void set_off_scale_policy(OffScalePolicy policy) { policy_ = policy; }
    OffScalePolicy off_scale_policy() const { return policy_; }

    int degree_count() const { return degree_count_; }

    int degree_semitone(int index) const {
        if (index < 0 || index >= degree_count_) return 0;
        return degrees_[static_cast<std::size_t>(index)];
    }

    DiatonicMapping map_midi(int midi_note, int interval_steps) const {
        DiatonicMapping out;
        out.input_midi = midi_note;
        if (degree_count_ <= 0) return out;
 
        const int steps = std::clamp(interval_steps, -kIntervalStepsMax, kIntervalStepsMax);
 
        // Pitch class within the key, and the octave block it sits in.
        const int rel = midi_note - root_;
        const int pc = ((rel % 12) + 12) % 12;
        const int octave = (rel - pc) / 12;
 
        // Snap to a degree.
        int degree = 0;
        bool chromatic = true;
        for (int i = 0; i < degree_count_; ++i) {
            if (degrees_[static_cast<std::size_t>(i)] == pc) {
                degree = i;
                chromatic = false;
                break;
            }
        }
        if (chromatic) degree = snap_degree(pc);
 
        out.chromatic = chromatic;
        out.degree = degree;
        out.snapped_midi =
            root_ + 12 * octave + degrees_[static_cast<std::size_t>(degree)];
 
        // Walk `steps` entries up (or down) the scale, wrapping the octave on
        // the degree COUNT — which is what makes a "third" a scale relation
        // rather than a fixed transposition.
        const int target_degree = degree + steps;
        const int wrapped = floor_div(target_degree, degree_count_);
        const int index = target_degree - wrapped * degree_count_;
        out.target_midi = root_ + 12 * (octave + wrapped) +
                          degrees_[static_cast<std::size_t>(index)];
 
        int shift = out.target_midi - out.snapped_midi;
        const int limited = std::clamp(shift, -kShiftSemitonesMax, kShiftSemitonesMax);
        out.clamped = limited != shift;
        out.shift_semitones = limited;
        return out;
    }

    DiatonicMapping map_hz(double f0_hz, int interval_steps) const {
        if (!(f0_hz > 0.0)) return DiatonicMapping{};
        const int midi = static_cast<int>(
            std::lround(units::hz_to_midi(f0_hz)));
        return map_midi(midi, interval_steps);
    }
 
private:
    static int floor_div(int a, int b) {
        const int q = a / b;
        return (a % b != 0 && ((a < 0) != (b < 0))) ? q - 1 : q;
    }

    int snap_degree(int pc) const {
        int lower = degree_count_ - 1;  // wraps to the degree below the octave
        for (int i = degree_count_ - 1; i >= 0; --i) {
            if (degrees_[static_cast<std::size_t>(i)] <= pc) {
                lower = i;
                break;
            }
        }
        if (policy_ != OffScalePolicy::nearest) return lower;
 
        const int upper = (lower + 1) % degree_count_;
        const int upper_pc = degrees_[static_cast<std::size_t>(upper)] +
                             (upper == 0 ? 12 : 0);
        const int down = pc - degrees_[static_cast<std::size_t>(lower)];
        const int up = upper_pc - pc;
        return up < down ? upper : lower;  // ties → lower
    }
 
    void rebuild() {
        const auto index = static_cast<std::size_t>(scale_);
        const std::uint16_t mask =
            index < static_cast<std::size_t>(kScaleCount) ? kScaleTable[index]
                                                          : kScaleTable[0];
        degree_count_ = 0;
        for (int k = 0; k < 12; ++k) {
            if ((mask >> k) & 1u) {
                degrees_[static_cast<std::size_t>(degree_count_)] = k;
                ++degree_count_;
            }
        }
        if (degree_count_ == 0) {  // a mask with no root is not a scale
            degrees_[0] = 0;
            degree_count_ = 1;
        }
    }
 
    int root_ = 0;
    ScaleType scale_ = ScaleType::major;
    OffScalePolicy policy_ = OffScalePolicy::nearest_lower;
    int degrees_[kMaxDegrees] = {};
    int degree_count_ = 0;
};
 
using DiatonicMap = DiatonicMapT<float>;
using DiatonicMap64 = DiatonicMapT<double>;
 
// ── Stages 3 & 4 — voices, detune, glide ──────────────────────────────────

template <typename SampleType = float>
class HarmonyEngineT {
public:
    static constexpr int kMaxVoices = 2;

    static constexpr double kCrossfadeMsDefault = 20.0;
    static constexpr double kCrossfadeMsMin = 10.0;
    static constexpr double kCrossfadeMsMax = 50.0;

    static constexpr double kGlideMsDefault = 60.0;
    static constexpr double kGlideMsMax = 500.0;

    static constexpr double kVoiceMuteMs = 10.0;

    static constexpr double kDetuneMaxCents = 50.0;

    static constexpr double kHumanizeDefault = 0.0;
    static constexpr double kHumanizeMaxCents = 15.0;

    static constexpr int kV1IntervalDefault = 2;
    static constexpr int kV2IntervalDefault = 4;

    static constexpr double kLevelMinDb = -60.0;
    static constexpr double kLevelMaxDb = 6.0;
    static constexpr double kLevelMaxLinear = 1.9952623149688795;

    static constexpr double kWorstCaseGain =
        kLevelMaxLinear *
        (1.0 + kMaxVoices * PitchShifterT<double>::kDcBlockerPeakGain);

    static constexpr double kMuteFloor = 1e-5;
 
    HarmonyEngineT() {
        for (int v = 0; v < kMaxVoices; ++v) {
            glide_[v].set_mode(SlewMode::linear);
            glide_[v].set_time_ms(kGlideMsDefault);
            interval_[v] = v == 0 ? kV1IntervalDefault : kV2IntervalDefault;
            // Independent drift streams: one seed for both voices would make
            // two "humanised" voices wander in lockstep, which is exactly the
            // decorrelation the feature exists to provide.
            drift_[v].set_seed(kDriftSeedBase ^
                               (0x9E3779B9u * static_cast<std::uint32_t>(v + 1)));
            drift_[v].set_depth_percent(0.0);
        }
        mute_.set_mode(SlewMode::linear);
        mute_.set_time_ms(kVoiceMuteMs);
        enabled_[0] = true;
        enabled_[1] = false;
    }

    void prepare(double sample_rate) {
        sample_rate_ = (std::isfinite(sample_rate) && sample_rate > 0.0)
                           ? sample_rate
                           : 48000.0;
        tracker_.prepare(sample_rate_);
 
        for (int v = 0; v < kMaxVoices; ++v) {
            auto& shifter = shifter_[static_cast<std::size_t>(v)];
            shifter.prepare(sample_rate_);
            shifter.set_shift_source(ShiftSource::direct);
            shifter.set_window_ms(crossfade_ms_);
            // The engine owns the glide, in cents, so the shifter's own
            // semitone-domain slew must be out of the way — two serial glides
            // would neither reach the target in `glide_ms` nor be linear in
            // cents.
            shifter.set_glide_ms(0.0, 0.0);
            shifter.set_mix(1.0);
            shifter.set_interp(interp_);
            glide_[v].prepare(sample_rate_);
            drift_[v].prepare(sample_rate_);
        }
        mute_.prepare(sample_rate_);
 
        update_alignment();
        // Sized for the WORST case, not the current one: `set_crossfade_ms` can
        // grow the shifter's latency after `prepare`, and with a high `kF0Min`
        // the shifter's window can exceed the tracker's. Sizing to the current
        // latency would leave a later crossfade change reading past the end of
        // the ring — which wraps rather than faults, so it would surface as a
        // mysterious wrong-sounding dry path rather than as a crash.
        const int worst_shifter_latency = static_cast<int>(
            std::lround(kCrossfadeMsMax * sample_rate_ / 2000.0));
        align_.prepare(std::max(tracker_.latency_samples(), worst_shifter_latency) + 1);
        reset();
    }
 
    void reset() {
        tracker_.reset();
        align_.reset();
        dc_.reset();
        mute_.set_immediate(0.0);
        for (int v = 0; v < kMaxVoices; ++v) {
            shifter_[static_cast<std::size_t>(v)].reset();
            drift_[v].reset();
            shift_semitones_[v] = 0;
            glide_[v].set_immediate(target_cents(v));
            apply_ratio(v, glide_[v].value());
        }
    }

    void discard_history() noexcept {
        tracker_.discard_history();
        align_.discard_history();
        dc_.reset();
        mute_.set_immediate(0.0);
        for (int v = 0; v < kMaxVoices; ++v) {
            shifter_[static_cast<std::size_t>(v)].discard_history();
            drift_[v].reset();
            shift_semitones_[v] = 0;
            glide_[v].set_immediate(target_cents(v));
            apply_ratio(v, glide_[v].value());
        }
    }
 
    // ── Control surface ───────────────────────────────────────────────────
 
    void set_key(int root_pitch_class) { map_.set_key(root_pitch_class); }
    void set_scale(ScaleType scale) { map_.set_scale(scale); }
    void set_off_scale_policy(OffScalePolicy policy) {
        map_.set_off_scale_policy(policy);
    }

    void set_voice_interval(int voice, int steps) {
        if (!valid(voice)) return;
        interval_[voice] = std::clamp(steps, -DiatonicMapT<SampleType>::kIntervalStepsMax,
                                      DiatonicMapT<SampleType>::kIntervalStepsMax);
    }
    int voice_interval(int voice) const { return valid(voice) ? interval_[voice] : 0; }
 
    void set_voice_detune_cents(int voice, double cents) {
        if (!valid(voice) || !std::isfinite(cents)) return;
        detune_cents_[voice] = clamp_finite(cents, -kDetuneMaxCents, kDetuneMaxCents);
    }
    double voice_detune_cents(int voice) const {
        return valid(voice) ? detune_cents_[voice] : 0.0;
    }
 
    void set_voice_level_db(int voice, double db) {
        if (!valid(voice) || !std::isfinite(db)) return;
        level_db_[voice] = clamp_finite(db, kLevelMinDb, kLevelMaxDb);
        level_[voice] = units::db_to_linear(level_db_[voice]);
    }
    double voice_level_db(int voice) const {
        return valid(voice) ? level_db_[voice] : kLevelMinDb;
    }
 
    void set_voice_enabled(int voice, bool on) {
        if (!valid(voice)) return;
        enabled_[voice] = on;
    }
    bool voice_enabled(int voice) const { return valid(voice) && enabled_[voice]; }
 
    void set_dry_level_db(double db) {
        if (!std::isfinite(db)) return;
        dry_db_ = clamp_finite(db, kLevelMinDb, kLevelMaxDb);
        dry_ = units::db_to_linear(dry_db_);
    }
    double dry_level_db() const { return dry_db_; }

    void set_glide_ms(double ms) {
        if (!std::isfinite(ms)) return;
        glide_ms_ = clamp_finite(ms, 0.0, kGlideMsMax);
        for (auto& g : glide_) g.set_time_ms(glide_ms_);
    }
    double glide_ms() const { return glide_ms_; }

    void set_humanize_cents(double cents) {
        if (!std::isfinite(cents)) return;
        humanize_cents_ = clamp_finite(cents, 0.0, kHumanizeMaxCents);
        // DriftT's depth is a percentage of a multiplier near 1, and
        // `cents = 1200·log2(m) ≈ (1200/ln2)·ε` for small ε — so the depth that
        // produces a 1σ of `humanize_cents_` is this, not a free constant.
        const double epsilon_sigma = humanize_cents_ / kCentsPerUnitEpsilon;
        for (auto& d : drift_) d.set_depth_percent(100.0 * epsilon_sigma);
    }
    double humanize_cents() const { return humanize_cents_; }

    void set_crossfade_ms(double ms) {
        if (!std::isfinite(ms)) return;
        crossfade_ms_ = clamp_finite(ms, kCrossfadeMsMin, kCrossfadeMsMax);
        for (auto& s : shifter_) s.set_window_ms(crossfade_ms_);
        update_alignment();
    }
    double crossfade_ms() const { return crossfade_ms_; }
 
    void set_interp(PitchInterp interp) {
        interp_ = interp;
        for (auto& s : shifter_) s.set_interp(interp);
    }
 
    // ── Reporting ─────────────────────────────────────────────────────────

    int latency_samples() const { return latency_; }

    int tracker_latency_samples() const { return tracker_.latency_samples(); }

    int shifter_latency_samples() const {
        return shifter_[0].latency_samples();
    }
 
    double tracked_f0_hz() const { return tracker_.f0_hz(); }
    bool voiced() const { return tracker_.voiced(); }

    const DiatonicMapping& voice_mapping(int voice) const {
        return mapping_[valid(voice) ? voice : 0];
    }

    double voice_cents(int voice) const {
        return valid(voice) ? glide_[voice].value() : 0.0;
    }
    double voice_ratio(int voice) const {
        return valid(voice) ? units::cents_to_ratio(glide_[voice].value()) : 1.0;
    }
    int voice_shift_semitones(int voice) const {
        return valid(voice) ? shift_semitones_[voice] : 0;
    }

    double mute_gain() const { return mute_.value(); }
 
    const YinTrackerT<SampleType>& tracker() const { return tracker_; }
    const DiatonicMapT<SampleType>& diatonic_map() const { return map_; }
 
    // ── Audio ─────────────────────────────────────────────────────────────
 
    SampleType process(SampleType x) {
        // An unprepared instance emits silence rather than faulting. The
        // contract says `prepare` must precede `process`, but `DelayLineT::push`
        // writes its buffer unguarded — so without this, a host that calls
        // `process` early (which real hosts do) gets an out-of-bounds write
        // rather than a diagnosable silence.
        if (latency_ <= 0) return SampleType{0};
 
        if (!std::isfinite(static_cast<double>(x))) {
            discard_history();
            return SampleType{0};
        }
 
        const double in = static_cast<double>(x);
 
        // Stage 1. The tracker sees a DC-blocked copy: a DC offset biases the
        // difference function toward long lags and costs low-note accuracy.
        if (tracker_.process(dc_.process(static_cast<SampleType>(in)))) {
            // A new estimate landed — re-run the mapping. Between hops the
            // mapper output is held and only the glide keeps moving.
            update_targets();
        }
 
        // The alignment line carries the input for BOTH legs: the dry is read
        // at the full reported latency, the shifters at the pre-delay so their
        // own throughput lands them on the same instant.
        align_.push(static_cast<SampleType>(in));
        const double dry = static_cast<double>(align_.read(latency_));
        const double shifter_in = static_cast<double>(align_.read(pre_delay_));
 
        const double gate = mute_.process(tracker_.voiced() ? 1.0 : 0.0);
 
        double wet_sum = 0.0;
        for (int v = 0; v < kMaxVoices; ++v) {
            if (!enabled_[v]) continue;
 
            // Stage 3 control cadence: the cents ramp advances every sample so
            // the glide is smooth between hops, and the ratio follows it.
            double cents = glide_[v].process(target_cents(v));
            if (humanize_cents_ > 0.0)
                cents += units::ratio_to_cents(static_cast<double>(drift_[v].next()));
            apply_ratio(v, cents);
 
            const double wet = static_cast<double>(
                shifter_[static_cast<std::size_t>(v)].process_wet(
                    static_cast<SampleType>(shifter_in)));
            if (gate > kMuteFloor) wet_sum += level_[v] * gate * wet;
        }
 
        return static_cast<SampleType>(dry_ * dry + wet_sum);
    }
 
private:
    static constexpr double kCentsPerUnitEpsilon = 1731.2340490667560888;

    static constexpr std::uint32_t kDriftSeedBase = 0x48524D4Eu;  // 'HRMN'
 
    static bool valid(int voice) { return voice >= 0 && voice < kMaxVoices; }
 
    static double clamp_finite(double v, double lo, double hi) {
        if (!std::isfinite(v)) return lo;
        return std::clamp(v, lo, hi);
    }

    double target_cents(int voice) const {
        return 100.0 * static_cast<double>(shift_semitones_[voice]) +
               detune_cents_[voice];
    }
 
    void apply_ratio(int voice, double cents) {
        shifter_[static_cast<std::size_t>(voice)].set_shift_semitones(cents / 100.0);
    }

    void update_targets() {
        if (!tracker_.voiced()) return;  // hold through unvoiced frames
        for (int v = 0; v < kMaxVoices; ++v) {
            mapping_[v] = map_.map_hz(tracker_.f0_hz(), interval_[v]);
            shift_semitones_[v] = mapping_[v].shift_semitones;
        }
    }
 
    void update_alignment() {
        const int tracker_latency = tracker_.latency_samples();
        const int shifter_latency = shifter_[0].latency_samples();
        latency_ = std::max(tracker_latency, shifter_latency);
        pre_delay_ = latency_ - shifter_latency;
    }
 
    // ── Composed primitives ───────────────────────────────────────────────
    YinTrackerT<SampleType> tracker_{};
    DiatonicMapT<SampleType> map_{};
    PitchShifterT<SampleType> shifter_[kMaxVoices]{};
    DelayLineT<SampleType> align_{};
    DcBlocker<SampleType> dc_{};
    ConstantTimeSlewLimiterT<double> glide_[kMaxVoices]{};
    ConstantTimeSlewLimiterT<double> mute_{};
    DriftT<double> drift_[kMaxVoices]{};
 
    // ── Configuration ─────────────────────────────────────────────────────
    double sample_rate_ = 48000.0;
    double crossfade_ms_ = kCrossfadeMsDefault;
    double glide_ms_ = kGlideMsDefault;
    double humanize_cents_ = kHumanizeDefault;
    double dry_db_ = 0.0;
    double dry_ = 1.0;
    PitchInterp interp_ = PitchInterp::linear;
    int interval_[kMaxVoices] = {kV1IntervalDefault, kV2IntervalDefault};
    double detune_cents_[kMaxVoices] = {0.0, 0.0};
    double level_db_[kMaxVoices] = {0.0, 0.0};
    double level_[kMaxVoices] = {1.0, 1.0};
    bool enabled_[kMaxVoices] = {true, false};
 
    // ── Running state ─────────────────────────────────────────────────────
    DiatonicMapping mapping_[kMaxVoices] = {};
    int shift_semitones_[kMaxVoices] = {0, 0};
    int latency_ = 0;
    int pre_delay_ = 0;
};
 
using HarmonyEngine = HarmonyEngineT<float>;
using HarmonyEngine64 = HarmonyEngineT<double>;
 
}  // namespace pulp::signal