polyblep_osc.hpp

Go to the documentation of this file.

// sbl/dsp/widgets/polyblep_osc.hpp — PolyBLEP oscillator (Audio Stack — Widgets)
//
// Band-limited saw, square, triangle, and pulse waveforms using polynomial
// band-limited step (PolyBLEP) correction at discontinuities. Zero flash
// cost — no lookup tables needed.
//
// Saw/Square/Pulse use 2nd-order PolyBLEP at step discontinuities.
// Triangle uses 4th-order Integrated PolyBLEP at slope discontinuities,
// directly band-limiting the naive triangle without a post-filter.
//
// Adapted from Mutable Instruments Plaits oscillator.h and Warps
// oscillator.cc (MIT license, Copyright 2014-2016 Emilie Gillet).
// PolyBLEP math from stmlib/dsp/polyblep.h.
//
// Usage:
//   sbl::dsp::widgets::PolyBlepOsc osc;
//   osc.set_waveform(Waveform::Saw);
//   osc.set_pitch_su(0.5f);
//   osc.process(buf, frames);
//
//   // With per-sample FM modulation:
//   osc.process(buf, frames, fm_buf, 1.0f);  // ±1 octave FM
 
#ifndef SBL_DSP_WIDGETS_POLYBLEP_OSC_HPP_
#define SBL_DSP_WIDGETS_POLYBLEP_OSC_HPP_
 
#include <cstdint>
 
#include <sbl/dsp/math/pitch.hpp>
#include <sbl/dsp/math/clamp.hpp>
#include <sbl/dsp/math/fast_math.hpp>
#include <sbl/dsp/types/fixed.hpp>
 
// Forward declaration for friend access from modules layer
namespace sbl::dsp::modules { template<uint16_t> class MorphOsc; }
 
namespace sbl::dsp::widgets {
 
struct PolyBlepOscTestAccess;
 
enum class Waveform : uint8_t {
    Saw,
    Square,
    Triangle,
    Pulse,
};
 
class PolyBlepOsc {
public:
    /// @note All public methods are ISR-safe — bounded computation, no I/O.
 
    /// Default constructor (Saw waveform).
    PolyBlepOsc() = default;
 
    /// Construct with a specific waveform.
    explicit PolyBlepOsc(Waveform w) : waveform_(w) {}
 
    // ── Configuration ───────────────────────────────────────────
 
    /** @brief Set active waveform (discrete — no _su suffix) */
    void set_waveform(Waveform w) { waveform_ = w; }
 
    // ── Signal-unit inputs ──────────────────────────────────────
 
    /**
     * @brief Set pitch from normalized [0,1] signal unit
     * @param su Unipolar [0,1] — maps to MIDI 24–120 (C1–C9)
     */
    void set_pitch_su(float su) {
        if (su < 0.0f) su = 0.0f;
        if (su > 1.0f) su = 1.0f;
        float note = kMinNote + su * kNoteRange;
        set_note(note);
    }
 
    /**
     * @brief Set pulse width from normalized [0,1] signal unit
     * @param su Unipolar [0,1] — linear map to 0.05–0.95
     */
    void set_pulse_width_su(float su) {
        if (su < 0.0f) su = 0.0f;
        if (su > 1.0f) su = 1.0f;
        set_pulse_width(0.05f + su * 0.9f);
    }
 
    /**
     * @brief Set output amplitude from normalized [0,1] signal unit
     * @param su Unipolar [0,1] — linear map to amplitude 0–1
     */
    void set_amplitude_su(float su) {
        if (su < 0.0f) su = 0.0f;
        if (su > 1.0f) su = 1.0f;
        set_amplitude(su);
    }
 
    /**
     * @brief Generate audio samples, optionally with per-sample FM modulation
     *
     * When fm_mod is provided, per-sample exponential FM is applied around
     * the stored frequency (set via set_pitch_su / set_frequency):
     *   freq = frequency_ * 2^(fm_mod[i] * fm_depth)
     *
     * @param out      Output buffer (float, [-1.0, 1.0])
     * @param frames   Number of samples
     * @param fm_mod   FM modulation signal [-1,1] (nullptr = no modulation)
     * @param fm_depth FM depth in octaves (e.g., 1.0 = ±1 octave)
     */
    void process(float* out, uint16_t frames,
                 const float* fm_mod = nullptr, float fm_depth = 0.0f) {
        if (fm_mod != nullptr) {
            // Per-sample FM — render one sample at a time
            float base_freq = frequency_;
            for (uint16_t i = 0; i < frames; ++i) {
                float freq = base_freq * math::fast_exp2f(fm_mod[i] * fm_depth);
                if (freq < 0.0f) freq = 0.0f;
                frequency_ = freq;
                render_float(out + i, 1);
            }
            frequency_ = base_freq;  // Restore base frequency
            if (amplitude_ != 1.0f) {
                for (uint16_t i = 0; i < frames; ++i) {
                    out[i] *= amplitude_;
                }
            }
        } else {
            // No modulation — block render
            uint16_t pos = 0;
            while (pos < frames) {
                uint16_t n = frames - pos;
                if (n > MAX_BLOCK) n = MAX_BLOCK;
                render_float(out + pos, n);
                if (amplitude_ != 1.0f) {
                    for (uint16_t i = 0; i < n; ++i) {
                        out[pos + i] *= amplitude_;
                    }
                }
                pos += n;
            }
        }
    }
 
