Setting up floating-point RNGs

modules/dsp/chowdsp_compressor/Compressor/chowdsp_CompressorLevelDetector.h

@@ -45,7 +45,9 @@ class CompressorLevelDetector
                 newAttackMs = detectors.modifyAttack (newAttackMs);
         }
 
-        computeBallisticCoeffs (newAttackMs, fs);
+        const auto coeffs = computeBallisticCoeffs (newAttackMs, fs);
+        a1_a = coeffs.a1;
+        b0_a = coeffs.b0;
     }
 
     /** Sets the release time in milliseconds */
@@ -67,7 +69,9 @@ class CompressorLevelDetector
                 newReleaseMs = detectors.modifyRelease (newReleaseMs);
         }
 
-        computeBallisticCoeffs (newReleaseMs, fs);
+        const auto coeffs = computeBallisticCoeffs (newReleaseMs, fs);
+        a1_r = coeffs.a1;
+        b0_r = coeffs.b0;
     }
 
     /** Sets the threshold level in Decibels */

modules/dsp/chowdsp_math/Math/chowdsp_RandomFloat.h

@@ -0,0 +1,138 @@
+#pragma once
+
+#include <cstdint>
+
+namespace chowdsp
+{
+/** Helpers for generating random numbers. */
+namespace RandomUtils
+{
+    template <typename T>
+    struct BasisType;
+
+    template <>
+    struct BasisType<float>
+    {
+        using type = uint32_t;
+    };
+
+    template <>
+    struct BasisType<double>
+    {
+        using type = uint64_t;
+    };
+
+#if ! CHOWDSP_NO_XSIMD
+    template <>
+    struct BasisType<xsimd::batch<float>>
+    {
+        using type = xsimd::batch<uint32_t>;
+    };
+
+    template <>
+    struct BasisType<xsimd::batch<double>>
+    {
+        using type = xsimd::batch<uint64_t>;
+    };
+#endif
+
+    template <typename T>
+    using BasisTypeT = typename BasisType<T>::type;
+
+    constexpr float rng_0_1 (uint32_t& seed)
+    {
+        seed = seed * 196314165UL + 907633515UL;
+        return static_cast<float> (seed >> 8) / 16777216.0f;
+    }
+
+    constexpr double rng_0_1 (uint64_t& seed)
+    {
+        seed = seed * 6364136223846793005ULL + 1442695040888963407ULL;
+        return static_cast<double> (seed >> 11) / 9007199254740992.0;
+    }
+
+#if ! CHOWDSP_NO_XSIMD
+    inline xsimd::batch<float> rng_0_1 (xsimd::batch<uint32_t>& seed)
+    {
+        seed = seed * 196314165UL + 907633515UL;
+        return xsimd::batch_cast<float> (seed >> 8) / 16777216.0f;
+    }
+
+    inline xsimd::batch<double> rng_0_1 (xsimd::batch<uint64_t>& seed)
+    {
+        seed = seed * 6364136223846793005ULL + 1442695040888963407ULL;
+        return xsimd::batch_cast<double> (seed >> 11) / 9007199254740992.0;
+    }
+#endif
+
+    constexpr float rng_m1_1 (uint32_t& seed)
+    {
+        seed = seed * 196314165UL + 907633515UL;
+        auto temp = static_cast<int32_t> (seed >> 7) - 16777216L;
+        return static_cast<float> (temp) / 16777216.0f;
+    }
+
+    constexpr double rng_m1_1 (uint64_t& seed)
+    {
+        seed = seed * 6364136223846793005ULL + 1442695040888963407ULL;
+        auto temp = static_cast<int64_t> (seed >> 10) - 9007199254740992LL;
+        return static_cast<double> (temp) / 9007199254740992.0;
+    }
+
+#if ! CHOWDSP_NO_XSIMD
+    inline xsimd::batch<float> rng_m1_1 (xsimd::batch<uint32_t>& seed)
+    {
+        seed = seed * 196314165UL + 907633515UL;
+        auto temp = xsimd::batch_cast<int32_t> (seed >> 7) - 16777216L;
+        return xsimd::batch_cast<float> (temp) / 16777216.0f;
+    }
+
+    inline xsimd::batch<double> rng_m1_1 (xsimd::batch<uint64_t>& seed)
+    {
+        seed = seed * 6364136223846793005ULL + 1442695040888963407ULL;
+        auto temp = xsimd::batch_cast<int64_t> (seed >> 10) - 9007199254740992LL;
+        return xsimd::batch_cast<double> (temp) / 9007199254740992.0;
+    }
+#endif
+} // namespace RandomUtils
+
+/**
+ * Wrapper class for generating random floating-point numbers.
+ *
+ * If you need numbers generated on a symmetric two-sided range,
+ * (e.g. [-1, 1]), you'll get higher quality random numbers using
+ * the two-sided distribution (see template arguments).
+ *
+ * Reference: https://audiodev.blog/random-numbers/
+ */
+template <typename T, bool twoSided = false>
+class RandomFloat
+{
+public:
+    using Basis = RandomUtils::BasisTypeT<T>;
+
+    explicit RandomFloat (Basis s)
+        : seed (s)
+    {
+    }
+
+    T operator()()
+    {
+        if constexpr (twoSided)
+            return RandomUtils::rng_m1_1 (seed);
+        else
+            return RandomUtils::rng_0_1 (seed);
+    }
+
+    template <bool TS = twoSided>
+    std::enable_if_t<TS, T> operator() (T min, T max)
+    {
+        return min + (max - min) * RandomUtils::rng_0_1 (seed);
+    }
+
+private:
+    Basis seed {};
+
+    JUCE_DECLARE_NON_COPYABLE_WITH_LEAK_DETECTOR (RandomFloat)
+};
+} // namespace chowdsp

