// ------------------------------------------------------------------------
// File:    633_awgn_hdd.c
// Author:  Robert Morelos-Zaragoza
// Date:    August 7, 2000
//
// Simulation of a binary linear (6,3,3) code with binary transmission
// over an AWGN 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[1024];
int hard[1024], corr[1024];
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);

// Codewords
//
int    codeword[8][6] = {  0,  0,  0,  0,  0,  0,                            
                           0,  0,  1,  1,  0,  1,                            
                           0,  1,  0,  0,  1,  1,                            
                           0,  1,  1,  1,  1,  0,                            
                           1,  0,  0,  1,  1,  0,                            
                           1,  0,  1,  0,  1,  1,                            
                           1,  1,  0,  1,  0,  1,                           
                           1,  1,  1,  0,  0,  0  };             

// Parity-check matrix
//
int H[3][6] = { 1, 0, 1, 1, 0, 0,
                1, 1, 0, 0, 1, 0,
                0, 1, 1, 0, 0, 1 };

// Look-up table (LUT):    syndrome --> error pattern
int   LUT[8][6] = { 0, 0, 0, 0, 0, 0,           // s[] = 000
                    0, 0, 0, 0, 0, 1,           // s[] = 001
                    0, 0, 0, 0, 1, 0,           // s[] = 010
                    0, 1, 0, 0, 0, 0,           // s[] = 011
                    0, 0, 0, 1, 0, 0,           // s[] = 100
                    0, 0, 1, 0, 0, 0,           // s[] = 101
                    1, 0, 0, 0, 0, 0,           // s[] = 110
                    1, 0, 0, 0, 0, 1 };         // s[] = 111

int syndrome[3];
int flag;

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 random bits

        	   bpsk_awgn();

                   // BPSK demodulation

                   for (i=0; i<6; i++)
                     if (recd[i] >= 0.0) hard[i] = 0;
                     else                hard[i] = 1;

// hard[0] ^= 0;
// hard[1] ^= 0;
// hard[2] ^= 0;
// hard[3] ^= 0;
// hard[4] ^= 0;
// hard[5] ^= 1;

		   decode();

// printf("codeword  = ");
// for (i=0; i<6; i++)
// printf("%d ", codeword[data][i]);
// printf("\n");
// printf("hard dec  = ");
// for (i=0; i<6; i++)
// printf("%d ", hard[i]);
// printf("\n");
// printf("corrected = ");
// for (i=0; i<6; i++)
// printf("%d ", corr[i]);
// printf("\n");

                   flag = 0;
                   for (i=0; i<6; i++)
                     if (codeword[data][i] != corr[i])
                       flag = 1;

		   if (flag)
		     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()
{
int i,j;
int index;

  // Compute syndrome
  for (i=0; i<3; i++)
    {
    syndrome[i] = 0;
    for (j=0; j<6; j++)
      syndrome[i] ^= (hard[j] & H[i][j]);
    }

  // Transform to an integer to index the LUT
  index = 0;
  for (i=0; i<3; i++)
    index = 2*index + syndrome[i];

// printf("index = %d\n", index);


  // Look-up error pattern and add it to the received vector
  for (i=0; i<6; i++)
    corr[i] = hard[i] ^ LUT[index][i];
}


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] = -(2.0*codeword[data][i]-1.0) + (noise/amp);
    }
}