Merge pull request #286 from jpcima/gauss

Fast gauss distribution

benchmarks/BM_random.cpp

@@ -65,7 +65,30 @@ BENCHMARK_DEFINE_F(RandomFill, FastRandomBipolar)(benchmark::State& state) {
     }
 }
 
+BENCHMARK_DEFINE_F(RandomFill, StdNormal)(benchmark::State& state) {
+    std::minstd_rand prng;
+    std::normal_distribution<float> dist(0.0f, 0.25f);
+
+    for (auto _ : state)
+    {
+        for (float &out : output)
+            out = dist(prng);
+    }
+}
+
+BENCHMARK_DEFINE_F(RandomFill, FastNormal)(benchmark::State& state) {
+    fast_gaussian_generator<float, 4> generator(0.0f, 0.25f);
+
+    for (auto _ : state)
+    {
+        for (float &out : output)
+            out = generator();
+    }
+}
+
 BENCHMARK_REGISTER_F(RandomFill, StdRandom)->RangeMultiplier(4)->Range(1 << 2, 1 << 12);
 BENCHMARK_REGISTER_F(RandomFill, StdRandomBipolar)->RangeMultiplier(4)->Range(1 << 2, 1 << 12);
 BENCHMARK_REGISTER_F(RandomFill, FastRandom)->RangeMultiplier(4)->Range(1 << 2, 1 << 12);
 BENCHMARK_REGISTER_F(RandomFill, FastRandomBipolar)->RangeMultiplier(4)->Range(1 << 2, 1 << 12);
+BENCHMARK_REGISTER_F(RandomFill, StdNormal)->RangeMultiplier(4)->Range(1 << 2, 1 << 12);
+BENCHMARK_REGISTER_F(RandomFill, FastNormal)->RangeMultiplier(4)->Range(1 << 2, 1 << 12);

src/sfizz/MathHelpers.h

@@ -619,3 +619,47 @@ private:
 namespace Random {
 static fast_rand randomGenerator;
 } // namespace Random
+
+/**
+ * @brief Generate normally distributed noise.
+ *
+ * This sums the output of N uniform random generators.
+ * The higher the N, the better is the approximation of a normal distribution.
+ */
+template <class T, unsigned N = 4>
+class fast_gaussian_generator {
+    static_assert(N > 1, "Invalid quality setting");
+
+public:
+    explicit fast_gaussian_generator(float mean, float variance, uint32_t initialSeed = Random::randomGenerator())
+    {
+        mean_ = mean;
+        gain_ = variance / std::sqrt(N / 3.0);
+        seed(initialSeed);
+    }
+
+    void seed(uint32_t s)
+    {
+        seeds_[0] = s;
+        for (unsigned i = 1; i < N; ++i) {
+            s += s * 1664525u + 1013904223u;
+            seeds_[i] = s;
+        }
+    }
+
+    float operator()() noexcept
+    {
+        float sum = 0;
+        for (unsigned i = 0; i < N; ++i) {
+            uint32_t next = seeds_[i] * 1664525u + 1013904223u;
+            seeds_[i] = next;
+            sum += static_cast<int32_t>(next) * (1.0f / (1ll << 31));
+        }
+        return mean_ + gain_ * sum;
+    }
+
+private:
+    std::array<uint32_t, N> seeds_ {};
+    float mean_ { 0 };
+    float gain_ { 0 };
+};

src/sfizz/Voice.cpp