    /** @brief Hard sync — reset phase to zero */
    void sync() {
        phase_ = 0.0f;
        next_sample_ = 0.0f;
        high_ = true;
    }
 
    /** @brief Reset phase to a specific value (0.0–1.0) */
    void sync(float phase) {
        phase_ = phase;
        next_sample_ = 0.0f;
    }
 
    /** @brief Current phase as uint32_t */
    uint32_t phase() const {
        return static_cast<uint32_t>(phase_ * 4294967296.0f);
    }
 
    /** @brief Current normalized frequency (freq_hz / sample_rate) */
    float frequency() const { return frequency_; }
 
private:
    friend struct PolyBlepOscTestAccess;
    template<uint16_t> friend class sbl::dsp::modules::MorphOsc;
 
    // ── Pitch mapping constants ─────────────────────────────────
    static constexpr float kMinNote  = 24.0f;   // C1
    static constexpr float kMaxNote  = 120.0f;  // C9
    static constexpr float kNoteRange = kMaxNote - kMinNote;  // 96 semitones
 
    // ── Engineering-unit setters (private — use _su inputs) ─────
 
    void set_frequency(float freq_hz) {
        frequency_ = freq_hz / types::SAMPLE_RATE_F;
    }
 
    void set_note(float midi_note) {
        frequency_ = math::note_to_frequency(midi_note) / types::SAMPLE_RATE_F;
    }
 
    void set_pulse_width(float pw) {
        pw_ = (pw < 0.05f) ? 0.05f : (pw > 0.95f) ? 0.95f : pw;
    }
 
    void set_amplitude(float amp) { amplitude_ = amp; }
 
    // ── State ───────────────────────────────────────────────────
    static constexpr uint16_t MAX_BLOCK = 48;
 
    float phase_ = 0.0f;
    float frequency_ = 0.0f;
    float pw_ = 0.5f;
    float next_sample_ = 0.0f;
    bool high_ = true;
    Waveform waveform_ = Waveform::Saw;
    float amplitude_ = 1.0f;
 
    // -------------------------------------------------------------------------
    // Dispatch to waveform-specific float render
    // -------------------------------------------------------------------------
    void render_float(float* out, uint16_t frames) {
        switch (waveform_) {
            case Waveform::Saw:      render_saw_f(out, frames); break;
            case Waveform::Square:   render_square_f(out, frames, 0.5f); break;
            case Waveform::Triangle: render_triangle_f(out, frames); break;
            case Waveform::Pulse:    render_square_f(out, frames, pw_); break;
        }
    }
 
    // -------------------------------------------------------------------------
    // Saw — single discontinuity at phase wrap (Warps lines 141-155)
    // -------------------------------------------------------------------------
    void render_saw_f(float* out, uint16_t frames) {
 
        float phase = phase_; // TODO why aren't we using phase accumulator here?
        float next = next_sample_;
        float freq = frequency_;
        float gain = hf_gain(freq);
 
        for (uint16_t i = 0; i < frames; ++i) {
            float this_sample = next;
            next = 0.0f;
 
            phase += freq;
            if (phase >= 1.0f) {
                phase -= 1.0f;
                float t = phase / freq;
                this_sample -= this_blep(t);
                next -= next_blep(t);
            }
            next += phase;
            this_sample = 2.0f * this_sample - 1.0f;
 
            out[i] = math::clamp(this_sample, -1.0f, 1.0f) * gain;
        }
        phase_ = phase;
        next_sample_ = next;
    }
 
    // -------------------------------------------------------------------------
    // Square / Pulse — two edges per cycle
    // -------------------------------------------------------------------------
    void render_square_f(float* out, uint16_t frames, float pw) {
 
        float phase = phase_;
        float next = next_sample_;
        float freq = frequency_;
        bool high = high_;
        float gain = hf_gain(freq);
 
        float fall_at = pw;
 
        for (uint16_t i = 0; i < frames; ++i) {
            float this_sample = next;
            next = 0.0f;
 
            phase += freq;
 
            // Falling edge at duty cycle boundary
            if (high && phase >= fall_at) {
                float t = (phase - fall_at) / freq;
                this_sample -= this_blep(t);
                next -= next_blep(t);
                high = false;
            }
            // Rising edge at phase wrap
            if (phase >= 1.0f) {
                phase -= 1.0f;
                float t = phase / freq;
                this_sample += this_blep(t);
                next += next_blep(t);
                high = true;
            }
 
            next += high ? 1.0f : 0.0f;
            this_sample = 2.0f * this_sample - 1.0f;
 
            out[i] = math::clamp(this_sample, -1.0f, 1.0f) * gain;
        }
        phase_ = phase;
        next_sample_ = next;
        high_ = high;
    }
 
