Merge branch 'master' into next

Conflicts: gr-atsc/lib/qa_atsci_fake_single_viterbi.cc gr-atsc/lib/qa_atsci_fs_correlator.cc gr-atsc/lib/qa_atsci_single_viterbi.cc
author: Johnathan Corgan <johnathan@corganlabs.com> 2015-09-05 18:05:13 -0700
committer: Johnathan Corgan <johnathan@corganlabs.com> 2015-09-05 18:05:13 -0700
commit: 24c336cd99d8012681d40ccb00b1966e3c28450f (patch)
tree: 45cf228c11a1669f0e62bee13e3efbad17fcd4bd /gnuradio-runtime
parent: 273e5068d6d418eebb77dc0b5e15da20018bb90f (diff)
parent: 45c0feeaeffde561758a4d92bf4463d60b85633a (diff)
5 files changed, 282 insertions, 108 deletions
diff --git a/gnuradio-runtime/apps/evaluation_random_numbers.py b/gnuradio-runtime/apps/evaluation_random_numbers.py
new file mode 100644
index 0000000000..069493c73e
--- /dev/null
+++ b/gnuradio-runtime/apps/evaluation_random_numbers.py
@@ -0,0 +1,139 @@
+#!/usr/bin/env python
+#
+# Copyright 2015 Free Software Foundation, Inc.
+#
+# This file is part of GNU Radio
+#
+# GNU Radio is free software; you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation; either version 3, or (at your option)
+# any later version.
+#
+# GNU Radio is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with GNU Radio; see the file COPYING.  If not, write to
+# the Free Software Foundation, Inc., 51 Franklin Street,
+# Boston, MA 02110-1301, USA.
+#
+
+from gnuradio import gr
+import numpy as np
+from scipy.stats import norm, laplace, rayleigh
+from matplotlib import pyplot as plt
+
+# NOTE: scipy and matplotlib are optional packages and not included in the default gnuradio dependencies
+
+#*** SETUP ***#
+
+# Number of realisations per histogram
+num_tests = 1000000
+
+# Set number of bins in histograms
+uniform_num_bins = 31
+gauss_num_bins = 31
+rayleigh_num_bins = 31
+laplace_num_bins = 31
+
+rndm = gr.random() # instance of gnuradio random class (gr::random)
+
+print 'All histograms contain',num_tests,'realisations.'
+
+#*** GENERATE DATA ***#
+
+uniform_values = np.zeros(num_tests)
+gauss_values = np.zeros(num_tests)
+rayleigh_values = np.zeros(num_tests)
+laplace_values = np.zeros(num_tests)
+
+for k in range(num_tests):
+    uniform_values[k] = rndm.ran1()
+    gauss_values[k] = rndm.gasdev()
+    rayleigh_values[k] = rndm.rayleigh()
+    laplace_values[k] = rndm.laplacian()
+
+#*** HISTOGRAM DATA AND CALCULATE EXPECTED COUNTS ***#
+
+uniform_bins = np.linspace(0,1,uniform_num_bins)
+gauss_bins = np.linspace(-8,8,gauss_num_bins)
+laplace_bins = np.linspace(-8,8,laplace_num_bins)
+rayleigh_bins = np.linspace(0,10,rayleigh_num_bins)
+
+uniform_hist = np.histogram(uniform_values,uniform_bins)
+gauss_hist = np.histogram(gauss_values,gauss_bins)
+rayleigh_hist = np.histogram(rayleigh_values,rayleigh_bins)
+laplace_hist = np.histogram(laplace_values,laplace_bins)
+
+uniform_expected = np.zeros(uniform_num_bins-1)
+gauss_expected = np.zeros(gauss_num_bins-1)
+rayleigh_expected = np.zeros(rayleigh_num_bins-1)
+laplace_expected = np.zeros(laplace_num_bins-1)
+
+for k in range(len(uniform_hist[0])):
+    uniform_expected[k] = num_tests/float(uniform_num_bins-1)
+
+for k in range(len(gauss_hist[0])):
+    gauss_expected[k] = float(norm.cdf(gauss_hist[1][k+1])-norm.cdf(gauss_hist[1][k]))*num_tests
+
+for k in range(len(rayleigh_hist[0])):
+    rayleigh_expected[k] = float(rayleigh.cdf(rayleigh_hist[1][k+1])-rayleigh.cdf(rayleigh_hist[1][k]))*num_tests
+
+for k in range(len(laplace_hist[0])):
+    laplace_expected[k] = float(laplace.cdf(laplace_hist[1][k+1])-laplace.cdf(laplace_hist[1][k]))*num_tests
+
+#*** PLOT HISTOGRAMS AND EXPECTATIONS TAKEN FROM SCIPY ***#
+
+uniform_bins_center = uniform_bins[0:-1]+(uniform_bins[1]-uniform_bins[0])/2.0
+gauss_bins_center = gauss_bins[0:-1]+(gauss_bins[1]-gauss_bins[0])/2.0
+rayleigh_bins_center = rayleigh_bins[0:-1]+(rayleigh_bins[1]-rayleigh_bins[0])/2.0
+laplace_bins_center = laplace_bins[0:-1]+(laplace_bins[1]-laplace_bins[0])/2.0
+
+plt.figure(1)
+
+plt.subplot(2,1,1)
+plt.plot(uniform_bins_center,uniform_hist[0],'s--',uniform_bins_center,uniform_expected,'o:')
+plt.xlabel('Bins'), plt.ylabel('Count'), plt.title('Uniform: Distribution')
+plt.legend(['histogram gr::random','calculation scipy'],loc=1)
+
+plt.subplot(2,1,2)
+plt.plot(uniform_bins_center,uniform_hist[0]/uniform_expected,'rs--')
+plt.xlabel('Bins'), plt.ylabel('Relative deviation'), plt.title('Uniform: Relative deviation to scipy')
+
+plt.figure(2)
+
+plt.subplot(2,1,1)
+plt.plot(gauss_bins_center,gauss_hist[0],'s--',gauss_bins_center,gauss_expected,'o:')
+plt.xlabel('Bins'), plt.ylabel('Count'), plt.title('Gauss: Distribution')
+plt.legend(['histogram gr::random','calculation scipy'],loc=1)
+
+plt.subplot(2,1,2)
+plt.plot(gauss_bins_center,gauss_hist[0]/gauss_expected,'rs--')
+plt.xlabel('Bins'), plt.ylabel('Relative deviation'), plt.title('Gauss: Relative deviation to scipy')
+
+plt.figure(3)
+
+plt.subplot(2,1,1)
+plt.plot(rayleigh_bins_center,rayleigh_hist[0],'s--',rayleigh_bins_center,rayleigh_expected,'o:')
+plt.xlabel('Bins'), plt.ylabel('Count'), plt.title('Rayleigh: Distribution')
+plt.legend(['histogram gr::random','calculation scipy'],loc=1)
+
+
+plt.subplot(2,1,2)
+plt.plot(rayleigh_bins_center,rayleigh_hist[0]/rayleigh_expected,'rs--')
+plt.xlabel('Bins'), plt.ylabel('Relative deviation'), plt.title('Rayleigh: Relative deviation to scipy')
+
+plt.figure(4)
+
+plt.subplot(2,1,1)
+plt.plot(laplace_bins_center,laplace_hist[0],'s--',laplace_bins_center,laplace_expected,'o:')
+plt.xlabel('Bins'), plt.ylabel('Count'), plt.title('Laplace: Distribution')
+plt.legend(['histogram gr::random','calculation scipy'],loc=1)
+
+plt.subplot(2,1,2)
+plt.plot(laplace_bins_center,laplace_hist[0]/laplace_expected,'rs--')
+plt.xlabel('Bins'), plt.ylabel('Relative deviation'), plt.title('Laplace: Relative deviation to scipy')
+
+plt.show()
diff --git a/gnuradio-runtime/include/gnuradio/random.h b/gnuradio-runtime/include/gnuradio/random.h
index e01fcb7be7..e95f36d1c2 100644
--- a/gnuradio-runtime/include/gnuradio/random.h
+++ b/gnuradio-runtime/include/gnuradio/random.h
@@ -1,6 +1,6 @@
 /* -*- c++ -*- */
 /*
- * Copyright 2002 Free Software Foundation, Inc.
+ * Copyright 2002, 2015 Free Software Foundation, Inc.
  *
  * This file is part of GNU Radio
  *
@@ -26,17 +26,9 @@
 #include <gnuradio/api.h>
 #include <gnuradio/gr_complex.h>
 
-// While rand(3) specifies RAND_MAX, random(3) says that the output
-// ranges from 0 to 2^31-1 but does not specify a macro to denote
-// this.  We define RANDOM_MAX for cleanliness.  We must omit the
-// definition for systems that have made the same choice.  (Note that
-// random(3) is from 4.2BSD, and not specified by POSIX.)
-
-#ifndef RANDOM_MAX
-static const int RANDOM_MAX = 2147483647; // 2^31-1
-#endif /* RANDOM_MAX */
-
 #include <stdlib.h>
