DWG (Digital Waveguide) inclusion in chowdsp_utils ?

[ghost]

I don't know if you've already have this idea ? Long time ago, I've work on my own implementation of a waveguide (bi-directional, / multi-dimensionnal) that is an extension of the Julis Orion Smith historical implementation. It's not very modern but it's efficient and extensible. The code is not framework ready but that can be a starting point if you want to integration this kind of feature in your framework.

The main class is waveguide.hpp and consist of a multi-dimensionnal bi-directional two-rails (delaylines) with random access and two terminations :

#pragma once

#include <algorithm>
#include <vector>

#include <delay.hpp>
#include <lagrange.hpp>

namespace dwg {

template <typename T, size_t dims = 1, typename interpolation = details::lagrange<T, 5>>
class waveguide
{
public:
    waveguide(double size)
        : upper(size * 0.5)
        , lower(size * 0.5)
    {}

    constexpr size_t size()
    {
        return upper.size();
    }

    constexpr void set(double position, T value, size_t dim = 0)
    {
        auto up = upper.size(dim) * position;
        auto lp = lower.size(dim) - up;
        upper.set(up, value * 0.5, dim);
        lower.set(lp, value * 0.5, dim);
    }

    constexpr T get(double position, size_t dim = 0)
    {
        auto up = upper.size(dim) * position;
        auto lp = lower.size(dim) - up;
        return dc(upper.get(up, dim)
                + lower.get(lp, dim));
    }

    constexpr void fill(std::vector<T> values, size_t dim = 0)
    {
        auto normalized = std::transform(values.begin(), values.end(), [](auto& x) { return x * 0.5; });
        upper.fill(normalized, dim);
        std::reverse(values.begin(), values.end());
        lower.fill(normalized, dim);
    }

    constexpr void move()
    {
        #pragma simd
        for (auto dim = 0; dim < dims; ++dim)
        {
            auto up = upper.out(dim);
            auto lo = lower.out(dim);
            upper.in( LT ( lo, dim ), dim);
            lower.in( RT ( up, dim ), dim);
        }
    }

protected:
    details::delay<T, dims, interpolation, details::forward>  upper; // fractional delay
    details::delay<T, dims, interpolation, details::backward> lower; // fractional delay

private:
    struct dcremover
    {
        inline T operator()(T value)
        {
            static T n = 1;
            static T mean = 0;
            mean += value;
            return value - (mean / n++);
        }
    };
    dcremover dc;

private:
    virtual T LT (T value, size_t dim) = 0; //  left termination
    virtual T RT (T value, size_t dim) = 0; // right termination
};

} // namespace dwg

the delay.hpp is a simple old-school shared aligned memory with moving-head-pointer, the code include the well-know crossfade-trick that allow fast & click-free resize of the delay as experimental feature (fixed period / samplerate in constructor).

#pragma once

#include <crossfade.hpp>
#include <lagrange.hpp>
#include <memory.hpp>
#include <pointer.hpp>

namespace dwg {

namespace details {

template <typename T, size_t dims, typename interpolation = lagrange<T, 5>, size_t direction = forward>
class delay : public interpolation, public crossfade<T, dims>
{
public:
    using Ptr = pointer<T, interpolation, direction>; // read/write head (oldschool moving-pointer trick)

public:
    delay(double size, size_t max_size = 1024)
        : crossfade<T>(std::round(44.1 * 28)) // crossfade of 28 ms @ 44100 hz
    {
        for (auto dim = 0; dim < dims; ++dim)
        {
            data[dim] = head[dim] = (T*)mem::aligned_alloc(1024, max_size * sizeof(T));
            std::memset(data[dim], 0, max_size * sizeof(T));

            cur_ptr[dim] = new Ptr(data[dim], size, max_size, order);
            nxt_ptr[dim] = new Ptr(data[dim], size, max_size, order);
        }
    }

    ~delay()
    {
        for (auto dim = 0; dim < dims; ++dim)
        {
            delete cur_ptr[dim];
            delete nxt_ptr[dim];

            mem::free(data[dim]);
        }
    }

    constexpr size_t size(size_t dim = 0)
    {
        return cur_ptr[dim]->size;
    }

    constexpr void resize(double new_size, size_t dim = 0, bool smooth = true)
    {
        if (smooth)
        {
            nxt_ptr[dim]->resize(new_size);
            crossfade<T>::reset(dim);
        }
        else
            cur_ptr[dim]->resize(new_size);
    }

