git.gag.com Git - debian/gnuradio/blob - gnuradio-core/src/lib/general/gr_random.cc

   1 /* -*- c++ -*- */
   2 /*
   3  * Copyright 2002 Free Software Foundation, Inc.
   4  *
   5  * This file is part of GNU Radio
   6  *
   7  * GNU Radio is free software; you can redistribute it and/or modify
   8  * it under the terms of the GNU General Public License as published by
   9  * the Free Software Foundation; either version 2, or (at your option)
  10  * any later version.
  11  *
  12  * GNU Radio is distributed in the hope that it will be useful,
  13  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  15  * GNU General Public License for more details.
  16  *
  17  * You should have received a copy of the GNU General Public License
  18  * along with GNU Radio; see the file COPYING.  If not, write to
  19  * the Free Software Foundation, Inc., 51 Franklin Street,
  20  * Boston, MA 02110-1301, USA.
  21  */
  22
  23 /*
  24  *  Copyright 1997 Massachusetts Institute of Technology
  25  *
  26  *  Permission to use, copy, modify, distribute, and sell this software and its
  27  *  documentation for any purpose is hereby granted without fee, provided that
  28  *  the above copyright notice appear in all copies and that both that
  29  *  copyright notice and this permission notice appear in supporting
  30  *  documentation, and that the name of M.I.T. not be used in advertising or
  31  *  publicity pertaining to distribution of the software without specific,
  32  *  written prior permission.  M.I.T. makes no representations about the
  33  *  suitability of this software for any purpose.  It is provided "as is"
  34  *  without express or implied warranty.
  35  *
  36  */
  37
  38 #include <math.h>
  39 #include <gr_random.h>
  40
  41 #define IA 16807
  42 #define IM 2147483647
  43 #define AM (1.0/IM)
  44 #define IQ 127773
  45 #define IR 2836
  46 #define NDIV (1+(IM-1)/NTAB)
  47 #define EPS 1.2e-7
  48 #define RNMX (1.0-EPS)
  49
  50
  51 gr_random::gr_random (long seed)
  52 {
  53   reseed (seed);
  54 }
  55
  56 void
  57 gr_random::reseed (long seed)
  58 {
  59   d_seed = seed;
  60   d_iy = 0;
  61   for (int i = 0; i < NTAB; i++)
  62     d_iv[i] = 0;
  63   d_iset = 0;
  64   d_gset = 0;
  65 }
  66
  67 /*
  68  * This looks like it returns a uniform random deviate between 0.0 and 1.0
  69  * It looks similar to code from "Numerical Recipes in C".
  70  */
  71 float gr_random::ran1()
  72 {
  73   int j;
  74   long k;
  75   float temp;
  76
  77   if (d_seed <= 0 || !d_iy)  {
  78     if (-d_seed < 1)
  79       d_seed=1;
  80     else
  81       d_seed = -d_seed;
  82     for (j=NTAB+7;j>=0;j--) {
  83       k=d_seed/IQ;
  84       d_seed=IA*(d_seed-k*IQ)-IR*k;
  85       if (d_seed < 0)
  86         d_seed += IM;
  87       if (j < NTAB)
  88         d_iv[j] = d_seed;
  89     }
  90     d_iy=d_iv[0];
  91   }
  92   k=(d_seed)/IQ;
  93   d_seed=IA*(d_seed-k*IQ)-IR*k;
  94   if (d_seed < 0)
  95     d_seed += IM;
  96   j=d_iy/NDIV;
  97   d_iy=d_iv[j];
  98   d_iv[j] = d_seed;
  99   temp=AM * d_iy;
 100   if (temp > RNMX)
 101     temp = RNMX;
 102   return temp;
 103 }
 104
 105 /*
 106  * Returns a normally distributed deviate with zero mean and variance 1.
 107  * Also looks like it's from "Numerical Recipes in C".
 108  */
 109 float gr_random::gasdev()
 110 {
 111   float fac,rsq,v1,v2;
 112   d_iset = 1 - d_iset;
 113   if  (d_iset) {
 114     do {
 115       v1=2.0*ran1()-1.0;
 116       v2=2.0*ran1()-1.0;
 117       rsq=v1*v1+v2*v2;
 118     } while (rsq >= 1.0 || rsq == 0.0);
 119     fac= sqrt(-2.0*log(rsq)/rsq);
 120     d_gset=v1*fac;
 121     return v2*fac;
 122   }
 123   return d_gset;
 124 }
 125
 126 /*
 127  * Copied from The KC7WW / OH2BNS Channel Simulator
 128  * FIXME Need to check how good this is at some point
 129  */
 130
 131 float gr_random::laplacian()
 132 {
 133   float z = ran1();
 134   if (z < 0.5)
 135     return log(2.0 * z) / M_SQRT2;
 136   else
 137     return -log(2.0 * (1.0 - z)) / M_SQRT2;
 138 }
 139
 140 /*
 141  * Copied from The KC7WW / OH2BNS Channel Simulator
 142  * FIXME Need to check how good this is at some point
 143  */
 144
 145   // 5 => scratchy, 8 => Geiger
 146
 147 float gr_random::impulse(float factor = 5)
 148 {
 149   float z = -M_SQRT2 * log(ran1());
 150   if (fabsf(z) <= factor)
 151     return 0.0;
 152   else
 153     return z;
 154 }
 155
 156 /*
 157  * Complex rayleigh is really gaussian I and gaussian Q
 158  * It can also be generated by real rayleigh magnitude and
 159  * uniform random angle
 160  * Adapted from The KC7WW / OH2BNS Channel Simulator
 161  * FIXME Need to check how good this is at some point
 162  */
 163
 164 gr_complex gr_random::rayleigh_complex()
 165 {
 166   return gr_complex(gasdev(),gasdev());
 167 }
 168
 169 /*   Other option
 170      mag = rayleigh();
 171      ang = 2.0 * M_PI * RNG();
 172      *Rx = rxx * cos(z);
 173      *Iy = rxx * sin(z);
 174 */
 175
 176
 177 float gr_random::rayleigh()
 178 {
 179   return sqrt(-2.0 * log(ran1()));
 180 }