WavetableOscillator.h

// DSPark - Professional Audio DSP Framework
// Copyright (c) 2026 Cristian Moresi - MIT License
 
#pragma once
 
#include "DspMath.h"
#include "AudioSpec.h"
#include "AudioBuffer.h"
#include "Phasor.h"
#include "Interpolation.h"
#include "FFT.h"
 
#include <array>
#include <algorithm>
#include <cassert>
#include <cmath>
#include <vector>
 
namespace dspark {
 
template <FloatType T>
class WavetableOscillator
{
public:
    static constexpr int kTableSize = 2048;
    static constexpr int kTableMask = kTableSize - 1; // Used for ultra-fast wrapping
    static constexpr int kMaxMipLevels = 20;
 
    WavetableOscillator() = default;
 
    void prepare(double sampleRate)
    {
        assert(sampleRate > 0.0);
        if (!(sampleRate > 0.0)) return;
 
        sampleRate_ = sampleRate;
        phasor_.prepare(sampleRate);
        updateMipCutoffs();
        updateMipSelection();
    }
 
    void prepare(const AudioSpec& spec) 
    { 
        prepare(spec.sampleRate); 
    }
 
    // -- Built-in waveform generators (OFFLINE ONLY) ----------------------------
 
    void buildSaw()
    {
        buildFromHarmonics([](int harmonic) -> T {
            return (harmonic % 2 == 0) ? T(-1.0 / harmonic) : T(1.0 / harmonic);
        });
    }
 
    void buildSquare()
    {
        buildFromHarmonics([](int harmonic) -> T {
            return (harmonic % 2 == 0) ? T(0) : T(1.0 / harmonic);
        });
    }
 
    void buildTriangle()
    {
        buildFromHarmonics([](int harmonic) -> T {
            if (harmonic % 2 == 0) return T(0);
            T sign = ((harmonic / 2) % 2 == 0) ? T(1) : T(-1);
            return sign / static_cast<T>(harmonic * harmonic);
        });
    }
 
    void buildSine()
    {
        numMipLevels_ = 1;
        mipData_.assign(kTableSize, T(0)); // Contiguous layout
        levelHarmonics_.assign(1, 1);
 
        for (int i = 0; i < kTableSize; ++i)
        {
            const double phase = twoPi<double> * static_cast<double>(i) / static_cast<double>(kTableSize);
            mipData_[static_cast<size_t>(i)] = static_cast<T>(std::sin(phase));
        }
 
        updateMipCutoffs();
        updateMipSelection();
    }
 
    template <typename HarmonicFunc>
    void buildFromHarmonics(HarmonicFunc harmonicFunc)
    {
        std::vector<double> cosAmp(kMaxHarmonics + 1, 0.0);
        std::vector<double> sinAmp(kMaxHarmonics + 1, 0.0);
        for (int h = 1; h <= kMaxHarmonics; ++h)
            sinAmp[static_cast<size_t>(h)] = static_cast<double>(harmonicFunc(h));
        buildLevels(cosAmp, sinAmp, kMaxHarmonics);
    }
 
    void loadWavetable(const T* data, int size)
    {
        if (!data || size <= 0) return;
 
        // Perform DFT on raw input data to extract true harmonics without
        // interpolation loss. The DC term is discarded (an oscillator should
        // not reproduce offset), and for even sizes the Nyquist bin is
        // excluded too: the 2/N single-sided scaling below would count it
        // twice ((size - 1) / 2 stops one bin short of it).
        int rawHarmonicLimit = (size - 1) / 2;
        int maxHarmonics = std::min(rawHarmonicLimit, kMaxHarmonics);
        if (maxHarmonics < 1) return; // size < 3: no extractable harmonics
 
        // Analysis in double whatever T is (setup-time work).
        std::vector<double> cosCoeffs(static_cast<size_t>(maxHarmonics + 1), 0.0);
        std::vector<double> sinCoeffs(static_cast<size_t>(maxHarmonics + 1), 0.0);
 
        for (int h = 1; h <= maxHarmonics; ++h)
        {
            double sumCos = 0.0, sumSin = 0.0;
            for (int i = 0; i < size; ++i)
            {
                // Integer phase reduction keeps the argument small (exact).
                const auto k = static_cast<long long>(h) * i % size;
                const double phase = twoPi<double> * static_cast<double>(k) / static_cast<double>(size);
                sumCos += static_cast<double>(data[i]) * std::cos(phase);
                sumSin += static_cast<double>(data[i]) * std::sin(phase);
            }
            cosCoeffs[static_cast<size_t>(h)] = sumCos * 2.0 / static_cast<double>(size);
            sinCoeffs[static_cast<size_t>(h)] = sumSin * 2.0 / static_cast<double>(size);
        }
 
        buildLevels(cosCoeffs, sinCoeffs, maxHarmonics);
    }
 
