Gaussian Distributed Random Numbers

Status
Not open for further replies.

gfvalvo

Well-known member
Hi all.
Wondering if any one knows of a Teensy-optimized function for generating random numbers from a Gaussian distribution. I'm looking at implementing the Box-Muller Transform myself. But, I'm not too proud to use an existing solution.
Thanks.
 
I have some code which implements the polar form of the Box_Muller transform. The original code was from here. It is not optimized except that I use it on a T3.6 which has an FPU and therefore the code is fast enough, for me anyway :).

Pete
 
Thanks, but that ftp site is asking for username / password. Do you have a link to your implementation, or can you post it?
 
Give it username anonymous and an email address as the password.

If that fails, I'll post my version.

Pete
 
If you're looking for this to run optimized, might be a good idea to change sqrt() to sqrtf() and log() to logf().

Normal sqrt() & log() without the extra "f" use 64 bit double, which is much slower than 32 bit float.
 
The Box-Muller algorithm is useful because it can be use to generate many kind of distributions aside from the Normal law. However, for Gaussian random variables, I think the best method is to use the one described in Numerical recipes in C++ which use a "ratio of uniform" method. It is a very elegant solution that generate perfect gaussian r.v. (up to the precision of the floating type used) and usually outperform Box-Muller...

Below is a simple example (using floats). On a Teensy 3.6, this outputs 2 000 000 r.v. per second. The RNG used here is not very good but it is very fast and will still pass many simple statistical test. If you want better stat. properties, you should consider using another RNG such as XorGen or a Mersenne Twister and make double precision computations but it will impact performance dramatically.



Code:
/**
* Very fast RNG. (for test purpose: no guarantee of good statistical properties). 
* The generator always use the same deterministic seed.
**/
class FastRNG
{
public:

    /* type of integer returned by the generator */
    typedef uint32_t result_type;

    /* Default constructor. Always initialize with the same seed. */
    FastRNG() : _gen_x(123456789), _gen_y(362436069), _gen_z(521288629) { }

    /* min value */
    static constexpr result_type min() { return 0; }

    /* max value */
    static constexpr result_type max() { return 4294967295UL; }

    /* return a random number */
    inline uint32_t operator()()
        {
        uint32_t t;
        _gen_x ^= _gen_x << 16; _gen_x ^= _gen_x >> 5; _gen_x ^= _gen_x << 1;
        t = _gen_x; _gen_x = _gen_y; _gen_y = _gen_z; _gen_z = t ^ _gen_x ^ _gen_y;
        return _gen_z;
        }

    /* return a uniform number on [0,1).  */
    inline float unif() { return ((float)(operator()())) / (4294967296.0); }

private:

    uint32_t _gen_x, _gen_y, _gen_z;        // state of the generator
};



/* generate N(0,1) random variable using "ratio of uniform" method 
   taken from "Numerical recipies in C++" 
 */
template<class random_t> inline float normalLaw(random_t & gen)
    {
    float u, v, x, y, q;
    do {
        u = gen.unif(); v = 1.7156*(gen.unif() - 0.5);
        x = u - 0.449871; y = abs(v) + 0.386595;
        q = (x*x) + y * (0.19600*y - 0.25472*x);
        } 
    while ((q > 0.27597) && (q > 0.27846 || (v*v) > -4.*logf(u)*(u*u)));
    return (v / u);
    }


/* test */
void setup()
    {
    delay(1000);
    while (!Serial) { yield(); }
    const int32_t N = 1000000;
    FastRNG gen;
    float mean = 0.0, var = 0.0; 
    elapsedMillis timer = 0;
    for (int32_t i = 0; i < N; i++)
        {
        const float x = normalLaw(gen); 
        mean += x;
        var += x*x; 
        }
    const uint32_t t = timer;
    Serial.print("mean : "); Serial.println(mean/N, 5);
    Serial.print("standart dev: "); Serial.println(sqrt(var/(N-1)), 5);
    Serial.print("generation took "); Serial.print(t); Serial.println("ms.");
    }


void loop()
{

}
 
Thanks for the replies.

Turns out, the straight forward code at the ftp site linked by @el_supremo has perfectly adequate performance for my application. No need for anything better in this case. As they say, “better is the enemy of good-enough”.
 
I needed a Gaussian generator for the Audio Vector Network Analyzer project. After getting that together, it seemed like others might be able to use this as an Teensy Audio class object. Thus this post. These are available for any use, and if someone sees bugs, I'd like to know :)

This uses the Knuth and Lewis quick uniform generator along with a Central Limit Theorem sum of 12 to transform to Gaussian. I have not tried all methods, but this seemed to be both fast and statistically valid for use in audio and audio testing. Also included is an optional 2-pole Butterworth (maximally flay) low pass filter to limit the bandwidth of the noise. This is implemented as a single stage of Biquad IIR. There are functions to set the standard deviation, filter cutoff frequency and the seed. See the .h file for detail.

Here are plots of a 1000 samples with and without filtering. Both data sets had the same seed, so the effect of filtering can be seen by comparing the two.
GaussianNoise1000.gif

The two files for object generation
View attachment synth_GaussianWhiteNoise.h
View attachment synth_GaussianWhiteNoise.cpp

and a test INO that I have run on T3.6 and T4.0
View attachment gaussianNoiseAndSignal.ino

I might note that the white noise generator in the Audio Library produces a uniform distribution. That gives a different sound than the Gaussian types. More important for some purposes is that many natural processes produce a Gaussian (sometimes called Normal or Bell-shaped) distribution.

Bob
 
Status
Not open for further replies.
Back
Top