    constexpr void in(T value, size_t dim = 0)
    {
        if (crossfade<T>::running())
        {
            interpolation::deinterpolate(cur_ptr[dim], head[dim], 0, value);
            interpolation::deinterpolate(nxt_ptr[dim], head[dim], 0, value);
        }
        else
        {
            interpolation::deinterpolate(cur_ptr[dim], head[dim], 0, value);
        }
    }

    constexpr T out(size_t dim = 0)
    {
        head[dim] = cur_ptr[dim]->move(head[dim]);

        return (crossfade<T>::running())
            ? crossfade<T>::process
             (interpolation::interpolate(cur_ptr[dim], head[dim], cur_ptr[dim]->size),
              interpolation::interpolate(nxt_ptr[dim], head[dim], nxt_ptr[dim]->size))
            : interpolation::interpolate(cur_ptr[dim], head[dim], cur_ptr[dim]->size);
    }

    constexpr void set(double position, T value, size_t dim = 0)
    {
        if (crossfade<T>::running())
        {
            interpolation::deinterpolate(cur_ptr[dim], head[dim], cur_ptr[dim]->size * position, value);
            interpolation::deinterpolate(nxt_ptr[dim], head[dim], nxt_ptr[dim]->size * position, value);
        }
        else
        {
            interpolation::deinterpolate(cur_ptr[dim], head[dim], cur_ptr[dim]->size * position, value);
        }
    }

    constexpr T get(double position, size_t dim = 0)
    {
        return (crossfade<T>::running(dim))
            ? crossfade<T>::process
             (interpolation::interpolate(cur_ptr[dim], head[dim], cur_ptr[dim]->size * position),
              interpolation::interpolate(nxt_ptr[dim], head[dim], nxt_ptr[dim]->size * position), dim)
            : interpolation::interpolate(cur_ptr[dim], head[dim], cur_ptr[dim]->size * position);
    }

    constexpr void fill(std::vector<T> const& values, size_t dim = 0)
    {
        int size = values.size();
        for (auto i = 0; i < size; ++i)
            cur_ptr[dim]->set(head[dim], i, values[i]);
    }

private:
    T* data[dims]; // shared aligned memory
    T* head[dims]; // shared pointer

    // Legato trick (smooth resize)
    Ptr* cur_ptr[dims]; // current note
    Ptr* nxt_ptr[dims]; // next note to crossfade

    // alternate crossfade trick (see. crossfade.hpp)
    void complete(size_t dim = 0) override
    {
        // just swap pointers
        auto temp    = cur_ptr[dim];
        cur_ptr[dim] = nxt_ptr[dim];
        nxt_ptr[dim] = temp;
    }

    constexpr void move(int dim = -1)
    {
        if constexpr (dim >= 0)
        {
            head[dim] = cur_ptr[dim]->move(head[dim]);
        }
        else
        {
            #pragma simd
            for (auto d = 0; d < dims; ++d)
                head[d] = cur_ptr[d]->move(head[d]);
        }
    }
};

} // namespace details

} // namespace dwg

the lagrange.hpp is a underrated technique of interpolation/deinterpolation that allow the random access over the delay line at any fractional point. This is a N-order Lagrange Interpolation based on faust code with some trick from some scientific litterature. Default is 5th order initialized by the template parameters.

#pragma once

#include <pointer.hpp>

namespace dwg {

namespace details {

template <typename T, size_t N>
class lagrange
{
public:
    using Ptr = pointer<T, lagrange<T, N>>;

public:
    static inline T fractional(T size) 
    {
        auto o = (double(N) - 1.00001) * 0.5; // ~ center FIR interpolator
        auto dmo = size - o; // assumed >=0 [size > (N-1)/2]
        return o + (dmo - std::floor(dmo)); // fractional part
    }

    static inline void coefficients(double d, std::vector<T>& h)
    {
        h.resize(N+1);
        std::fill(h.begin(), h.end(), 1);
        std::vector<size_t> n(N); 
        std::iota(n.begin(), n.end(), 0);
        for (auto k = 0; k < N+1; ++k)
        {
            auto idx = std::find_if(n.begin(), n.end(), [k](auto x) { return x != k; });
            for (auto* i : idx)
                h[i] = h[i] * (d - k) / (n[i] - k);
        }
    }

    static inline T interpolate(Ptr* ptr, T* head, int index)
    {
        auto out = T(0);
        for (auto i = 0; i < N+1; ++i)
            out += ptr->get(head, index-i) * ptr->fir(i);
        return out;
    }

    static inline void deinterpolate(Ptr* ptr, T* head, int index, T value)
    {
        for (auto i = 0; i < N+1; ++i)
            ptr->add(head, index-(N-i), value * ptr->fir(i));
    }
};

} // namespace details

} // namespace dwg

