// ------------------------------------------------------------------------
// File: osd_order5.c
//
// Simulation of ordered statistics decoding, reprocessing order = 5
// ------------------------------------------------------------------------
// 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 <time.h>
#include "def.h"

#define MAXRAND LONG_MAX 

char name1[40], name2[40];
int G[K][N];
int info[N], code[N], deco[N];
double R[N], D[N];

int ord;

double INIT, FINAL, INC;
double RATE;
double snr;
double sigma;
long seed;
long num_sim, error, sim;

void channel();
void encoder();
void order5();


main(int argc, char *argv[])
{
  register int i,j;
  FILE *fp;

  // Command line processing
  if (argc != 8)
    {
      printf("Usage: %s file_matrix file_ber init_snr final_snr snr_inc num_sim seed\n", argv[0]);
    }

  sscanf(argv[1],"%s", name1);
  sscanf(argv[2],"%s", name2);
  sscanf(argv[3],"%lf", &INIT);
  sscanf(argv[4],"%lf", &FINAL);
  sscanf(argv[5],"%lf", &INC);
  sscanf(argv[6],"%ld",&num_sim);
  sscanf(argv[7],"%ld",&seed);

  RATE = (double)K/(double)N;

  printf("\nSimulation of soft-decision decoding with ordered statistics\n");
  printf("This is a (%d,%d,%d) linear code with generator matrix:\n",
            N, K, D_H);
  printf("(IMPORTANT NOTE: The matrix MUST be in systematic form)\n\n");

  fp = fopen(name1,"r");
  for (i=0; i<K; i++)
    for (j=0; j<N; j++)
      fscanf(fp,"%1d ", &G[i][j]);
  fclose(fp);

  for (i=0; i<K; i++) {
    for (j=0; j<N; j++)
      printf("%1d", G[i][j]); printf("\n"); }
  printf("\n");

  fp = fopen(name2,"w");

  snr = INIT;
  srandom(seed);

  // ========================== SNR LOOP ===============================

  while ( snr < (FINAL+1.0e-6) ) {

    // USE THIS FOR Es/No SIMULATIONS
    sigma = sqrt( 1.0 / ( 2.0 * pow(10.0, (snr/10.0)) ) ); 

    // UNCOMMENT THE LINE BELOW FOR Eb/No SIMULATIONS 
    //sigma = sqrt( 1.0 / ( 2.0 * RATE * pow(10.0, (snr/10.0)) ) ); 

    error = 0; 
    sim = 0; 

    // ####################  SIMULATION LOOP  #####################

    while (sim < num_sim) {
      sim++;
      for (i=0; i<K; i++) info[i] = (random()>>10)&01;

      encoder();
      channel();
      order5(G,R,D); 

      for (i=0; i<N; i++)
        if (D[i] == -1.0) deco[i] = 0;
        else              deco[i] = 1;

      for (i=0; i<N; i++)
        if (deco[i]^code[i]) error+=1;

//printf("%10ld %10ld\n", sim, error);

      }

    printf("%7.3lf  %10.7e\n", snr, ((double)error/(double)(N*num_sim)) );
    fprintf(fp,"%7.3lf  %10.7e\n", snr, ((double)error/(double)(N*num_sim)) );

    snr += INC;
  }

}



void encoder()
{
register int i,j;

  for (j=0; j<K; j++)
   code[j] = info[j];

  for (j=K; j<N; j++) {
    code[j] = 0;
    for (i=0; i<K; i++)
      code[j] ^= (info[i] & G[i][j]);
  }             
}




void channel() {
// -------------------------------------------------------------------------
//                 BPSK modulation and AWGN addition  
// -------------------------------------------------------------------------

  int i;
  double rndm, u1, u2, s, noise;

  for (i=0; i<N; i++) {
    do {
        rndm = (double)(random())/MAXRAND;
        u1 = rndm * 2 - 1.0;
        rndm = (double)(random())/MAXRAND;
        u2 = rndm * 2 - 1.0;
        s = u1 * u1 + u2 * u2;
      } while( s >= 1 );
    noise = u1 * sqrt( (-2.0*log(s))/s );
    if (code[i])
      R[i] = -1.0 + noise*sigma;
    else
      R[i] = +1.0 + noise*sigma;
  }

}