    // -------------------------------------------------------------------------
    // Triangle — Integrated PolyBLEP on slope discontinuities
    // -------------------------------------------------------------------------
    void render_triangle_f(float* out, uint16_t frames) {
 
        float phase = phase_;
        float next = next_sample_;
        float freq = frequency_;
        bool high = high_;
        float gain = hf_gain(freq);
 
        constexpr float slope_up = 2.0f;
        constexpr float slope_down = 2.0f;
        const float discontinuity = (slope_up + slope_down) * freq;
 
        for (uint16_t i = 0; i < frames; ++i) {
            float this_sample = next;
            next = 0.0f;
 
            phase += freq;
 
            // Slope change at midpoint (rising → falling)
            if (high ^ (phase < 0.5f)) {
                float t = (phase - 0.5f) / freq;
                this_sample -= this_integrated_blep(t) * discontinuity;
                next -= next_integrated_blep(t) * discontinuity;
                high = phase < 0.5f;
            }
            // Slope change at wrap (falling → rising)
            if (phase >= 1.0f) {
                phase -= 1.0f;
                float t = phase / freq;
                this_sample += this_integrated_blep(t) * discontinuity;
                next += next_integrated_blep(t) * discontinuity;
                high = true;
            }
 
            // Naive triangle: rises [0, 0.5) then falls [0.5, 1.0)
            next += high
                ? phase * slope_up
                : 1.0f - (phase - 0.5f) * slope_down;
 
            out[i] = math::clamp(2.0f * this_sample - 1.0f, -1.0f, 1.0f) * gain;
        }
        phase_ = phase;
        next_sample_ = next;
        high_ = high;
    }
 
    // -------------------------------------------------------------------------
    // High-frequency gain fadeout (FDP-049)
    //
    // Smoothly attenuate output as normalized frequency approaches Nyquist.
    // Instead of hard-clamping frequency (which creates FM sweep discontinuities)
    // or allowing incoherent above-Nyquist output, the oscillator tracks its
    // modulator faithfully — it just gets quiet where aliasing lives.
    //
    // Smoothstep curve (3t² - 2t³): zero derivative at both boundaries, so
    // entering and leaving the fade zone creates no spectral splatter from
    // abrupt amplitude modulation. Analog VCOs fade naturally into silence
    // above the audio band — this is the digital equivalent.
    // -------------------------------------------------------------------------
    static float hf_gain(float freq) {
        constexpr float FADE_START = 0.40f;  // 19.2 kHz at 48 kHz
        constexpr float FADE_END   = 0.48f;  // 23.0 kHz at 48 kHz
        if (freq <= FADE_START) return 1.0f;
        if (freq >= FADE_END)   return 0.0f;
        // Normalize to [0, 1] within fade range
        constexpr float INV_RANGE = 1.0f / (FADE_END - FADE_START);  // 12.5
        float t = (freq - FADE_START) * INV_RANGE;
        // Smoothstep: 1 - (3t² - 2t³) = 1 - t²(3 - 2t)
        return 1.0f - t * t * (3.0f - 2.0f * t);
    }
 
    // -------------------------------------------------------------------------
    // PolyBLEP correction functions (private implementation)
    //
    // Polynomial band-limited step corrections at discontinuities.
    // Adapted from MI stmlib/dsp/polyblep.h (MIT, Emilie Gillet).
    // -------------------------------------------------------------------------
 
    /// 2nd-order step correction — current sample
    static float this_blep(float t) { return 0.5f * t * t; }
 
    /// 2nd-order step correction — next sample
    static float next_blep(float t) { t = 1.0f - t; return -0.5f * t * t; }
 
    /// 4th-order slope correction — next sample
    static float next_integrated_blep(float t) {
        const float t1 = 0.5f * t;
        const float t2 = t1 * t1;
        const float t4 = t2 * t2;
        return 0.1875f - t1 + 1.5f * t2 - t4;
    }
 
    /// 4th-order slope correction — current sample
    static float this_integrated_blep(float t) {
        return next_integrated_blep(1.0f - t);
    }
};
 
} // namespace sbl::dsp::widgets
 
#endif // SBL_DSP_WIDGETS_POLYBLEP_OSC_HPP_