Oversamping in C++

[Allen]

Has anybody done oversampling in C++ nodes?

Since the oversampling nodes are not for poly nodes so I'm trying to directly oversample inside c++ now

[Dan Korneff]

@Allen I'm accessing it directly in my C++ nodes with

juce::dsp::Oversampling<float>;

[griffinboy]

@Allen

I made this custom one at one point when I needed stronger alias rejection.
It's old though and I only used it once. just posting an alternative in case it's useful.

Something you'll need to be careful of when oversampling is that the antialising filters usually introduce latency.
So if you have things that are oversampled mixed with things that aren't, make sure to phase align them by compensating latency.

// FILE: src/Oversampler2xKaiser.h
/*
    2x Kaiser-windowed FIR oversampler/downsamlper 
    Linear-phase, unity DC gain. Default taps=63, beta=8.0.

    USAGE:
      - Oversampler2xKaiser os; os.prepare(63, 8.0); os.reset();
      - double u0,u1; os.upsample(x, u0, u1);                   // generate 2x stream
      - double y = os.downsample(v0, v1);                        // consume pair from 2x domain
*/

#pragma once
#include <vector>
#include <cmath>
#include <cstring>
#include "Utils_General.h"


#define OS2X_HAVE_XSIMD 1


class Oversampler2xKaiser
{
public:
    Oversampler2xKaiser() noexcept = default;

    void prepare(int taps = 63, double beta = 8.0)
    {
        if ((taps & 1) == 0) ++taps;
        N = taps;
        B = beta;
        buildCoeffs();

        upHist.assign((size_t)H, 0.0);
        dEven.assign((size_t)L0, 0.0);
        dOdd.assign((size_t)L1, 0.0);
    }

    void reset() noexcept
    {
        std::fill(upHist.begin(), upHist.end(), 0.0);
        std::fill(dEven.begin(), dEven.end(), 0.0);
        std::fill(dOdd.begin(), dOdd.end(), 0.0);
    }

    inline void upsample(double x, double& y0, double& y1) noexcept
    {
        shift_in(upHist.data(), H, x);
        y0 = dotSIMD(p0.data(), upHist.data(), L0);
        y1 = dotSIMD(p1.data(), upHist.data(), L1);
    }

    inline double downsample(double v0, double v1) noexcept
    {
        shift_in(dEven.data(), L0, v0);
        shift_in(dOdd.data(), L1, v1);
        const double a = dotSIMD(p0.data(), dEven.data(), L0);
        const double b = dotSIMD(p1.data(), dOdd.data(), L1);
        return a + b;
    }

private:
    // PURPOSE: Kaiser window, halfband ideal (wc=pi/2) + normalization. Build polyphase strides p0/p1.
    void buildCoeffs()
    {
        const int M = (N - 1) / 2;
        std::vector<double> h((size_t)N);
        const double kPi = gdsp::PI;
        const double i0B = i0(B);

        for (int n = 0; n < N; ++n)
        {
            const int m = n - M;
            const double r = (double)m / (double)M;
            const double w = i0(B * std::sqrt(std::max(0.0, 1.0 - r * r))) / i0B;

            double sinc;
            if (m == 0) sinc = 0.5;
            else         sinc = std::sin(0.5 * kPi * m) / (kPi * m);

            h[(size_t)n] = w * sinc;
        }
        double s = 0.0; for (double v : h) s += v;
        const double invS = (s != 0.0) ? (1.0 / s) : 1.0;
        for (double& v : h) v *= invS;

        L0 = (N + 1) / 2;
        L1 = N / 2;
        H = std::max(L0, L1);

        p0.resize((size_t)L0);
        p1.resize((size_t)L1);

        for (int k = 0; k < L0; ++k) p0[(size_t)k] = h[(size_t)(2 * k)];
        for (int k = 0; k < L1; ++k) p1[(size_t)k] = h[(size_t)(2 * k + 1)];
    }

    static inline void shift_in(double* hist, int len, double x) noexcept
    {
        std::memmove(hist + 1, hist, sizeof(double) * (size_t)(len - 1));
        hist[0] = x;
    }

    // PURPOSE: SIMD dot with scalar fallback.
    static inline double dotSIMD(const double* a, const double* b, int n) noexcept
    {
#if OS2X_HAVE_XSIMD
        using vd = xsimd::batch<double>;
        constexpr int W = (int)vd::size;
        int i = 0;
        vd acc = vd::broadcast(0.0);
        for (; i + W <= n; i += W)
            acc = xsimd::fma(xsimd::load_aligned(a + i), xsimd::load_aligned(b + i), acc);
        double s = xsimd::reduce_add(acc);
        for (; i < n; ++i) s += a[i] * b[i];
        return s;
#else
        double s = 0.0;
        for (int i = 0; i < n; ++i) s += a[i] * b[i];
        return s;
#endif
    }

