// ------------------------------------------------------------------------
// File:    633_rayleigh.c
// Author:  Robert Morelos-Zaragoza
// Date:    May 14, 2001
//
// Simulation of a binary linear (6,3,3) code with binary transmission
// over an Rayleigh fading channel. HARD-DECISION decoding.
// ------------------------------------------------------------------------
// 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 <math.h>
#include <stdio.h>
#include <float.h>
#include <limits.h>
#include <stdlib.h>
#define MAX_RANDOM LONG_MAX    // Maximum value of random() 

int i;
int n, k;
double recd[6];
int data,est;
FILE *fp2;

double rate;
float init_snr;
float final_snr;
float snr_increment;
double snr;
double num_sim;
double sim;
double amp;
long seed;
double error;

char filename[40], name2[40];

void decode(void);
void bpsk_awgn(void);

double codeword[8][6] = { +1, +1, +1, +1, +1, +1,                // 0b000000 
                          +1, +1, -1, -1, +1, -1,                // 0b001101 
                          +1, -1, +1, +1, -1, -1,                // 0b010011 
                          +1, -1, -1, -1, -1, +1,                // 0b011110 
                          -1, +1, +1, -1, -1, +1,                // 0b100110 
                          -1, +1, -1, +1, -1, -1,                // 0b101011 
                          -1, -1, +1, -1, +1, -1,                // 0b110101 
                          -1, -1, -1, +1, +1, +1             };  // 0b111000 


main(int argc, char *argv[])
{
  	// Command line processing
  	if (argc != 7)
    	{
      	printf("Usage: %s init_snr final_snr snr_inc num_sim output_file seed\n", 
                      argv[0]);
	printf(" - init_snr is the initial value of Eb/No (dB)\n");
	printf(" - final_snr is the final value of Eb/No (dB)\n");
	printf(" - snr_inc is the increment in Eb/No (dB)\n");
	printf(" - num_sim is the number of simulations per Eb/No value\n");
	printf(" - output_file is the name of a file with Eb/No and BER\n");
	printf(" - seed is the value used in srandom\n");
      	exit(0);
    	}

  	sscanf(argv[1],"%f", &init_snr);
  	sscanf(argv[2],"%f", &final_snr);
  	sscanf(argv[3],"%f", &snr_increment);
  	sscanf(argv[4],"%lf",&num_sim);
  	sscanf(argv[5],"%s", name2);
  	sscanf(argv[6],"%lf",&seed);

  	fp2 = fopen(name2,"w");

	srandom(seed);

  	rate = 0.5;

	snr = init_snr;

	while ( snr < (final_snr+0.001) )
	  { 
  		amp = sqrt(2.0*rate*pow(10.0,(snr/10.0)));
  		error = 0.0;
  		sim = 0.0;

		while (sim < num_sim)
		{
		   data = ( random() >> 10) % 8;   // 3 bits

        	   bpsk_awgn();
		   decode();

		   if (data != est)
		     error += 1.0;
        	   sim += 1.0;
            }
    	  printf("%f %8.0f %13.8e\n", snr, error, (error/sim)); 
    	  fflush(stdout);
    	  fprintf(fp2, "%f %13.8e\n", snr, (error/sim) );
    	  fflush(fp2);
    	  snr += snr_increment;
	}
}




void decode()
{
double aux1, aux2;
double M1,M2,M3,M4,M12,M34;
int dec1,dec2,dec3,dec4,dec12,dec34;

  aux1 = recd[0] + recd[1];
  aux2 = recd[0] - recd[1];

  M1 = aux1 + recd[2];
  M2 = aux1 - recd[2];
  M3 = aux2 + recd[2];
  M4 = aux2 - recd[2];

  if (M1 >= 0.0) dec1 = 0;   else { dec1 = 7; M1 = -M1; };
  if (M2 >= 0.0) dec2 = 1;   else { dec2 = 6; M2 = -M2; };
  if (M3 >= 0.0) dec3 = 2;   else { dec3 = 5; M3 = -M3; };
  if (M4 >= 0.0) dec4 = 3;   else { dec4 = 4; M4 = -M4; };

  M1 +=   ( recd[3] + recd[4] + recd[5]);
  M2 +=   (-recd[3] + recd[4] - recd[5]);
  M3 +=   ( recd[3] - recd[4] - recd[5]);
  M4 +=   (-recd[3] - recd[4] + recd[5]);

  if (M1 > M2) 
    {
    dec12 = dec1;
    M12 = M1;
    }
  else 
    {
    dec12 = dec2;
    M12 = M2;
    }

  if (M3 > M4) 
    {
    dec34 = dec3;
    M34 = M3;
    }
  else 
    {
    dec34 = dec4;
    M34 = M4;
    }

  if (M12 > M34) est = dec12;   else est = dec34;

}



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())/MAX_RANDOM;
      u1 = rndm * 2.0 - 1.0;
      rndm = (double)(random())/MAX_RANDOM;
      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);
}


void bpsk_awgn()
//
// BPSK map, AWGN add and BPSK detect
// 
{
  double u1,u2,s,noise,randmum;
  int i;
  for (i=0; i<6; i++)
    	{
      		do {
            		randmum = (double)(random())/MAX_RANDOM;
            		u1 = randmum*2.0 - 1.0;
            		randmum = (double)(random())/MAX_RANDOM;
            		u2 = randmum*2.0 - 1.0;
            		s = u1*u1 + u2*u2;
            	} while( s >= 1);
      		noise = u1 * sqrt( (-2.0*log(s))/s );
      		recd[i] = rayleigh() * codeword[data][i] + (noise/amp);
    	}
}