diff options
| author | Stefan <stefan.wunsch@student.kit.edu> | 2015-09-01 15:09:17 +0200 |
|---|---|---|
| committer | Stefan <stefan.wunsch@student.kit.edu> | 2015-09-01 15:09:17 +0200 |
| commit | e172d55e2bfe63e5b99e7d432f9386c37bff8c84 (patch) | |
| tree | 1a45aea0aabb5c948c3c6a0b0f9969999871431f /gnuradio-runtime/lib/math/random.cc | |
| parent | ebad2162b40b2b144449ff2927dbe89c683c4972 (diff) | |
add boost.random as random number generator
Diffstat (limited to 'gnuradio-runtime/lib/math/random.cc')
| -rw-r--r-- | gnuradio-runtime/lib/math/random.cc | 106 |
1 files changed, 40 insertions, 66 deletions
diff --git a/gnuradio-runtime/lib/math/random.cc b/gnuradio-runtime/lib/math/random.cc index ecf70e0f9a..76de3bdde5 100644 --- a/gnuradio-runtime/lib/math/random.cc +++ b/gnuradio-runtime/lib/math/random.cc @@ -44,90 +44,67 @@ namespace gr { -#define IA 16807 -#define IM 2147483647 -#define AM (1.0/IM) -#define IQ 127773 -#define IR 2836 -#define NDIV (1+(IM-1)/NTAB) -#define EPS 1.2e-7 -#define RNMX (1.0-EPS) + random::random(unsigned int seed) + { + d_gauss_stored = false; // set gasdev (gauss distributed numbers) on calculation state + + // Setup random number generator + d_rng = new boost::mt19937; + reseed(seed); // set seed for random number generator + d_uniform = new boost::uniform_real<float>; + d_generator = new boost::variate_generator<boost::mt19937&, boost::uniform_real<float> > (*d_rng,*d_uniform); // create number generator in [0,1) from boost.random + } - random::random(long seed) + random::~random() { - reseed(seed); + delete d_rng; + delete d_uniform; + delete d_generator; } + /* + * Seed is initialized with time if the given seed is 0. Otherwise the seed is taken directly. Sets the seed for the random number generator. + */ void - random::reseed(long seed) + random::reseed(unsigned int seed) { - d_seed = seed; - d_iy = 0; - for(int i = 0; i < NTAB; i++) - d_iv[i] = 0; - d_iset = 0; - d_gset = 0; + if(seed==0) d_seed = static_cast<unsigned int>(std::time(0)); // FIXME: add seed method correctly + else d_seed = seed; + d_rng->seed(d_seed); } /* - * This looks like it returns a uniform random deviate between 0.0 and 1.0 - * It looks similar to code from "Numerical Recipes in C". + * Returns uniformly distributed numbers in [0,1) taken from boost.random using a Mersenne twister */ float random::ran1() { - int j; - long k; - float temp; - - if(d_seed <= 0 || !d_iy) { - if(-d_seed < 1) - d_seed=1; - else - d_seed = -d_seed; - for(j=NTAB+7;j>=0;j--) { - k=d_seed/IQ; - d_seed=IA*(d_seed-k*IQ)-IR*k; - if(d_seed < 0) - d_seed += IM; - if(j < NTAB) - d_iv[j] = d_seed; - } - d_iy=d_iv[0]; - } - k=(d_seed)/IQ; - d_seed=IA*(d_seed-k*IQ)-IR*k; - if(d_seed < 0) - d_seed += IM; - j=d_iy/NDIV; - d_iy=d_iv[j]; - d_iv[j] = d_seed; - temp=AM * d_iy; - if(temp > RNMX) - temp = RNMX; - return temp; + return (*d_generator)(); } /* * Returns a normally distributed deviate with zero mean and variance 1. - * Also looks like it's from "Numerical Recipes in C". + * Used is the Marsaglia polar method. + * Every second call a number is stored because the transformation works only in pairs. Otherwise half calculation is thrown away. */ float random::gasdev() { - float fac,rsq,v1,v2; - d_iset = 1 - d_iset; - if(d_iset) { - do { - v1=2.0*ran1()-1.0; - v2=2.0*ran1()-1.0; - rsq=v1*v1+v2*v2; - } while(rsq >= 1.0 || rsq == 0.0); - fac= sqrt(-2.0*log(rsq)/rsq); - d_gset=v1*fac; - return v2*fac; - } - return d_gset; + if(d_gauss_stored){ // just return the stored value if available + d_gauss_stored = false; + return d_gauss_value; + } + else{ // generate a pair of gaussian distributed numbers + float x,y,s; + do{ + x = 2.0*ran1()-1.0; + y = 2.0*ran1()-1.0; + s = x*x+y*y; + }while(not(s<1.0)); + d_gauss_stored = true; + d_gauss_value = x*sqrt(-2.0*log(s)/s); + return y*sqrt(-2.0*log(s)/s); + } } float @@ -153,9 +130,6 @@ namespace gr { return z; } - /* - * Complex rayleigh is really gaussian I and gaussian Q - */ gr_complex random::rayleigh_complex() { |