+#include <boost/random.hpp>
+#include <ctime>
 
 namespace gr {
 
@@ -47,31 +39,51 @@ namespace gr {
   class GR_RUNTIME_API random
   {
   protected:
-    static const int NTAB = 32;
     long d_seed;
-    long d_iy;
-    long d_iv[NTAB];
-    int	d_iset;
-    float d_gset;
+    bool d_gauss_stored;
+    float d_gauss_value;
+
+    boost::mt19937 *d_rng; // mersenne twister as random number generator
+    boost::uniform_real<float> *d_uniform; // choose uniform distribution, default is [0,1)
+    boost::variate_generator<boost::mt19937&, boost::uniform_real<float> > *d_generator;
 
   public:
-    random(long seed=3021);
+    random(unsigned int seed=0);
+    ~random();
 
-    void reseed(long seed);
+    /*!
+     * \brief Change the seed for the initialized number generator. seed = 0 initializes the random number generator with the system time. Note that a fast initialization of various instances can result in the same seed.
+     */
+    void reseed(unsigned int seed);
 
     /*!
-     * \brief uniform random deviate in the range [0.0, 1.0)
+     * \brief Uniform random numbers in the range [0.0, 1.0)
      */
     float ran1();
 
     /*!
-     * \brief normally distributed deviate with zero mean and variance 1
+     * \brief Normally distributed random numbers (Gaussian distribution with zero mean and variance 1)
      */
     float gasdev();
 
+    /*!
+     * \brief Laplacian distributed random numbers with zero mean and variance 1
+     */
     float laplacian();
-    float impulse(float factor);
+
+    /*!
+     * \brief Rayleigh distributed random numbers (zero mean and variance 1 for the underlying Gaussian distributions)
+     */
     float rayleigh();
+
+    /*!
+     * \brief FIXME: add description
+     */
+    float impulse(float factor);
+
+    /*!
+     * \brief Normally distributed random numbers with zero mean and variance 1 on real and imaginary part. This results in a Rayleigh distribution for the amplitude and an uniform distribution for the phase.
+     */
     gr_complex rayleigh_complex();
   };
 
diff --git a/gnuradio-runtime/lib/math/random.cc b/gnuradio-runtime/lib/math/random.cc
index 7170f27edf..a2b2621307 100644
--- a/gnuradio-runtime/lib/math/random.cc
+++ b/gnuradio-runtime/lib/math/random.cc
@@ -1,6 +1,6 @@
 /* -*- c++ -*- */
 /*
- * Copyright 2002 Free Software Foundation, Inc.
+ * Copyright 2002, 2015 Free Software Foundation, Inc.
  *
  * This file is part of GNU Radio
  *
@@ -44,104 +44,75 @@
 
 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));
+    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);
+      }
   }
 
-  /*
-   * Copied from The KC7WW / OH2BNS Channel Simulator
-   * FIXME Need to check how good this is at some point
-   */
   float
   random::laplacian()
   {
-    float z = ran1();
-    if(z < 0.5)
-      return log(2.0 * z) / M_SQRT2;
-    else
-      return -log(2.0 * (1.0 - z)) / M_SQRT2;
+    float z = ran1()-0.5;
+    if(z>0) return -log(1-2*z);
+    else return log(1+2*z);
   }
 
   /*
@@ -159,26 +130,12 @@ namespace gr {
       return z;
   }
 
-  /*
-   * Complex rayleigh is really gaussian I and gaussian Q
-   * It can also be generated by real rayleigh magnitude and
-   * uniform random angle
-   * Adapted from The KC7WW / OH2BNS Channel Simulator
-   * FIXME Need to check how good this is at some point
-   */
   gr_complex
   random::rayleigh_complex()
   {
     return gr_complex(gasdev(),gasdev());
   }
 
-  /*   Other option
-       mag = rayleigh();
-       ang = 2.0 * M_PI * RNG();
-       *Rx = rxx * cos(z);
-       *Iy = rxx * sin(z);
-       */
-
   float
   random::rayleigh()
   {
diff --git a/gnuradio-runtime/python/gnuradio/gr/qa_random.py b/gnuradio-runtime/python/gnuradio/gr/qa_random.py
new file mode 100644
index 0000000000..83fee56181
--- /dev/null
+++ b/gnuradio-runtime/python/gnuradio/gr/qa_random.py
@@ -0,0 +1,64 @@
+#!/usr/bin/env python
+#
+# Copyright 2006,2007,2010,2015 Free Software Foundation, Inc.
+#
+# This file is part of GNU Radio
+#
+# GNU Radio is free software; you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation; either version 3, or (at your option)
+# any later version.
+#
+# GNU Radio is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with GNU Radio; see the file COPYING.  If not, write to
+# the Free Software Foundation, Inc., 51 Franklin Street,
+# Boston, MA 02110-1301, USA.
+#
+
+from gnuradio import gr, gr_unittest
+import numpy as np
+
+class test_random(gr_unittest.TestCase):
+
+    # NOTE: For tests on the output distribution of the random numbers, see gnuradio-runtime/apps/evaluation_random_numbers.py.
+
+    # Check for range [0,1) of uniform distributed random numbers
+    def test_1(self):
+        num_tests = 10000
+        values = np.zeros(num_tests)
+        rndm = gr.random()
+        for k in range(num_tests):
+            values[k] = rndm.ran1()
+        for value in values:
+            self.assertLess(value, 1)
+            self.assertGreaterEqual(value, 0)
+
+    # Check reseed method (init with time and seed as fix number)
+    def test_2(self):
+        num = 5
+
+        rndm0 = gr.random(42); # init with time
+        rndm1 = gr.random(42); # init with fix seed
+        for k in range(num):
+            x = rndm0.ran1();
+            y = rndm1.ran1();
+            self.assertEqual(x,y)
+
+        x = np.zeros(num)
+        y = np.zeros(num)
+        rndm0 = gr.random(42); # init with fix seed 1
+        for k in range(num):
+            x[k] = rndm0.ran1();
+        rndm1.reseed(43); # init with fix seed 2
+        for k in range(num):
+            y[k] = rndm0.ran1();
+        for k in range(num):
+            self.assertNotEqual(x[k],y[k])
+
+if __name__ == '__main__':
+    gr_unittest.run(test_random, "test_random.xml")
diff --git a/gnuradio-runtime/swig/runtime_swig.i b/gnuradio-runtime/swig/runtime_swig.i
index d4b55f134b..15cfebd300 100644
--- a/gnuradio-runtime/swig/runtime_swig.i
+++ b/gnuradio-runtime/swig/runtime_swig.i
@@ -60,6 +60,7 @@
 #include <gnuradio/top_block.h>
 #include <gnuradio/logger.h>
 #include <gnuradio/math.h>
+#include <gnuradio/random.h>
 %}
 
 %constant int sizeof_char 	= sizeof(char);
@@ -95,3 +96,4 @@
 %include "gr_ctrlport.i"
 
 %include "gnuradio/math.h"
+%include "gnuradio/random.h"
author	Johnathan Corgan <johnathan@corganlabs.com>	2015-09-05 18:05:13 -0700
committer	Johnathan Corgan <johnathan@corganlabs.com>	2015-09-05 18:05:13 -0700
commit	24c336cd99d8012681d40ccb00b1966e3c28450f (patch)
tree	45cf228c11a1669f0e62bee13e3efbad17fcd4bd /gnuradio-runtime
parent	273e5068d6d418eebb77dc0b5e15da20018bb90f (diff)
parent	45c0feeaeffde561758a4d92bf4463d60b85633a (diff)