    // PURPOSE: Modified Bessel I0 via series. Accurate for |x| up to ~12 within ~1e-12.
    static inline double i0(double x) noexcept
    {
        double t = 1.0;
        double y = (x * x) * 0.25;
        double s = 1.0;
        for (int k = 1; k < 50; ++k)
        {
            t *= y / (double)(k * k);
            s += t;
            if (t < 1e-12 * s) break;
        }
        return s;
    }

    int N{ 63 }, L0{ 0 }, L1{ 0 }, H{ 0 };
    double B{ 8.0 };

#if OS2X_HAVE_XSIMD
    template<typename T> using AAlloc = xsimd::aligned_allocator<T, 64>;
#else
    template<typename T> using AAlloc = std::allocator<T>;
#endif

    std::vector<double, AAlloc<double>> p0, p1;
    std::vector<double, AAlloc<double>> upHist, dEven, dOdd;
};

/*
    Utils_General.h

	General useful math for dsp programming
*/


#pragma once

#include <climits>
#include <cmath>

#if defined(_MSC_VER)
#pragma once
#pragma warning(4 : 4250)          // "Inherits via dominance" 
#endif

namespace gdsp
{

	//===================| Constants |===================

	// PI
	static const double  PI = 3.1415926535897932384626433832795;

    static const double TWOPI = 2 * PI;

	// FORCEINLINE
    #if defined(_MSC_VER)
    #define FORCEINLINE __forceinline
    #else
    #define FORCEINLINE inline
    #endif



    //===================| Fixed Point Math |===================

     /*----- 16-bit signed -----*/
#if defined(_MSC_VER)
    typedef __int16                 Int16;
#elif (defined(__MWERKS__) || defined(__GNUC__) || defined(__BEOS__)) && SHRT_MAX == 0x7FFF
    typedef short int               Int16;
#else
#error No signed 16-bit integer type defined for this compiler!
#endif

    /*----- 32-bit signed -----*/
#if defined(_MSC_VER)
    typedef __int32                 Int32;
#elif (defined(__MWERKS__) || defined(__GNUC__) || defined(__BEOS__)) && INT_MAX == 0x7FFFFFFFL
    typedef int                     Int32;
#else
#error No signed 32-bit integer type defined for this compiler!
#endif

    /*----- 64-bit signed -----*/
#if defined(_MSC_VER)
    typedef __int64                 Int64;
#elif defined(__MWERKS__) || defined(__GNUC__)
    typedef long long               Int64;
#elif defined(__BEOS__)
    typedef int64                   Int64;
#else
#error No 64-bit integer type defined for this compiler!
#endif

    /*----- 32-bit unsigned ----*/
#if defined(_MSC_VER)
    typedef unsigned __int32        UInt32;
#elif (defined(__MWERKS__) || defined(__GNUC__) || defined(__BEOS__)) && UINT_MAX == 0xFFFFFFFFUL
    typedef unsigned int            UInt32;
#else
#error No unsigned 32-bit integer type defined for this compiler!
#endif

	// Helpers for fixed-point math
    union Fixed3232
    {
        Int64  _all;
        class
        {
        public:
#if defined(__POWERPC__)            /* big-endian */
            Int32   _msw;
            UInt32  _lsw;
#else                               /* little-endian */
            UInt32  _lsw;
            Int32   _msw;
#endif
        } _part;
    }; 

     // simple branch-free min/max 
    template <typename T> FORCEINLINE T min(T a, T b) { return (a < b) ? a : b; }
    template <typename T> FORCEINLINE T max(T a, T b) { return (a > b) ? a : b; }

    // IEEE-safe rounding 
    FORCEINLINE int  round_int(double x) { return (x >= 0.0) ? int(x + 0.5) : int(x - 0.5); }
    FORCEINLINE long round_long(double x) { return (x >= 0.0) ? long(x + 0.5) : long(x - 0.5); }

    // bidirectional (signed) bit-shift helper
    template <typename T>
    FORCEINLINE T shift_bidi(T x, int s)
    {
        return (s >= 0) ? (x << s) : (x >> (-s));
    }



} // namespace gdsp

[griffinboy]

@weston-donald

Yes agreed.
Make sure to match your upsampling / downsampling filters.
Your upsample stage and downsample stage need to be designed together so that they don't cause bad ripple or other artefacts. Mismatched filters will cause aliasing even (going from 2x as many samples, down to half as many may create discontinuities if not smoothed / interpolated proper).

[Christoph Hart]

@griffinboy I think you're agreeing with an AI bot here.

I'm watching you closely, Weston, another generic post without any meaningful contribution and you're out of here!

[griffinboy]

@Christoph-Hart

Shoot you're right.
I didn't even think of that.
I wonder what the point of Ai botting a forum like this would be

[Christoph Hart]

@griffinboy probably an elaborated version of the bots that I'm throwing out here on a regular basis.

I would say four more posts from Weston about stuff with tangential overlap on the subject, then BAM! Link to an Indian Escort Service...

[griffinboy]

@Christoph-Hart ah, right. Yes looking at the profile it's very obvious.

[Allen]

@griffinboy
Thank you!! I'll take a look into it!