modules/dsp/chowdsp_math/chowdsp_math.h

@@ -37,6 +37,7 @@ JUCE_END_IGNORE_WARNINGS_GCC_LIKE
 #include "Math/chowdsp_MatrixOps.h"
 #include "Math/chowdsp_Polynomials.h"
 #include "Math/chowdsp_Power.h"
+#include "Math/chowdsp_RandomFloat.h"
 #include "Math/chowdsp_JacobiElliptic.h"
 #include "Math/chowdsp_Polylogarithm.h"
 #include "Math/chowdsp_TanhIntegrals.h"

tests/dsp_tests/chowdsp_math_test/CMakeLists.txt

@@ -15,4 +15,5 @@ target_sources(chowdsp_math_test
         PowApproxTest.cpp
         DecibelsApproxTest.cpp
         TrigApproxTest.cpp
+        RandomFloatTest.cpp
 )

tests/dsp_tests/chowdsp_math_test/RandomFloatTest.cpp

@@ -0,0 +1,100 @@
+#include <CatchUtils.h>
+#include <chowdsp_math/chowdsp_math.h>
+
+TEMPLATE_TEST_CASE ("Random Float Static Test", "[dsp][math][random]", float, double)
+{
+    using namespace chowdsp::RandomUtils;
+
+    std::array<int, 10> counters {};
+    static constexpr int N = 1'000'000;
+
+    using B = BasisTypeT<TestType>;
+    B seed = 444;
+
+    SECTION ("Range [0, 1]")
+    {
+        for (int n = 0; n < N; ++n)
+        {
+            const auto x = rng_0_1 (seed);
+            const auto bin = std::floor (x * (TestType) counters.size());
+            counters[(size_t) bin]++;
+        }
+
+        for (size_t i = 0; i < counters.size(); ++i)
+        {
+            const auto pct = (float) counters[i] * (float) counters.size() / (float) N;
+            REQUIRE (pct == Catch::Approx { 1.0f }.margin (0.1));
+        }
+    }
+
+    SECTION ("Range [-1, 1]")
+    {
+        for (int n = 0; n < N; ++n)
+        {
+            const auto x = rng_m1_1 (seed);
+            const auto bin = std::floor (std::abs (x) * (TestType) counters.size());
+            counters[(size_t) bin]++;
+        }
+
+        for (size_t i = 0; i < counters.size(); ++i)
+        {
+            const auto pct = (float) counters[i] * (float) counters.size() / (float) N;
+            REQUIRE (pct == Catch::Approx { 1.0f }.margin (0.1));
+        }
+    }
+}
+
+TEMPLATE_TEST_CASE ("Random SIMD Float Static Test", "[dsp][math][random]", float, double)
+{
+    using namespace chowdsp::RandomUtils;
+
+    std::array<int, 10> counters {};
+    static constexpr int N = 1'000'000;
+
+    using T = xsimd::batch<TestType>;
+    using B = BasisTypeT<T>;
+
+    B seed {};
+    if constexpr (T::size == 4)
+        seed = { 441, 442, 443, 444 };
+    else if constexpr (T::size == 2)
+        seed = { 440, 444 };
+
+    SECTION ("Range [0, 1]")
+    {
+        for (int n = 0; n < N; ++n)
+        {
+            const auto x = rng_0_1 (seed);
+            for (size_t k = 0; k < T::size; ++k)
+            {
+                const auto bin = std::floor (x.get (k) * (TestType) counters.size());
+                counters[(size_t) bin]++;
+            }
+        }
+
+        for (size_t i = 0; i < counters.size(); ++i)
+        {
+            const auto pct = (float) counters[i] * (float) counters.size() / (float) (N * T::size);
+            REQUIRE (pct == Catch::Approx { 1.0f }.margin (0.1));
+        }
+    }
+
+    SECTION ("Range [-1, 1]")
+    {
+        for (int n = 0; n < N; ++n)
+        {
+            const auto x = rng_m1_1 (seed);
+            for (size_t k = 0; k < T::size; ++k)
+            {
+                auto bin = std::floor (std::abs (x.get (k)) * (TestType) counters.size());
+                counters[(size_t) bin]++;
+            }
+        }
+
+        for (size_t i = 0; i < counters.size(); ++i)
+        {
+            const auto pct = (float) counters[i] * (float) counters.size() / (float) (N * T::size);
+            REQUIRE (pct == Catch::Approx { 1.0f }.margin (0.1));
+        }
+    }
+}