// ------------------------------------------------------------------------
// File: rm31_chase.c
//
// Simulation of RM(3,1), equivalent to the (8,4,4) extended Hamming code.
// Soft-decision decoding performed by Chase type-II algorithm.
// ------------------------------------------------------------------------
// 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>
#include <stdlib.h>

#define MAX_RANDOM LONG_MAX    // Maximum value of random() 
#define RATE 0.5               // Coding rate = 4/8 
#define INIT_SNR 3.0           // Initial value of Eb/N0 
#define FINAL_SNR  8.0         // Final value of Eb/N0 
#define SNR_INCREMENT 1.0      // Increment in Eb/N0 
#define NUMSIM 1000000        // Number of simulations (one per 4 bits) 
#define SWAP(a,b) itemp=(a); (a)=(b); (b)=itemp;

#define n 8 
#define k 4
#define nk 4                   // n-k

int L;

// Generator matric in systematic form
// Information positions: 1,2,3 and 4 (exchanged 4 and 5 in original)
//
int G[k][n] = { 1,0,0,0,1,1,1,0,
                0,1,0,0,1,1,0,1,
                0,0,1,0,1,0,1,1,
                0,0,0,1,0,1,1,1 };
// Original:
// int G[k][n] = { 1,0,0,1,0,1,1,0,
//                0,1,0,1,0,1,0,1,
//                0,0,1,1,0,0,1,1,
//                0,0,0,0,1,1,1,1 };



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      */



double sim, block_error;
double ber;
double amp;
double seed;
int error;
int data[k], codeword[n];
int data_int;
float snr;
float transmited[n];
float received[n];
int hard[n];
int decoded[n];


float reli[n];
float reli_ord[n+1];  // Indexed from 1 .. n
int perm[n+1];        // Indexed from 1 .. n


void initialize(void);
void awgn(void);
void encode(void);
void bpsk_dec(void);
void indexx(int ns, float arr[], int indx[]);
void dec2bin(int test, int err[]);
float metric(int vector[]);
void gen_err(int count[], int patt[]);



main()
{
  int i,j;
  unsigned long value, test;
  int err[n], patt[n];
  double max, temp;
  int count[n];
  int got_one;

  printf("Enter the value of L: "); scanf("%d", &L);

  // Compute 2^L - 1
  value = 1;
  for (i=0; i<L; i++) value = value << 1;
  value -= 1;

  printf("2^L-1 = %ld\n", value);

  snr = INIT_SNR;

  // *********** S I M U L A T I O N   L O O P  **************
  while ( snr < (FINAL_SNR+0.001) )
    {
      initialize();
      while (sim < NUMSIM) 
      { 

        // -----------        Random information
        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];

        // -----------        RM encode
        encode();

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

        // -----------        AWGN channel
        awgn();

        // -----------        Symbol-by-symbol decision --> h
        bpsk_dec();

        // -----------        Reliabilities of received symbols
        for (i=0; i<n; i++)
          reli[i] = fabs(received[i]);

        for (i=0; i<n; i++)
          reli_ord[i+1] = reli[i];

        // Order reliabilities in ascending order
        // -----------        Obtain permutation. NOTE: perm[1]..perm[n]
        indexx(n, reli_ord, perm);

#ifdef DEBUG
printf("prm = ");
for (i=0; i<n; i++) printf("%4d ",perm[i+1]);
printf("\n");
#endif

        test = 0;
        for (i=0; i<n; i++) patt[i] = hard[i];

        max = -DBL_MAX;      // For metric computation
        got_one = 0;          // Codeword found indication

        while (test <= value)
          {

#ifdef DEBUG
printf("test = %d pattern = ",test);
for (i=0; i<n;i++) printf("%4d ",patt[i]);
printf(" check = %d\n",check(patt));
#endif

            // -----------        Check that (h+e) is a codeword
            if (check(patt))
              {
              // -----------      Compute likelihood, compare & choose best
              got_one = 1;

#ifdef DEBUG
printf(" metric = %f\n", metric(patt));
#endif

              if ((temp=metric(patt)) > max)
                {
                  max = temp;
                  for (i=0; i<n; i++) decoded[i] = patt[i];
                }
              }

            // -----------      Increase counter over least L reliable positions
            test += 1;   // (Could do GRAY COUNTING here)
            dec2bin(test,count);

            // -----------      Generate error pattern, e
            gen_err(count,patt);

          }

        // ------------ Decode information
        i = 0;
        if (got_one)                   // If codeword was found
          for (j=0; j<k; j++)
            i = (i<<1) + decoded[j];
        else                           // else output hard decision
          for (j=0; j<k; j++)
            i = (i<<1) + hard[j];

        error = i ^ data_int;

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

        ber += (float) wh[error];

        sim+=1.0;
      }
    printf("%f %8.0f %8.0f %13.8e\n", snr, ber, (k*sim), (ber/(sim*k))); 
    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;
    }
