From 8c9b53915efe482f12f4d45828ad048db6ebda71 Mon Sep 17 00:00:00 2001 From: Jean Pierre Cimalando Date: Mon, 22 Jun 2020 02:00:37 +0200 Subject: [PATCH] Fast gauss distribution --- benchmarks/BM_random.cpp | 23 ++++++++++ src/sfizz/MathHelpers.h | 44 ++++++++++++++++++ src/sfizz/Voice.cpp | 4 +- src/sfizz/Voice.h | 2 +- tests/RandomHelpersT.cpp | 98 +++++++++++++++++++++++++++++++++++++++- 5 files changed, 166 insertions(+), 5 deletions(-) diff --git a/benchmarks/BM_random.cpp b/benchmarks/BM_random.cpp index 85890887..d5220322 100644 --- a/benchmarks/BM_random.cpp +++ b/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 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 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); diff --git a/src/sfizz/MathHelpers.h b/src/sfizz/MathHelpers.h index 70ce204a..b446d873 100644 --- a/src/sfizz/MathHelpers.h +++ b/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 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(next) * (1.0f / (1ll << 31)); + } + return mean_ + gain_ * sum; + } + +private: + std::array seeds_ {}; + float mean_ { 0 }; + float gain_ { 0 }; +}; diff --git a/src/sfizz/Voice.cpp b/src/sfizz/Voice.cpp index 7e384345..73c71a42 100644 --- a/src/sfizz/Voice.cpp +++ b/src/sfizz/Voice.cpp @@ -600,8 +600,8 @@ void sfz::Voice::fillWithGenerator(AudioSpan 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(); diff --git a/src/sfizz/Voice.h b/src/sfizz/Voice.h index d2b25ad8..54b03318 100644 --- a/src/sfizz/Voice.h +++ b/src/sfizz/Voice.h @@ -467,7 +467,7 @@ private: Voice* nextSisterVoice { this }; Voice* previousSisterVoice { this }; - std::normal_distribution noiseDist { 0, config::noiseVariance }; + fast_gaussian_generator noiseDist { 0.0f, config::noiseVariance }; ModifierArray> modifierSmoothers; Smoother gainSmoother; diff --git a/tests/RandomHelpersT.cpp b/tests/RandomHelpersT.cpp index 9c76fa80..89df6c05 100644 --- a/tests/RandomHelpersT.cpp +++ b/tests/RandomHelpersT.cpp @@ -8,6 +8,12 @@ #include "sfizz/MathHelpers.h" #include +template +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 -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 dist(min, max); @@ -41,10 +47,98 @@ TEST_CASE("[Random] Fast random generation") unsigned numGenerations = 1000; unsigned numDivisions = 128; - auto& test = randomTest>; + auto& test = uniformRandomTest>; 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 +static bool gaussianRandomTest(double mean, float variance, size_t numGen, size_t histSize, double maxAbsErr) +{ + fast_gaussian_generator gen(mean, variance); + + std::mt19937_64 prng; + std::normal_distribution dist(mean, variance); + + // the tested bounds + const double min = -1.0 + mean; + const double max = +1.0 + mean; + + // generate, quantify, count occurrences + std::vector 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(bin) < histSize) + counts[bin] += 1; + } + std::vector 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(bin) < histSize) + referenceCounts[bin] += 1; + } + + // bin probabilities + std::vector proba(histSize); + std::vector referenceProba(histSize); + for (size_t i = 0; i < histSize; ++i) { + proba[i] = static_cast(counts[i]) / numGen; + referenceProba[i] = static_cast(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 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)); +}