minblep.hpp

Go to the documentation of this file.

#pragma once
#include <dsp/common.hpp>
 
 
namespace rack {
namespace dsp {
 

void minBlepImpulse(int Z, int O, float* output);
 

template <int Z, int O>
struct MinBlep {
    float impulseReordered[O][1 + 2 * Z + 4] = {};
 
    MinBlep() {
        float impulse[2 * Z][O];
        dsp::minBlepImpulse(Z, O, &impulse[0][0]);
 
        // Fill impulseReordered with transposed and pre-subtracted data
        // Storage index = z + 1 (to allow z = -1 access at index 0)
        for (int o = 0; o < O; o++) {
            // z = -1: before impulse starts
            impulseReordered[o][0] = -1.f;
            // z = 0 to 2*Z-1: actual impulse data
            for (int z = 0; z < 2 * Z; z++) {
                impulseReordered[o][1 + z] = impulse[z][o] - 1.f;
            }
            // z = 2*Z to 2*Z+3: after impulse ends (for SIMD padding)
            for (int z = 2 * Z; z < 2 * Z + 4; z++) {
                impulseReordered[o][1 + z] = 0.f;
            }
        }
    }

    const float* getImpulse(int o, int z) const {
        z += o / O;
        o = o % O;
        if (o < 0) {
            z -= 1;
            o += O;
        }
        // assert(0 <= o && o < O);
        // assert(-1 <= z && z < 2 * Z + 4);
        return &impulseReordered[o][z + 1];
    }
 
    template <typename T = float>
    struct State {
        // Double buffer of length 2 * Z
        T buffer[4 * Z] = {};
        int32_t bufferIndex = 0;
    };

    template <typename T>
    void insertDiscontinuity(State<T>& s, float p, T x) const {
        if (!(0 < p && p <= 1))
            return;
 
        // Calculate impulse array index and fractional part
        float subsample = (1 - p) * O;
        int o = (int) subsample;
        float t = subsample - o;
 
        // For each zero crossing, cubic interpolate impulse response
        for (int z = 0; z < 2 * Z; z += 4) {
            simd::float_4 y0 = simd::float_4::load(getImpulse(o - 1, z));
            simd::float_4 y1 = simd::float_4::load(getImpulse(o, z));
            simd::float_4 y2 = simd::float_4::load(getImpulse(o + 1, z));
            simd::float_4 y3 = simd::float_4::load(getImpulse(o + 2, z));
            simd::float_4 y = catmullRomInterpolate(y0, y1, y2, y3, simd::float_4(t));
 
            // Write all 4 samples to contiguous buffer
            for (int zi = 0; zi < 4; zi++) {
                s.buffer[s.bufferIndex + z + zi] += y[zi] * x;
            }
        }
    }

    template <typename T>
    T process(State<T>& s) const {
        T v = s.buffer[s.bufferIndex];
        s.buffer[s.bufferIndex] = T(0);
        s.bufferIndex++;
        if (s.bufferIndex >= 2 * Z) {
            // Move second half of buffer to beginning
            std::memcpy(s.buffer, s.buffer + 2 * Z, 2 * Z * sizeof(T));
            std::memset(s.buffer + 2 * Z, 0, 2 * Z * sizeof(T));
            s.bufferIndex = 0;
        }
        return v;
    }
};
 

template <int Z, int O, typename T = float>
struct MinBlepGenerator {
    MinBlep<Z, O> minBlep;
    typename MinBlep<Z, O>::template State<T> state;

    void insertDiscontinuity(float p, T x) {
        minBlep.insertDiscontinuity(state, p + 1.f, x);
    }
 
    T process() {
        return minBlep.process(state);
    }
};
 
 
} // namespace dsp
} // namespace rack