@@ -600,8 +600,8 @@ void sfz::Voice::fillWithGenerator(AudioSpan<float> buffer) noexcept
     const auto rightSpan  = buffer.getSpan(1);
 
     if (region->sampleId.filename() == "*noise") {
-        absl::c_generate(leftSpan, [&](){ return noiseDist(Random::randomGenerator); });
-        absl::c_generate(rightSpan, [&](){ return noiseDist(Random::randomGenerator); });
+        absl::c_generate(leftSpan, noiseDist);
+        absl::c_generate(rightSpan, noiseDist);
     } else {
         const auto numFrames = buffer.getNumFrames();
 

src/sfizz/Voice.h

@@ -467,7 +467,7 @@ private:
     Voice* nextSisterVoice { this };
     Voice* previousSisterVoice { this };
 
-    std::normal_distribution<float> noiseDist { 0, config::noiseVariance };
+    fast_gaussian_generator<float> noiseDist { 0.0f, config::noiseVariance };
 
     ModifierArray<std::vector<Smoother>> modifierSmoothers;
     Smoother gainSmoother;

tests/RandomHelpersT.cpp

@@ -8,6 +8,12 @@
 #include "sfizz/MathHelpers.h"
 #include <vector>
 
+template <class T>
+static inline T squared(T x)
+{
+    return x * x;
+}
+
 /**
  * @brief Check the behavior of a uniform real random generator
  *
@@ -15,7 +21,7 @@
  * - ensure to have at least one element in every division of the result range
  */
 template <class Rand, class Dist>
-static bool randomTest(float min, float max, unsigned gen, unsigned div)
+static bool uniformRandomTest(float min, float max, unsigned gen, unsigned div)
 {
     fast_rand prng;
     fast_real_distribution<float> dist(min, max);
@@ -41,10 +47,98 @@ TEST_CASE("[Random] Fast random generation")
     unsigned numGenerations = 1000;
     unsigned numDivisions = 128;
 
-    auto& test = randomTest<fast_rand, fast_real_distribution<float>>;
+    auto& test = uniformRandomTest<fast_rand, fast_real_distribution<float>>;
 
     REQUIRE(test(0.0f, 1.0f, numGenerations, numDivisions));
     REQUIRE(test(-1.0f, 1.0f, numGenerations, numDivisions));
     REQUIRE(test(0.0f, 123.0f, numGenerations, numDivisions));
     REQUIRE(test(-123.0f, 0.0f, numGenerations, numDivisions));
 }
+
+/**
+ * @brief Check the behavior of a gaussian real random generator
+ *
+ * - compute a histogram of random generations
+ * - compare it against the standard library normal distribution
+ */
+template <unsigned Quality>
+static bool gaussianRandomTest(double mean, float variance, size_t numGen, size_t histSize, double maxAbsErr)
+{
+    fast_gaussian_generator<float, Quality> gen(mean, variance);
+
+    std::mt19937_64 prng;
+    std::normal_distribution<float> dist(mean, variance);
+
+    // the tested bounds
+    const double min = -1.0 + mean;
+    const double max = +1.0 + mean;
+
+    // generate, quantify, count occurrences
+    std::vector<size_t> counts(histSize);
+    for (size_t i = 0; i < numGen; ++i) {
+        double value = gen();
+        double normalized = (value - min) / (max - min);
+        long bin = std::lround(histSize * normalized);
+        if (bin >= 0 && static_cast<size_t>(bin) < histSize)
+            counts[bin] += 1;
+    }
+    std::vector<size_t> referenceCounts(histSize);
+    for (size_t i = 0; i < numGen; ++i) {
+        double value = dist(prng);
+        double normalized = (value - min) / (max - min);
+        long bin = std::lround(histSize * normalized);
+        if (bin >= 0 && static_cast<size_t>(bin) < histSize)
+            referenceCounts[bin] += 1;
+    }
+
+    // bin probabilities
+    std::vector<double> proba(histSize);
+    std::vector<double> referenceProba(histSize);
+    for (size_t i = 0; i < histSize; ++i) {
+        proba[i] = static_cast<double>(counts[i]) / numGen;
+        referenceProba[i] = static_cast<double>(referenceCounts[i]) / numGen;
+    }
+
+    // normalize to unity, for the sake of comparison
+    {
+        double sum = 0.0;
+        for (double x : referenceProba)
+            sum += x;
+        double scaleFactor = 1.0 / sum;
+        for (size_t i = 0; i < histSize; ++i) {
+            proba[i] *= scaleFactor;
+            referenceProba[i] *= scaleFactor;
+        }
+    }
+
+    auto normalizeToUnity = [](absl::Span<double> v) {
+        double s = 0;
+        for (double x : v)
+            s += x;
+        for (double& x : v)
+            x *= 1.0 / s;
+    };
+    normalizeToUnity(absl::MakeSpan(proba));
+    normalizeToUnity(absl::MakeSpan(referenceProba));
+
+    // compare
+    for (size_t i = 0; i < histSize; ++i) {
+        double absErr = std::abs(proba[i] - referenceProba[i]);
+        if (absErr > maxAbsErr)
+            return false;
+    }
+
+    return true;
+}
+
+TEST_CASE("[Random] Gaussian random generation")
+{
+    unsigned numGenerations = 16384;
+    unsigned numDivisions = 128;
+
+    const double maxAbsErr = 0.05; // PDF ±5%
+
+    REQUIRE(gaussianRandomTest<4>(0.0f, 0.25f, numGenerations, numDivisions, maxAbsErr));
+    REQUIRE(gaussianRandomTest<4>(0.0f, 0.50f, numGenerations, numDivisions, maxAbsErr));
+    REQUIRE(gaussianRandomTest<4>(0.0f, 0.75f, numGenerations, numDivisions, maxAbsErr));
+}