wavetable.h

// This Source Code Form is subject to the terms of the Mozilla Public
// License, v. 2.0. If a copy of the MPL was not distributed with this
// file, You can obtain one at https://mozilla.org/MPL/2.0/.
 
#include <stdlib.h>
#include <math.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
 
struct nec_wt {
    size_t mipmap_count;

    size_t variant_count;

    size_t sample_count;

    float data[];
};
 
// 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;
}
 
// TODO: Attributes, documentation
size_t nec_init_wt(
    void *wt,
    size_t mipmap_count,
    size_t variant_count,
    size_t sample_count
) {
    if (!mipmap_count) mipmap_count = 16;
    if (!variant_count) variant_count = 1;
    if (!sample_count) sample_count = 2048;
    const size_t size = sizeof(struct nec_wt) + sizeof(float[mipmap_count][variant_count][sample_count]);
    if (!wt) return size;
    struct nec_wt *it = wt;
    *it = (struct nec_wt){ mipmap_count, variant_count, sample_count };
    return size;
}
 
[[nodiscard]]
struct nec_wt *nec_alloc_wt(
    size_t mipmap_count,
    size_t variant_count,
    size_t sample_count
) [[clang::allocating]] {
    const size_t size = nec_init_wt(nullptr, mipmap_count, variant_count, sample_count);
    struct nec_wt *wt = malloc(size);
    if (!wt) return nullptr;
    nec_init_wt(wt, mipmap_count, variant_count, sample_count);
    return wt;
}
 
[[nodiscard, gnu::const, gnu::nonnull]]
static inline float nec_wt_get(
    struct nec_wt 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);
 
    // Mipmap indicies
    const auto m = normalized_frequency * (wt->mipmap_count - 1);
    const auto mi = floorf(m); // Integer part of m
    const auto mf = m - mi; // Fractional part of m
    const auto m0 = (size_t)(mi);
    const auto m1 = 1 + m0;
 
    // Variant indicies
    const auto v = variant * (wt->variant_count - 1);
    const auto vi = floorf(v); // Integer part of v
    const auto vf = v - vi; // Fractional part of v
    const auto v0 = (size_t)(vi);
    const auto v1 = 1 + v0;
 
    // Sample indicies
    const auto s = phase * wt->sample_count;
    const auto si = floorf(s); // Integer part of s
    const auto sf = s - si; // Fractional part of s
    const auto s0 = (size_t)(si);
    const auto s1 = (1 + s0) % wt->sample_count;
 
    float (*w)[wt->mipmap_count][wt->variant_count][wt->sample_count] = (typeof(w))(&wt->data);
    return nec_trilerp(
        (*w)[m0][v0][s0],
        (*w)[m0][v0][s1],
        (*w)[m0][v1][s0],
        (*w)[m0][v1][s1],
        (*w)[m1][v0][s0],
        (*w)[m1][v0][s1],
        (*w)[m1][v1][s0],
        (*w)[m1][v1][s1],
        mf, vf, sf
    );
}