// ------------------------------------------------------------------------
// File: rm13_awgn.c
//
// Simulation of RM(1,3), equivalent to the (8,4,4) extended Hamming code.
// Soft-decision decoding performed by the Green machine 
// ------------------------------------------------------------------------
// 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 MAX_RANDOM LONG_MAX    /* Maximum value of random() */
#define RATE 0.5               /* Coding rate = 4/8 */
#define INIT_SNR 7.5           /* Initial value of Eb/N0 */
#define FINAL_SNR 8.0          /* Final value of Eb/N0 */
#define SNR_INCREMENT 0.5      /* Increment in Eb/N0 */
#define NUMSIM 50000000         /* Number of simulations (one per 4 bits) */

int wh[16] = { 0, 1, 1, 2, 1, 2, 2, 3,       /* Hamming weight function: */
               1, 2, 2, 3, 2, 3, 3, 4 };     /* wh[i] = weight of i      */
int n = 8; 
int k = 4;
int G[4][8] = { 1,1,1,1,1,1,1,1,
                0,1,0,1,0,1,0,1,
                0,0,1,1,0,0,1,1,
                0,0,0,0,1,1,1,1 };

float sim, block_error;
float ber;
float amp;
double seed;
int error;
int data[4], codeword[8];
int data_int;
float snr;
float transmited[8];
float received[8];
int estiword[8], estidata[4];

void initialize(void);
void awgn(void);
void encode(void);
void HD_decode(void);
int green_machine(void);

main()
{
  int i,j;

  snr = INIT_SNR;
  while ( snr < (FINAL_SNR+0.001) )
    {
      initialize();
      while (sim < NUMSIM) 
      { 
        for (i=0; i<k; i++)
          data[i] = (random()>>10) & 0x01;
        /* convert data[] to integer for error computation purposes */
        data_int = 0;
        for (i=0; i<k; i++)
          data_int = (data_int << 1) ^ data[i];

        encode();

        for (i=0; i<n; i++)
          if (codeword[i]) transmited[i] = -1.0;
          else transmited[i] = 1.0;

        awgn();

        /* HD_decode(); /* */

        i = green_machine();

        error = i ^ data_int;
        /* if (error) block_error+=1.0;    /* block error rate */

        ber += (float) wh[error];

        sim+=1.0;
      }
    printf("%f %13.8e\n", snr, (ber/(sim*4.0))); /* 1 sim = 4 bits! */
    fflush(stdout);
    snr += SNR_INCREMENT;
  }
}

void encode()
{
  int i,j;
  for (j=0; j<n; j++)
    {
    codeword[j] = 0;
    for (i=0; i<k; i++)
      codeword[j] ^= ( data[i] * G[i][j] ) & 0x01;
    }
}

void HD_decode()
{

}

int green_machine()
{
  float r1[8], r2[8], r3[8];
  int i, est_data;
  float max;

  /* First stage */
  r1[0] = received[0] + received[1];
  r1[1] = received[0] - received[1];
  r1[2] = received[2] + received[3];
  r1[3] = received[2] - received[3];
  r1[4] = received[4] + received[5];
  r1[5] = received[4] - received[5];
  r1[6] = received[6] + received[7];
  r1[7] = received[6] - received[7];

  /* Second stage */
  r2[0] = r1[0] + r1[2];
  r2[1] = r1[0] - r1[2];
  r2[2] = r1[1] + r1[3];
  r2[3] = r1[1] - r1[3];
  r2[4] = r1[4] + r1[6];
  r2[5] = r1[4] - r1[6];
  r2[6] = r1[5] + r1[7];
  r2[7] = r1[5] - r1[7];

  /* Third stage */
  r3[0] = r2[0] + r2[4];
  r3[1] = r2[0] - r2[4];
  r3[2] = r2[1] + r2[5];
  r3[3] = r2[1] - r2[5];
  r3[4] = r2[2] + r2[6];
  r3[5] = r2[2] - r2[6];
  r3[6] = r2[3] + r2[7];
  r3[7] = r2[3] - r2[7];

  /* Find coordinate with largest magnitude */
  max = 0.0;
  for (i=0; i<n; i++)
    if (fabs(r3[i]) > max)
      {
      est_data = i;
      max = fabs(r3[i]);
      }

  /* Decode the sign bit */
  if (r3[est_data] < 0.0) est_data += 8;

  return(est_data);
}

void awgn()
{
  double u1,u2,s,noise,randmum;
  int i;
  for (i=0; i<n; 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 );
      received[i] = transmited[i] + noise/amp;
#ifdef NO_NOISE
      received[i] = transmited[i];
#endif
    }
}

void initialize()
{
  time(&seed);
  srandom(seed);
  amp = sqrt(2.0*RATE*pow(10.0,(snr/10.0)));
  block_error = 0.0;
  ber = 0.0;
  sim = 0.0;
}