HISE Logo Forum
    • Categories
    • Register
    • Login

    Oversamping in C++

    Scheduled Pinned Locked Moved C++ Development
    9 Posts 4 Posters 598 Views
    Loading More Posts
    • Oldest to Newest
    • Newest to Oldest
    • Most Votes
    Reply
    • Reply as topic
    Log in to reply
    This topic has been deleted. Only users with topic management privileges can see it.
    • A
      Allen
      last edited by

      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 KorneffD griffinboyG 2 Replies Last reply Reply Quote 0
      • Dan KorneffD
        Dan Korneff @Allen
        last edited by

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

        juce::dsp::Oversampling<float>;
        

        Dan Korneff - Producer / Mixer / Audio Nerd

        1 Reply Last reply Reply Quote 1
        • griffinboyG
          griffinboy @Allen
          last edited by 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
          
          A 1 Reply Last reply Reply Quote 2
          • griffinboyG
            griffinboy
            last edited by 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 HartC 1 Reply Last reply Reply Quote 0
            • Christoph HartC
              Christoph Hart @griffinboy
              last edited by

              @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!

              griffinboyG 1 Reply Last reply Reply Quote 3
              • griffinboyG
                griffinboy @Christoph Hart
                last edited by

                @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 HartC 1 Reply Last reply Reply Quote 0
                • Christoph HartC
                  Christoph Hart @griffinboy
                  last edited by

                  @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...

                  griffinboyG 1 Reply Last reply Reply Quote 3
                  • griffinboyG
                    griffinboy @Christoph Hart
                    last edited by griffinboy

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

                    1 Reply Last reply Reply Quote 0
                    • A
                      Allen @griffinboy
                      last edited by

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

                      1 Reply Last reply Reply Quote 0
                      • First post
                        Last post

                      21

                      Online

                      2.0k

                      Users

                      12.6k

                      Topics

                      109.5k

                      Posts