#ifdef DEBUG
printf("\n-------\ncod = ");
for (j=0; j<n; j++)
printf("%4d ",codeword[j]);
printf("\n");
#endif
}








void bpsk_dec()
{
  int j;
  for (j=0; j<n; j++)
    {
    if (received[j] <= 0.0)
      hard[j] = 1;
    else
      hard[j] = 0;
    }
#ifdef DEBUG
printf("rec = ");
for (j=0; j<n; j++)
printf("%4.2lf ",received[j]);
printf("\n");
printf("hrd = ");
for (j=0; j<n; j++)
printf("%4d ",hard[j]);
printf("\n");
#endif
}











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 indexx(int ns, float arr[], int indx[])
{
// Indexes an array arr[1..n]: Outputs indx[1..n] such that
// arr[indx[j]] is in ascending order. Neither n nor arr are changed.
//
{
  int i, indxt, ir=ns, itemp, j,ks,l=1;
  int jstack = 0, istack[50];
  float a;
  for (i=1; i <= ns; i++) 
    indx[i] = i;
  for (;;)
    {
      if ( (ir-l) < 7) {
        for (j=l+1; j <= ir; j++)
          {
            indxt = indx[j];
            a = arr[indxt];
            for (i=j-1; i >= 1; i--)
              {
                if (arr[indx[i]] <= a) break;
                indx[i+1] = indx[i];
              }
            indx[i+1] = indxt;
          }
        if (jstack == 0) break;
        ir = istack[jstack--];
        l = istack[jstack--];
      }
      else {
        ks = (l+ir) >> 1;
        SWAP(indx[ks],indx[l+1]);
        if (arr[indx[l]]>arr[indx[ir]]) {
          SWAP(indx[l],indx[ir]);
	}
        if (arr[indx[l+1]]>arr[indx[ir]]) {
          SWAP(indx[l+1],indx[ir]);
	}
        if (arr[indx[l]]>arr[indx[l+1]]) {
          SWAP(indx[l],indx[l+1]);
	}
        i = l+1;
        j = ir;
        indxt = indx[l+1];
        a = arr[indxt];
        for (;;)
          {
            do i++; while (arr[indx[i]] < a);
            do j--; while (arr[indx[j]] > a);
            if (j < i) break;
            SWAP(indx[i],indx[j]);
          }
        indx[l+1] = indx[j];
        indx[j] = indxt;
        jstack += 2;
	if (jstack > 50) { printf("\nError in indxx\n"); exit(0); }
        if ((ir-i+1) >= j-1) {
          istack[jstack] = ir;
          istack[jstack-1] = i;
          ir = j-1;
	}
        else {
          istack[jstack] = j-1;
          istack[jstack-1] = l;
          l = i;
	}
      }
    }
  }
}




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




int check(int patt[])
{
// Compute syndromes
int i,j;
int syn[nk];
  for (i=0; i<nk; i++)
    {
      syn[i] = 0;
      for (j=0; j<k; j++)
        syn[i] ^= ( patt[j] * G[j][i+k] ) & 0x01;
      syn[i] ^= patt[k+i];
      if (syn[i])
        return(0);    // If syndrome nonzero test fails
    }
  return(1);
}



void dec2bin(int test, int count[])
{
int i;
  for (i=0; i<L; i++)
    {
      count[i] = test%2;
      test /=2;
    }
}




float metric(int vector[])
{
int i;
float temp;
  temp = 0.0;
  for (i=0; i<n; i++)
    temp += -(2.0*(float)vector[i]+1.0)*received[i];
  return(temp);
}





void gen_err(int count[], int patt[])
{
// Generate error pattern from vector representing the combination of
// least reliable positions
int i;
  for (i=0; i<n; i++) patt[i] = hard[i];
  for (i=0; i<L; i++)
    if (count[i])
      patt[perm[i+1]-1] ^= 1;

}