// ------------------------------------------------------------------------
//        File: qpsk_rm31.c
//        Date: April 2, 2002
// Description: Simulation of coded QPSK with an extended Hamming (8,4,4) code
// ------------------------------------------------------------------------
// 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 explicitely 
// 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   0.0         /* Initial value of Eb/N0 */
#define FINAL_SNR 10.0         /* Final value of Eb/N0 */
#define SNR_INCREMENT 1.0      /* Increment in Eb/N0 */
#define NUMSIM  100000         /* 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 };

double sim, block_error;
double ber;
double amp;
double seed;
int error;
int data[4], codeword[8];
int data_int;
double snr;
double transmitted_I[4],transmitted_Q[4];
double received_I[4],received_Q[4];

double max, metric;
int winner;

void initialize(void);
double awgn(void);
void encode(void);

main()
{
  int i,j,code_index[4];

  double qpsk_I[4] = { M_SQRT1_2,  M_SQRT1_2, -M_SQRT1_2, -M_SQRT1_2};
  double qpsk_Q[4] = { M_SQRT1_2, -M_SQRT1_2,  M_SQRT1_2, -M_SQRT1_2};

  int message[16][4] = {
    0,0,0,0,
    1,0,0,0,
    0,1,0,0,
    0,0,1,0,
    0,0,0,1,
    1,1,0,0,
    1,0,1,0,
    1,0,0,1,
    0,1,1,0,
    0,1,0,1,
    0,0,1,1,
    1,1,1,0,
    1,1,0,1,
    1,0,1,1,
    0,1,1,1,
    1,1,1,1,
  };

  double reference_I[16][4] = { 
    qpsk_I[0], qpsk_I[0], qpsk_I[0], qpsk_I[0],
    qpsk_I[3], qpsk_I[3], qpsk_I[3], qpsk_I[3],
    qpsk_I[2], qpsk_I[2], qpsk_I[2], qpsk_I[2],
    qpsk_I[0], qpsk_I[3], qpsk_I[0], qpsk_I[3],
    qpsk_I[0], qpsk_I[0], qpsk_I[3], qpsk_I[3],
    qpsk_I[1], qpsk_I[1], qpsk_I[1], qpsk_I[1],
    qpsk_I[3], qpsk_I[0], qpsk_I[3], qpsk_I[0],
    qpsk_I[3], qpsk_I[3], qpsk_I[0], qpsk_I[0],
    qpsk_I[2], qpsk_I[1], qpsk_I[2], qpsk_I[1],
    qpsk_I[2], qpsk_I[2], qpsk_I[1], qpsk_I[1],
    qpsk_I[0], qpsk_I[3], qpsk_I[3], qpsk_I[0],
    qpsk_I[1], qpsk_I[2], qpsk_I[1], qpsk_I[2],
    qpsk_I[1], qpsk_I[1], qpsk_I[2], qpsk_I[2],
    qpsk_I[3], qpsk_I[0], qpsk_I[0], qpsk_I[3],
    qpsk_I[2], qpsk_I[1], qpsk_I[1], qpsk_I[2],
    qpsk_I[1], qpsk_I[2], qpsk_I[2], qpsk_I[1] 
  };

  double reference_Q[16][4] = { 
    qpsk_Q[0], qpsk_Q[0], qpsk_Q[0], qpsk_Q[0],
    qpsk_Q[3], qpsk_Q[3], qpsk_Q[3], qpsk_Q[3],
    qpsk_Q[2], qpsk_Q[2], qpsk_Q[2], qpsk_Q[2],
    qpsk_Q[0], qpsk_Q[3], qpsk_Q[0], qpsk_Q[3],
    qpsk_Q[0], qpsk_Q[0], qpsk_Q[3], qpsk_Q[3],
    qpsk_Q[1], qpsk_Q[1], qpsk_Q[1], qpsk_Q[1],
    qpsk_Q[3], qpsk_Q[0], qpsk_Q[3], qpsk_Q[0],
    qpsk_Q[3], qpsk_Q[3], qpsk_Q[0], qpsk_Q[0],
    qpsk_Q[2], qpsk_Q[1], qpsk_Q[2], qpsk_Q[1],
    qpsk_Q[2], qpsk_Q[2], qpsk_Q[1], qpsk_Q[1],
    qpsk_Q[0], qpsk_Q[3], qpsk_Q[3], qpsk_Q[0],
    qpsk_Q[1], qpsk_Q[2], qpsk_Q[1], qpsk_Q[2],
    qpsk_Q[1], qpsk_Q[1], qpsk_Q[2], qpsk_Q[2],
    qpsk_Q[3], qpsk_Q[0], qpsk_Q[0], qpsk_Q[3],
    qpsk_Q[2], qpsk_Q[1], qpsk_Q[1], qpsk_Q[2],
    qpsk_Q[1], qpsk_Q[2], qpsk_Q[2], qpsk_Q[1] 
  };

  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();

        code_index[0] = codeword[1]*2 + codeword[0];
        code_index[1] = codeword[3]*2 + codeword[2];
        code_index[2] = codeword[5]*2 + codeword[4];
        code_index[3] = codeword[7]*2 + codeword[6];

        for (i=0; i<4; i++) {
          transmitted_I[i] = qpsk_I[code_index[i]];
          transmitted_Q[i] = qpsk_Q[code_index[i]];

          received_I[i] = transmitted_I[i] + awgn();
          received_Q[i] = transmitted_Q[i] + awgn();
        }

        // Decoding

        max = -999999.0;
       
        for (i=0; i<16; i++) {
          metric = 0;
          for (j=0;j<4;j++) {
            metric += (received_I[j] * reference_I[i][j]);
            metric += (received_Q[j] * reference_Q[i][j]);
          }

          if (metric > max) {
            max = metric;
            winner = i;
          }
        }

	i =  (message[winner][3]   + message[winner][2]*2);
	i += (message[winner][1]*4 + message[winner][0]*8);

        error = i ^ data_int;

        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;
    }
}


double awgn()
{
  double u1,u2,s,noise,randmum;
  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 ))/(amp*sqrt(2.0));
  return(noise);
}

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;
}