the pointer.hpp is a naive implementation of a bidirectional moving pointer allowing to selected the 'true' propagation direction. This is the most efficient way to simulate signal move inside a delay-line and this is a special case (bidirectional) for me to keep the code logical and coherent with the theory.

#pragma once

#include <vector>

namespace dwg {

namespace details {

enum { forward, backward };

template <typename T, typename interpolation, size_t direction = forward>
class pointer
{
public:
    pointer(T* shared_memory, double size, size_t max_size)
        : data(shared_memory)
        , ends(shared_memory + max_size-1)
        , full(max_size)
    {
        resize(size);
    }

    void resize(double size)
    {
        this->frac = interpolation::fractional(size);
        interpolation::coefficients(size, coeffs);
        this->size = int(size - frac);
    }

    constexpr void in (T* head, T value)
    {
        *(head) = value;
    }

    constexpr T* move (T* head)
    {
        if constexpr (direction == backward)
            return mod(--head);
        else
            return mod(++head);
    }

    constexpr T out (T* head)
    {
        if constexpr (direction == backward)
            return *(mod(head + size-1));
        else
            return *(mod(head - size-1));
    }

    constexpr T get (T* head, size_t index)
    {
        if constexpr (direction == backward)
            return *(mod(head + index));
        else
            return *(mod(head - index));
    }

    constexpr void set (T* head, size_t index, T value)
    {
        if constexpr (direction == backward)
            *(mod(head + index)) = value;
        else
            *(mod(head - index)) = value;
    }

    constexpr void add (T* head, size_t index, T value)
    {
        if constexpr (direction == backward)
            *(mod(head + index)) += value;
        else
            *(mod(head - index)) += value;
    }

    constexpr void sub (T* head, size_t index, T value)
    {
        if constexpr (direction == backward)
            *(mod(head + index)) -= value;
        else
            *(mod(head - index)) -= value;
    }

    constexpr void mul (T* head, int index, T value)
    {
        if constexpr (direction == backward)
            *(mod(head + index)) *= value;
        else
            *(mod(head - index)) *= value;
    }

private:
    T* data; // shared memory block (left bound)
    T* ends; // shared memory block (right bound)

    size_t full; // memory block size

    inline T* mod (T* ptr)
    {
        while (ptr <  data)
               ptr += full;
        while (ptr >  ends)
               ptr -= full;
        return ptr;
    }

public:
    size_t size; // integer part
    double frac; // fractional part

    size_t N; // FIR filter order
    std::vector<T> coeffs; // FIR coefficients

    inline T fir (size_t index)
    {
        return coeffs[index];
    }
};

} // namespace details

} // namespace dwg

the crossfade.hpp is using to state-crossfader moving around a lookup-table of a S-shape curve (sigmoid) and is used for the 'legato' trick.

#pragma once

#include <vector>

namespace dwg {

namespace details {

template <typename T, size_t dims = 1>
class crossfade
{
public:
    crossfade(size_t period)
        : counter(0)
    {
        auto sigmoid = [](T x, T a = 0.787) -> T
        {
            constexpr T epsilon = 0.0001;
            constexpr T min_param_a = 0.0 + epsilon;
            constexpr T max_param_a = 1.0 - epsilon;
            a = std::max(min_param_a, std::min(max_param_a, a));
            a = (T(1) / (T(1) - a) - T(1));

            // "Logistic" Sigmoid function
            T A = 1.0 / (1.0 + exp(0 - ((x - 0.5)*a*2.0)));
            T B = 1.0 / (1.0 + exp(a));
            T C = 1.0 / (1.0 + exp(0 - a));
            T y = (A - B) / (C - B);

            return y;
        };

        // optimized S-shaped curve
        curve.resize(period+1);
        std::fill(curve.begin(), curve.end(), 1);
        auto ratio = 1.0 / period;
        for (auto i = 0; i < period; ++i)
            curve[i] = sigmoid(ratio * i);
    }

    inline void reset(size_t dim = 0)
    {
        counter[dim] = curve.size()-1;
    }

    inline bool running(size_t dim = 0)
    {
        return counter[dim] > 0;
    }

    inline T process(T A, T B, size_t dim = 0)
    {
        auto ratio = curve[step(dim)];
        return (A *        ratio)
             + (B * (1.0 - ratio));
    }

    virtual void complete(size_t dim = 0) = 0;

private:
    std::vector<T> curve; // S-shaped curve
    int counter[dims]; // crossfade evolution counter

    inline size_t step(size_t dim = 0)
    {
        if (--counter[dim] < 0)
            complete();

        return (counter[dim] < 0) ? 0 : counter[dim];
    }
};

} // namespace details

} // namespace dwg