    // -- Playback (REAL-TIME SAFE) ----------------------------------------------
 
    void setFrequency(T frequencyHz) noexcept
    {
        // NaN is ignored, mirroring the internal Phasor: otherwise the pitch
        // would keep the old value while the mip selection went to the
        // dullest level (inconsistent timbre).
        if (frequencyHz != frequencyHz) return;
        const bool changed = frequencyHz != frequency_;
        frequency_ = frequencyHz;
        phasor_.setFrequency(frequencyHz);
        if (changed) updateMipSelection();
    }
 
    [[nodiscard]] T getFrequency() const noexcept { return frequency_; }
 
    [[nodiscard]] inline T getSample() noexcept
    {
        T phase = phasor_.advance();
        return readTable(phase);
    }
 
    void processBlock(T* output, int numSamples) noexcept
    {
        for (int i = 0; i < numSamples; ++i)
            output[i] = getSample();
    }
 
    void generateBlock(AudioBufferView<T> buffer) noexcept
    {
        const int nCh = buffer.getNumChannels();
        const int nS  = buffer.getNumSamples();
        for (int i = 0; i < nS; ++i)
        {
            const T s = getSample();
            for (int ch = 0; ch < nCh; ++ch)
                buffer.getChannel(ch)[i] = s;
        }
    }
 
    void reset(T phase = T(0)) noexcept
    {
        phasor_.reset(phase);
    }
 
private:
 
    [[nodiscard]] inline T readFromLevel(T phase, int level) const noexcept
    {
        T pos = phase * static_cast<T>(kTableSize);
        int i1 = static_cast<int>(pos); // Fast truncation replacing std::floor
        T frac = pos - static_cast<T>(i1);
 
        // Bitwise mask wrapping ensures bounds without modulo operator penalty
        int i0 = (i1 - 1) & kTableMask;
        i1 = i1 & kTableMask;
        int i2 = (i1 + 1) & kTableMask;
        int i3 = (i1 + 2) & kTableMask;
 
        size_t offset = static_cast<size_t>(level * kTableSize);
        const T* table = &mipData_[offset];
 
        return interpolateHermite(table[i0], table[i1], table[i2], table[i3], frac);
    }
 
    [[nodiscard]] inline T readTable(T phase) const noexcept
    {
        if (mipData_.empty()) return T(0);
 
        const T s0 = readFromLevel(phase, selLevel_);
        if (selNext_ == selLevel_ || selWeight_ >= T(1))
            return s0;
        const T s1 = readFromLevel(phase, selNext_);
        return s1 + selWeight_ * (s0 - s1);
    }
 
    void updateMipSelection() noexcept
    {
        selLevel_ = 0;
        selNext_ = 0;
        selWeight_ = T(1);
        if (numMipLevels_ <= 1 || safeFreq_.size() < static_cast<size_t>(numMipLevels_))
            return;
 
        const T f = std::abs(frequency_);
        const int last = numMipLevels_ - 1;
        int level = 0;
        while (level < last && f > safeFreq_[static_cast<size_t>(level)])
            ++level;
 
        selLevel_ = level;
        selNext_ = std::min(level + 1, last);
        if (selNext_ == level) return;
 
        // Bottom of this level's range: the previous level's safe limit, or,
        // for level 0, the point one level-spacing below its own limit.
        const T hi = safeFreq_[static_cast<size_t>(level)];
        const T lo = (level > 0)
            ? safeFreq_[static_cast<size_t>(level - 1)]
            : hi * static_cast<T>(levelHarmonics_[1]) / static_cast<T>(levelHarmonics_[0]);
        const T span = hi - lo;
        selWeight_ = (span > T(0)) ? std::clamp((hi - f) / span, T(0), T(1)) : T(1);
    }
 
