wavetable.h

#include <stdlib.h>
#include "vendor/prelude/hint.h"
 
// TODO: Anti-aliasing filter. Should probably be in its own header
// since it's not strictly for use in the wavetables and should also be
// useful for real-time signals. Putting the reading links here first
// though since it will be used to filter the mipmaps.
// https://ldesoras.fr/doc/articles/resampler-en.pdf
// https://www.kvraudio.com/forum/viewtopic.php?t=458831
// https://www.earlevel.com/main/2020/01/04/further-thoughts-on-wave-table-oscillators
// https://community.native-instruments.com/discussion/411/anti-aliasing-strategies-for-a-phase-distortion-oscillator
 
typedef struct nec_wavetable {
    size_t mipmap_count; // Spectral axis (a mipmap of variants)
    size_t variant_count; // Feature axis (a variant of cycles)
    size_t sample_count; // Temporal axis (a cycle of frames)
    float data[]; // float[mipmap_count][variant_count][sample_count]
} nec_wt_t;
 
// TODO: Why isn't this autovectorized?
[[nodiscard, gnu::const]]
static inline float nec_trilerp(
    float c000,
    float c001,
    float c010,
    float c011,
    float c100,
    float c101,
    float c110,
    float c111,
    float x, float y, float z
) [[unsequenced]] {
    assume(0 <= x && x <= 1);
    assume(0 <= y && y <= 1);
    assume(0 <= z && z <= 1);
    const auto z1 = 1.f - z;
    const auto c00 = c000 * z1 + c001 * z;
    const auto c01 = c010 * z1 + c011 * z;
    const auto c10 = c100 * z1 + c101 * z;
    const auto c11 = c110 * z1 + c111 * z;
    const auto y1 = 1.f - y;
    const auto c0 = c00 * y1 + c01 * y;
    const auto c1 = c10 * y1 + c11 * y;
    return c0 * (1.f - x) + c1 * x;
}
 
[[nodiscard]]
nec_wt_t *nec_make_wt(
    size_t mipmap_count,
    size_t variant_count,
    size_t sample_count
) [[clang::allocating]] {
    if (!mipmap_count) { mipmap_count = 16; }
    if (!variant_count) { variant_count = 1; }
    if (!sample_count) { sample_count = 2048; }
    const size_t data_len = mipmap_count * variant_count * sample_count;
    return calloc(1, sizeof(nec_wt_t) + sizeof(float) * data_len);
}
 
[[nodiscard, gnu::const, gnu::nonnull]]
static inline float nec_wt_get(
    nec_wt_t wt[const restrict static 1],
    float normalized_frequency,
    float variant,
    float phase
) [[unsequenced]] {
    assume(0 <= normalized_frequency && normalized_frequency < 1);
    assume(0 <= variant && variant <= 1);
    assume(0 <= phase && phase < 1);
 
    const auto samples_per_mipmap = wt->sample_count * wt->variant_count;
    const auto mipmap_idxr = normalized_frequency * (wt->mipmap_count - 1);
    const auto mipmap_idxi = ceilf(mipmap_idxr);
    const auto mipmap_idxf = mipmap_idxr - mipmap_idxi;
    const auto mipmap_idx1 = (size_t)(mipmap_idxi);
    const auto mipmap_idx0 = mipmap_idx1 - 1;
    const auto m0 = mipmap_idx0 * samples_per_mipmap;
    const auto m1 = mipmap_idx1 * samples_per_mipmap;
 
    const auto variant_idxr = variant * (wt->variant_count - 1);
    const auto variant_idxi = floorf(variant_idxr);
    const auto variant_idxf = variant_idxr - variant_idxi;
    const auto variant_idx0 = (size_t)(variant_idxi);
    const auto variant_idx1 = 1 + variant_idx0;
    const auto v0 = variant_idx0 * wt->sample_count;
    const auto v1 = variant_idx1 * wt->sample_count;
 
    const auto sample_idxr = phase * wt->sample_count;
    const auto sample_idxi = floorf(sample_idxr);
    const auto sample_idxf = sample_idxr - sample_idxi;
    const auto sample_idx0 = (size_t)(sample_idxi);
    const auto sample_idx1 = (1 + sample_idx0) % wt->sample_count;
 
    return nec_trilerp(
        wt->data[m0 + v0 + sample_idx0],
        wt->data[m0 + v0 + sample_idx1],
        wt->data[m0 + v1 + sample_idx0],
        wt->data[m0 + v1 + sample_idx1],
        wt->data[m1 + v0 + sample_idx0],
        wt->data[m1 + v0 + sample_idx1],
        wt->data[m1 + v1 + sample_idx0],
        wt->data[m1 + v1 + sample_idx1]
    );
}