wer_mc_union.c

error correction code

// ------------------------------------------------------------------------
// File: WER_MC_union.c
//
// UNION BOUND on the error performance of a linear code with
// BPSK modulation over a flat Rayleigh fading channel. NO CSI.
//
// The expected value of the bit error probability is computed by
// Monte Carlo simulation.
// ------------------------------------------------------------------------
// This program is complementary material for the book:
//
// R.H. Morelos-Zaragoza, The Art of Error Correcting Coding, Wiley, 2002.
//
// ISBN 0471 49581 6
//
// This and other programs are available at http://the-art-of-ecc.com
//
// You may use this program for academic and personal purposes only. 
// If this program is used to perform simulations whose results are 
// published in a journal or book, please refer to the book above.
//
// The use of this program in a commercial product requires explicit
// written permission from the author. The author is not responsible or 
// liable for damage or loss that may be caused by the use of this program. 
//
// Copyright (c) 2002. Robert H. Morelos-Zaragoza. All rights reserved.
// ------------------------------------------------------------------------

#include <stdio.h>
#include <math.h>
#include <float.h>
#include <limits.h>

#define  MAXRAND     LONG_MAX                // for random number generation
//#define  ITERATIONS  100000                  // iterations for expected value

double Q(double x);
double rayleigh( void );

main(int argc, char *argv[])
{

  int    iter,i,i1,j,n,k,dimension,n_max,level;
  double ii1,init_snr,final_snr,inc_snr,eb_no_db,es_no,P1,P2;
  double sum,sum1,sum2,sumj,aux,beta,rate,c1,jaux,d12;
  double seed, temp;
  long int ITERATIONS;
   
  int d_hi[128];
  double d_h[128],n_d[128],fade[128];
  char name1[40],name2[40];
   
  double delta1,delta2,average,sed,prod;
   
  FILE *fp,*fp1;

  // Command line processing
  if (argc != 9)
    {
      printf("Usage: %s length(n) dimension(k) file_weight_dist bound_file init_snr final_snr inc_snr ITERATIONS \n", argv[0]);
      exit(0);
    }

  sscanf(argv[1],"%d", &n);
  sscanf(argv[2],"%d", &k);
  sscanf(argv[3],"%s", name1);
  sscanf(argv[4],"%s", name2);
  sscanf(argv[5],"%lf", &init_snr);
  sscanf(argv[6],"%lf", &final_snr);
  sscanf(argv[7],"%lf", &inc_snr);
  sscanf(argv[8],"%ld", &ITERATIONS);

  rate = (float) k / (float) n;

  time(&seed);
  srandom(seed);
  
  fp1 = fopen(name1,"r");

  n_max = 0;
  while(fscanf(fp1,"%d %lf\n",&d_hi[n_max],&n_d[n_max])!=EOF)
    n_max++;
  fclose(fp1);

  printf("Weight distribution:\n");
  for (i=0; i<n_max; i++)
    printf("%d %lf\n", d_hi[i], n_d[i]);
  printf("\n");

  fp = fopen(name2,"w");
  for (i=0;i<n_max;i++)
    d_h[i] = (double) d_hi[i];

  // Union bound without CSI
  for (eb_no_db = init_snr; eb_no_db<=final_snr; eb_no_db+=inc_snr)
    {
      es_no = pow(10.0,(eb_no_db/10.0));   /* snr per bit */
      es_no = es_no*rate;                  /* snr per symbol */
      sum = 0.0;
      for (i=0;i<n_max;i++)
      {
        average = 0.0;
        for (iter = 0; iter<ITERATIONS; iter++)
        {
          sum1 = 0.0;
          for (j=0; j<d_hi[i]; j++)
            {
            temp = rayleigh();
            // sum1 += (temp*temp);
            sum1 += rayleigh();
            }
          // average += Q ( sqrt(2.0*es_no*sum1) );
          average += Q ( sqrt(2.0*es_no*sum1*sum1/d_h[i]) );
        }
        sum += ( n_d[i] * (average/ITERATIONS) );
      }
      printf("%f\t%e\n",eb_no_db,sum);
      fprintf(fp, "%f\t%e\n",eb_no_db,sum);
      fflush(stdout);
      fflush(fp);
    }
 }

double rayleigh(void)
//
// Generate a Rayleigh distributed random number X with E{X^2}=1
//
{
  double rndm, u1, u2, s, x1, x2, aux, br;
  do
    {
      rndm = (double)(random())/MAXRAND;
      u1 = rndm * 2.0 - 1.0;
      rndm = (double)(random())/MAXRAND;
      u2 = rndm * 2.0 - 1.0;
      s = u1 * u1 + u2 * u2;
    } while( s >= 1.0 );
  x1 = u1 * M_SQRT2 * sqrt( (-log(s))/s );  /* Gaussian E{x1^2} = 1 */
  x2 = u2 * M_SQRT2 * sqrt( (-log(s))/s );  /* Gaussian E{x2^2} = 1 */
  // Rayleigh E{br^2} = 1
  br =  sqrt(x1*x1 + x2*x2) / M_SQRT2;
  return(br);
}

double Q(double a)
{
  double erfc( double x );
  return(0.5*erfc(a/sqrt(2.0)));
}