awgn.c

根据LDPC码中码子的构造方法中的PEG算法

#include <iostream.h>
#include <math.h>

#include "AWGN.h"

AWGN::AWGN(long seed, double sigma) {
  (*this).seed = seed;
  (*this).sigma = sigma;

  // be sure to prevent seed to be non-positive
  if (seed < 0) {
    seed = -(seed);
  }
  if (seed == 0) {
    seed = 1;
  }
    
  // load a shuffle table
  // (used to be sure that there is less correlation
  //  between sequential random numbers)
  for (int j = NTAB + 7; j >= 0; j--) {
    long k = (seed)/IQ;
    seed = IA*(seed - k*IQ) - IR*k;
    if (seed < 0) {
      seed += IM;
    }
    if (j < NTAB) {
      iv[j] = seed;
    }
  }
  iy = iv[0];
  iset = 0;
    
}

void AWGN::setSigma(double sigma) {
  (*this).sigma = sigma;
}

double AWGN::nextValue(void) {
  double v1, v2, fac, rsq;

  if (iset == 0) {
    do {
      v1 = 2.0*random(&seed) - 1.0;
      v2 = 2.0*random(&seed) - 1.0;
      rsq = v1*v1 + v2*v2;
    } while ((rsq >= 1.0) || (rsq == 0.0));
    fac = sqrt(-2.0*log(rsq)/rsq);
    gset = v1*fac;
    iset = 1;
    return (v2*fac*sigma);
  } else {
    iset = 0;
    return (gset*sigma);
  }
}

inline long AWGN::getSeed(void) {
  return (seed);
}

double AWGN::random(long *idum) {
  int j;
  long k;
  double temp;

  // compute idum by Schrage's method
  k = (*idum)/IQ;
  *idum = IA*(*idum - k*IQ) - IR*k;
  if (*idum < 0) {
    *idum += IM;
  }
  j = iy/NDIV;
  iy = iv[j];
  iv[j] = *idum;
  temp = iy*AM;
  if (temp > RNMX) {
    return (RNMX);
  } else {
    return (temp);
  }
}

void AWGN::test(void) {
  AWGN *awgn; // pointer to an AWGN object
  awgn = new AWGN(1, 1);  // instantiate a new AWGN object
  int count[100];
  double moment[9];
  int index;
  double value;
  double tmp;

  for (int i = 0; i < 100; i++) {
    count[i] = 0;
  }

  for (i = 0; i < 9; i++) {
    moment[i] = 0;
  }

  for (i = 0; i < NOV; i++) {
    value = (*awgn).nextValue();
    index = (int)(value*10) + 50;
    if (index < 0) {
      index = 0;
    }
    if (index > 99) {
      index = 99;
    }
    count[index]++;
    tmp = value;
    for (int j = 0; j < 9; j++) {
      moment[j] += tmp;
      tmp *= value;
    }
  }

  for (i = 0; i < 100; i++) {
    cout << count[i]/(1.0*NOV) << "\n";
  }

  for (i = 0; i < 9; i++) {
    cout << moment[i]/(1.0*NOV) << "\n";
  }

  delete awgn; // delete the AWGN object...
  awgn = NULL; // ... and set the pointer to NULL
}