To finish a memory.hpp helpers is here to workaround the aligned_memory for microsoft visual studio compiler. Don't know if it's always pertinent since latest VS2019 version includes the std::aligned_memory method.

#pragma once

#include <memory>
#include <type_traits>

// Since Visual Studio C++17 doesn't implement std::aligned_alloc
// just hack the std namespace for the Microsoft equivalent methods

#ifdef _MSC_VER 
namespace mem
{
    inline void* aligned_alloc(std::size_t alignment, std::size_t size)
    {
        return _aligned_malloc(size, alignment);
    }

    inline void free(void* ptr)
    {
        _aligned_free(ptr);
    }
}
#else
namespace mem
{
    inline void* aligned_alloc(std::size_t alignment, std::size_t size)
    {
        return std::aligned_malloc(alignment, size);
    }

    inline void free(void* ptr)
    {
        std::free(ptr);
    }
}
#endif // MSVC

Now, how to use the classes ? The most simples example is the JOS digital waveguide example (string.hpp) that is a simple bidirectional waveguide with a simple one-pole lowpass at right termination.

#pragma once

#include <waveguide.hpp>

namespace dwg {

template <typename T, size_t dims = 1>
class string : public waveguide<T, dims>
{
public:
   string(double size) : waveguide<T, dims>(size)
   {
      // clear lowpass filter states
      for (auto d = 0; d < dims; ++d)
         state[d] = T(0);
   }

private:
   T state[dims]; // filter memory

   inline T LT (T value, size_t) override
   {
      return -value;
   }

   inline T RT (T value, size_t dim) override
   {
      // one-pole lowpass filter with feedback
      auto output = (state[dim] * 0.5) + (value * 0.5);
      state[dim] = output; // z-1
      return -output;
   }
};

} // namespace dwg

And this is a very basic code to use it (main.cpp):

#include <iostream>

#include <string.hpp>

template <typename T>
struct quadratic
{
    quadratic(T x1, T y1, T x2, T y2, T x3, T y3)
        : x1(x1), y1(y1)
        , x2(x2), y2(y2)
        , x3(x3), y3(y3)
    {}

    inline std::tuple<T, T> operator()(T perc)
    {
        auto xa = pt(x1, x2, perc);
        auto ya = pt(y1, y2, perc);
        auto xb = pt(x2, x3, perc);
        auto yb = pt(y2, y3, perc);

        auto x = pt(xa, xb, perc);
        auto y = pt(ya, yb, perc);

        return { x, y };
    }

private:
    T x1; T y1; // P0
    T x2; T y2; // P1
    T x3; T y3; // P2

    inline T pt(T n1, T n2, T perc)
    {
        auto diff = n2 - n1;
        return n1 + (diff * perc);
    }
};

template <typename T>
static std::vector<T> get_initial_displacement(size_t rail_length, double pick, double amplitude = 1.0)
{
    auto quad = quadratic<T>(0.f, 0.f, pick, 2.f, 1.f, 0.f);
    float ratio = 1.f / rail_length;
    std::vector<T> initial_shape(rail_length);
    for (auto i = 0; i < rail_length; ++i)
        initial_shape[i] = std::get<1>(quad(ratio * i) * amplitude);
    return initial_shape;
}

inline int samples_for_ms(double ms, double fs = 44100)
{
    return (fs/ 1000) * ms;
}

static std::vector<double> test_string(size_t num_samples, double position)
{
    // theorical dual-polarization string
    dwg::string<double, 2> str (330.5409);

    // just generate a initial displacement centred at 28% left of the string size
    auto buffer = get_initial_displacement<double>(str.size(), position);

    // TODO: change the vertical polarization to be effective
    str.fill(buffer, 0);
    str.fill(buffer, 1);

    // main process
    std::vector<double> buffer;
    for (auto i = 0; i < num_samples; ++i)
    {
        auto hor = str.get(position, 0);
        auto ver = str.get(position, 1);
                   str.move();
        buffer.push_back(hor + ver);
    }
    return buffer;
}

int main()
{
    auto samples = test_string(samples_for_ms(5000), 0.28);

    return 0;
}

OK, this is a old code that take dust in my hard-drive since long date. I hope that will inspire you, at least for future plugins, at best for a rewrite and an inclusion in your framework. In any case, it's shared and it makes me happy.

Sincerely, Max

[jatinchowdhury18]

Thanks for sharing, yeah this looks super useful! I'd love to integrate the delay line and interpolation code into what I've already got in this repo, so it may take a little time to work out exactly how to integrate everything, but I do plan to do it.

Thanks, Jatin