wavefold_aa.hpp

Go to the documentation of this file.

// sbl/dsp/comp/wavefold_aa.hpp — ADAA Wavefolder (Audio Stack — Compositions)
//
// Second-order antialiased wavefolder for real-time audio.
// Reflects a signal at +/-1.0 boundaries, creating harmonics through folding
// rather than filtering. The core West Coast synthesis technique — a sine
// through a folder is the Buchla 259 sound.
//
// Second-order ADAA (Esqueda/Bilbao, DAFx 2016) integrates over the input
// change twice, giving ~12 dB/octave alias rejection vs ~6 dB for first-order.
// Cost is ~3-4x the naive version — negligible on M7.
//
// For the stateless mathematical fold without antialiasing, see Folder in dsp/folder.hpp.
//
// Depth is applied externally: folder.process(input * depth)
//   depth 1.0 = clean passthrough (for input in [-1,1])
//   depth 2.0 = first fold
//   depth 4.0+ = heavy folding, rich harmonics
//
// Usage:
//   sbl::dsp::comp::WavefoldAA folder;
//   folder.process(buf, frames, 4.0f);              // constant depth
//   folder.process(buf, frames, 0.0f, depth_buf);   // per-sample depth
 
#ifndef SBL_DSP_COMP_WAVEFOLD_AA_HPP_
#define SBL_DSP_COMP_WAVEFOLD_AA_HPP_
 
#include <cmath>
#include <cstdint>
 
#include <sbl/dsp/primitives/folder.hpp>
 
namespace sbl::dsp::comp {
 
/// First antiderivative of the wavefold triangle-wave function (F1).
///
/// The fold function f(x) on its fundamental period [-1, 3):
///   [-1, 1): f(x) = x         (slope +1)
///   [ 1, 3): f(x) = 2 - x    (slope -1)
///
/// F1(x), continuous and periodic (period 4):
///   [-1, 1): F1(x) = x²/2 + 1
///   [ 1, 3): F1(x) = 2x - x²/2
inline float wavefold_antideriv(float x) {
    // Reduce to fundamental period [-1, 3)
    float q = x + 1.0f;
    q = q - 4.0f * floorf(q * 0.25f);  // mod 4, [0, 4)
    float xr = q - 1.0f;               // [-1, 3)
 
    if (xr < 1.0f) {
        return xr * xr * 0.5f + 1.0f;
    } else {
        return 2.0f * xr - xr * xr * 0.5f;
    }
}
 
/// Second antiderivative of the wavefold triangle-wave function (F2).
///
/// F2(x) = integral of F1(x) from 0 to x.
/// On the fundamental period [-1, 3):
///   [-1, 1): F2(x) = x³/6 + x
///   [ 1, 3): F2(x) = x² - x³/6 + 1/3
///
/// F1 has nonzero mean (3/2) over one period, so F2 grows by 6 per period.
/// For arbitrary x: F2(x) = F2_core(x_reduced) + 6 * periods.
inline float wavefold_antideriv2(float x) {
    float q = x + 1.0f;
    float k = floorf(q * 0.25f);       // Number of complete periods
    q = q - 4.0f * k;                  // Reduce to [0, 4)
    float xr = q - 1.0f;               // [-1, 3)
 
    float core;
    if (xr < 1.0f) {
        // integral of (x²/2 + 1) = x³/6 + x
        core = xr * xr * xr * (1.0f / 6.0f) + xr;
    } else {
        // integral of (2x - x²/2) = x² - x³/6 + 1/3
        core = xr * xr - xr * xr * xr * (1.0f / 6.0f) + (1.0f / 3.0f);
    }
 
    return core + 6.0f * k;
}
 
/// Second-order ADAA wavefolder for real-time audio.
///
/// Uses the integrated linear interpolation formulation (Bilbao et al.):
///
///   d1(a, b) = (F2(a) - F2(b)) / (a - b)      when |a - b| > ε
///            = F1(a)                             when |a - b| ≤ ε
///
///   y[n] = 2 * (d1(x[n], x[n-1]) - d1(x[n-1], x[n-2])) / (x[n] - x[n-2])
///
/// The outer denominator spans two samples (x[n] - x[n-2]) for centered
/// differencing.  When well-conditioned, this gives ~12 dB/octave alias
/// rejection — double first-order.
///
/// Near signal peaks where dx is small, the outer double-division suffers
/// catastrophic cancellation.  The fallback is first-order ADAA (single
/// division of F1), which stays numerically stable and still antialiases.
class WavefoldAA {
public:
    /// @note All public methods are ISR-safe — bounded computation, no I/O.
 
