gr-osmosdr/lib/runtime/math/random.cc

/* -*- c++ -*- */
/*
 * Copyright 2002 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.
 */

/*
 *  Copyright 1997 Massachusetts Institute of Technology
 *
 *  Permission to use, copy, modify, distribute, and sell this software and its
 *  documentation for any purpose is hereby granted without fee, provided that
 *  the above copyright notice appear in all copies and that both that
 *  copyright notice and this permission notice appear in supporting
 *  documentation, and that the name of M.I.T. not be used in advertising or
 *  publicity pertaining to distribution of the software without specific,
 *  written prior permission.  M.I.T. makes no representations about the
 *  suitability of this software for any purpose.  It is provided "as is"
 *  without express or implied warranty.
 *
 */

#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include <math.h>
#include <gnuradio/random.h>

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(long seed)
  {
    reseed(seed);
  }

  void
  random::reseed(long seed)
  {
    d_seed = seed;
    d_iy = 0;
    for(int i = 0; i < NTAB; i++)
      d_iv[i] = 0;
    d_iset = 0;
    d_gset = 0;
  }

  /*
   * 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".
   */
  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;
  }

  /*
   * Returns a normally distributed deviate with zero mean and variance 1.
   * Also looks like it's from "Numerical Recipes in C".
   */
  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;
  }

  /*
   * 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;
  }

  /*
   * Copied from The KC7WW / OH2BNS Channel Simulator
   * FIXME Need to check how good this is at some point
   */
  // 5 => scratchy, 8 => Geiger
  float
  random::impulse(float factor = 5)
  {
    float z = -M_SQRT2 * log(ran1());
    if(fabsf(z) <= factor)
      return 0.0;
    else
      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()
  {
    return sqrt(-2.0 * log(ran1()));
  }

} /* namespace gr */