wt_lofi.hpp

Go to the documentation of this file.

#pragma once

 
#include <pulp/signal/windowed_sinc_design.hpp>
 
#include <algorithm>
#include <cmath>
#include <cstddef>
#include <span>
#include <utility>
#include <vector>
 
namespace pulp::signal::osc {

class LofiWtOscillator {
public:
    static constexpr int kDefaultBitDepth = 8;
    static constexpr int kFullPrecision = 0;

    static constexpr int kDefaultOversampleFactor = 16;
    static constexpr double kReconstructionStopbandDb = 96.0;

    bool set_tables(std::vector<std::vector<double>> tables,
                    int bit_depth = kDefaultBitDepth) {
        std::size_t length = tables.empty() ? 0 : tables.front().size();
        for (const auto& t : tables) {
            if (t.size() != length) return false; // Mismatched: reject atomically.
        }
 
        bit_depth_ = bit_depth;
        tables_.clear();
        table_length_ = length;
        if (table_length_ > 0) {
            tables_.reserve(tables.size());
            for (auto& src : tables) {
                std::vector<double> table(table_length_, 0.0);
                for (std::size_t i = 0; i < table_length_; ++i)
                    table[i] = quantize(src[i], bit_depth_);
                tables_.push_back(std::move(table));
            }
        }
        current_table_ = std::min<int>(current_table_,
                                       tables_.empty() ? 0
                                       : static_cast<int>(tables_.size()) - 1);
        return true;
    }

    void prepare(double sample_rate) {
        if (sample_rate > 0.0) sample_rate_ = sample_rate;
    }

    void set_scan(double pos) {
        if (tables_.size() < 2) { current_table_ = 0; return; }
        const double clamped = std::clamp(pos, 0.0, 1.0);
        const double scaled = clamped * static_cast<double>(tables_.size() - 1);
        current_table_ = static_cast<int>(std::lround(scaled));
        current_table_ = std::clamp(current_table_, 0,
                                    static_cast<int>(tables_.size()) - 1);
    }

    void set_reconstruction(bool enabled) {
        reconstruction_ = enabled;
        rebuild_filter();
    }

    void set_oversample_factor(int factor) {
        oversample_ = factor > 1 ? factor : 1;
        rebuild_filter();
    }

    void reset() {
        table_phase_ = 0.0;
        std::fill(ring_.begin(), ring_.end(), 0.0);
        ring_pos_ = 0;
    }

    double next(double increment) {
        if (tables_.empty() || table_length_ == 0) return 0.0;
        const double step = increment * static_cast<double>(table_length_);
 
        if (!reconstruction_active()) {
            const double s = sample_zoh();
            advance_phase(step);
            return s;
        }
 
        // Oversample the staircase, push each fine sample through the
        // reconstruction FIR, and keep one decimated sample. Only the retained
        // output's convolution is computed (filter-then-decimate), so the cost is
        // one dot product per host sample regardless of the oversample factor.
        const double sub_step = step / static_cast<double>(oversample_);
        double out = 0.0;
        for (int m = 0; m < oversample_; ++m) {
            push_ring(sample_zoh());
            advance_phase(sub_step);
            if (m == oversample_ - 1) out = convolve_ring();
        }
        return out;
    }
 
    std::size_t table_length() const { return table_length_; }
    std::size_t table_count() const { return tables_.size(); }
    int current_table() const { return current_table_; }
    int bit_depth() const { return bit_depth_; }
    std::span<const double> stored_table(std::size_t i) const {
        return i < tables_.size() ? std::span<const double>(tables_[i])
                                  : std::span<const double>();
    }

    double sample_rate() const { return sample_rate_; }
    bool reconstruction_enabled() const { return reconstruction_; }
    int oversample_factor() const { return oversample_; }
 
private:
    static double quantize(double x, int bits) {
        if (bits <= 0 || bits >= 53) return x;
        const double levels = std::ldexp(1.0, bits - 1) - 1.0;
        return std::round(x * levels) / levels;
    }
 
    bool reconstruction_active() const {
        return reconstruction_ && oversample_ > 1 && !coeffs_.empty();
    }

    double sample_zoh() const {
        const auto& table = tables_[static_cast<std::size_t>(current_table_)];
        std::size_t idx = static_cast<std::size_t>(table_phase_);
        if (idx >= table_length_) idx = table_length_ - 1; // guard fp edge only.
        return table[idx];
    }
 
    void advance_phase(double step) {
        table_phase_ += step;
        const double L = static_cast<double>(table_length_);
        if (table_phase_ >= L || table_phase_ < 0.0)
            table_phase_ -= L * std::floor(table_phase_ / L);
        if (table_phase_ >= L) table_phase_ = 0.0; // fp snap at the boundary.
    }
 
    void push_ring(double x) {
        ring_[ring_pos_] = x;
        ring_pos_ = ring_pos_ + 1 == ring_.size() ? 0 : ring_pos_ + 1;
    }
 
    double convolve_ring() const {
        const std::size_t n = coeffs_.size();
        double acc = 0.0;
        std::size_t idx = ring_pos_ == 0 ? n - 1 : ring_pos_ - 1; // newest sample.
        for (std::size_t j = 0; j < n; ++j) {
            acc += coeffs_[j] * ring_[idx];
            idx = idx == 0 ? n - 1 : idx - 1;
        }
        return acc;
    }

    void rebuild_filter() {
        coeffs_.clear();
        ring_.clear();
        ring_pos_ = 0;
        if (!reconstruction_ || oversample_ <= 1) return;
 
        const double host_nyquist_norm = 0.5 / static_cast<double>(oversample_);
        const double cutoff = host_nyquist_norm;
        const double transition = 0.20 * host_nyquist_norm;
        const double beta = kaiser_beta_for_stopband(kReconstructionStopbandDb);
        const std::size_t taps =
            kaiser_length_for_transition(kReconstructionStopbandDb, transition);
        coeffs_ = design_windowed_sinc(taps, cutoff, beta);
        ring_.assign(coeffs_.size(), 0.0);
    }
 
    std::vector<std::vector<double>> tables_;
    std::size_t table_length_ = 0;
    int bit_depth_ = kDefaultBitDepth;
    int current_table_ = 0;
 
    double sample_rate_ = 48000.0;
    double table_phase_ = 0.0;
 
    bool reconstruction_ = true;
    int oversample_ = kDefaultOversampleFactor;
    std::vector<double> coeffs_;
    std::vector<double> ring_;
    std::size_t ring_pos_ = 0;
};
 
} // namespace pulp::signal::osc