    void reset() {
        x1_ = 0.0f;
        x2_ = 0.0f;
        F2_x1_ = wavefold_antideriv2(0.0f);
        F1_x1_ = wavefold_antideriv(0.0f);
        D1_prev_ = wavefold_antideriv(0.0f);   // d1(0, 0) → F1(0) = 1.0
    }
 
    /// Process a single sample (pre-multiplied by depth)
    float process(float x) {
        float F2_x = wavefold_antideriv2(x);
        float F1_x = wavefold_antideriv(x);
 
        // Inner: d1(x[n], x[n-1])
        float dx1 = x - x1_;
        float D1;
        constexpr float EPS = 1e-5f;
        if (fabsf(dx1) > EPS) {
            D1 = (F2_x - F2_x1_) / dx1;
        } else {
            D1 = F1_x;  // Limit: F2'(x) = F1(x)
        }
 
        // Outer: y = 2 * (D1[n] - D1[n-1]) / (x[n] - x[n-2])
        // Use a larger threshold here because the double division amplifies
        // float cancellation errors quadratically.  When ill-conditioned,
        // fall back to first-order ADAA which only divides once.
        float dx2 = x - x2_;
        float y;
        constexpr float OUTER_EPS = 1e-2f;
        if (fabsf(dx2) > OUTER_EPS) {
            y = 2.0f * (D1 - D1_prev_) / dx2;
        } else if (fabsf(dx1) > EPS) {
            // First-order ADAA fallback: single division, always stable
            y = (F1_x - F1_x1_) / dx1;
        } else {
            y = primitives::Folder::fold(x);  // True degenerate: DC input
        }
 
        x2_ = x1_;
        x1_ = x;
        F2_x1_ = F2_x;
        F1_x1_ = F1_x;
        D1_prev_ = D1;
        return y;
    }
 
    /// Process a block in-place with constant or per-sample fold depth
    ///
    /// @param buf       Float audio buffer (modified in-place)
    /// @param frames    Number of samples
    /// @param depth     Constant fold depth (used when depth_buf is nullptr)
    /// @param depth_buf Per-sample depth buffer (nullptr = use constant depth)
    void process(float* buf, uint16_t frames, float depth,
                 const float* depth_buf = nullptr) {
        if (depth_buf != nullptr) {
            for (uint16_t i = 0; i < frames; ++i) {
                buf[i] = process(buf[i] * depth_buf[i]);
            }
        } else {
            for (uint16_t i = 0; i < frames; ++i) {
                buf[i] = process(buf[i] * depth);
            }
        }
    }
 
private:
    float x1_ = 0.0f;                                  // x[n-1]
    float x2_ = 0.0f;                                  // x[n-2]
    float F2_x1_ = wavefold_antideriv2(0.0f);          // Cached F2(x[n-1])
    float F1_x1_ = wavefold_antideriv(0.0f);           // Cached F1(x[n-1])
    float D1_prev_ = wavefold_antideriv(0.0f);         // d1(x[n-1], x[n-2])
};
 
} // namespace sbl::dsp::comp
 
#endif // SBL_DSP_COMP_WAVEFOLD_AA_HPP_