    void updateMipCutoffs()
    {
        if (numMipLevels_ <= 0) return;
        safeFreq_.resize(static_cast<size_t>(numMipLevels_));
 
        // Highest harmonic frequency whose alias still lands above 20 kHz
        // (never below Nyquist itself: at low rates there is no free band).
        const double fold = std::max(sampleRate_ * 0.5, sampleRate_ - 20000.0);
        for (int level = 0; level < numMipLevels_; ++level)
            safeFreq_[static_cast<size_t>(level)] = static_cast<T>(
                fold / static_cast<double>(levelHarmonics_[static_cast<size_t>(level)]));
    }
 
    void buildLevels(const std::vector<double>& cosAmp, const std::vector<double>& sinAmp,
                     int available)
    {
        numMipLevels_ = kMaxMipLevels;
        mipData_.assign(static_cast<size_t>(numMipLevels_ * kTableSize), T(0));
        levelHarmonics_.assign(static_cast<size_t>(numMipLevels_), 1);
 
        FFTReal<double> fft(static_cast<size_t>(kTableSize));
        std::vector<double> spec(static_cast<size_t>(kTableSize + 2));
        std::vector<double> cycle(static_cast<size_t>(kTableSize));
        const double half = 0.5 * static_cast<double>(kTableSize);
 
        for (int level = 0; level < numMipLevels_; ++level)
        {
            const int budget = std::min(kLevelHarmonics[static_cast<size_t>(level)], available);
            levelHarmonics_[static_cast<size_t>(level)] = std::max(budget, 1);
 
            // x[n] = sum a_h cos(2 pi h n / N) + b_h sin(2 pi h n / N) has the
            // one-sided spectrum X[h] = (N/2) (a_h - j b_h) under the inverse
            // FFT's 1/N scaling.
            std::fill(spec.begin(), spec.end(), 0.0);
            for (int h = 1; h <= budget; ++h)
            {
                spec[static_cast<size_t>(2 * h)]     =  half * cosAmp[static_cast<size_t>(h)];
                spec[static_cast<size_t>(2 * h + 1)] = -half * sinAmp[static_cast<size_t>(h)];
            }
            fft.inverse(spec.data(), cycle.data());
 
            T* dst = &mipData_[static_cast<size_t>(level * kTableSize)];
            for (int i = 0; i < kTableSize; ++i)
                dst[i] = static_cast<T>(cycle[static_cast<size_t>(i)]);
        }
 
        updateMipCutoffs();
 
        // Normalise ALL levels with ONE global factor (the peak of the richest
        // level). Per-level peak normalisation made the level/timbre jump
        // slightly at every mip crossfade, because the Gibbs overshoot varies
        // with the harmonic count.
        normalizeAllLevelsGlobally();
        updateMipSelection();
    }
 
    void normalizeAllLevelsGlobally()
    {
        T maxVal = T(0);
        for (const T v : mipData_)
            maxVal = std::max(maxVal, std::abs(v));
 
        if (maxVal > T(0))
        {
            const T invMax = T(1) / maxVal;
            for (T& v : mipData_)
                v *= invMax;
        }
    }
 
    static constexpr int kMaxHarmonics = kTableSize / 2 - 1;
 
    static constexpr std::array<int, kMaxMipLevels> kLevelHarmonics = {
        kMaxHarmonics, 724, 512, 362, 256, 181, 128, 90, 64, 45,
        32, 22, 16, 11, 8, 5, 4, 3, 2, 1 };
 
    double sampleRate_ = 48000.0;
    T frequency_ = T(440);
 
    Phasor<T> phasor_;
 
    int numMipLevels_ = 0;
    // Cache-friendly contiguous flat layout (Level0 + Level1 + ...)
    std::vector<T> mipData_;
    std::vector<int> levelHarmonics_;   
    std::vector<T> safeFreq_;           
 
    // Level pair and crossfade weight for the current pitch
    // (updateMipSelection(); read per sample by readTable()).
    int selLevel_ = 0;
    int selNext_ = 0;
    T selWeight_ = T(1);
};
 
} // namespace dspark