// ------------------------------------------------------------------------
// file: pe_hdd_awgn.c
//
// Probability of decoding error (block error) for a binary code with
// binary transmission over an AWGN channel and 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 <stdio.h>
#include <math.h>

#define NMAX 64            // Maximum code length 

main(argc,argv)
int argc;
char **argv;
{
  int n, k;
  int t, dmin;
  double a[NMAX+1];        // Weight distribution
  double sum, p, Pe;
  double rate;
  int i, n_max, ki;

  int d_hi[64];
  double d_h[64],n_d[64];

  double prob_sphere(int n, int k, int j, double p);
  double fact(double a);
  double comb(double a, double b);
  double Q(double a);

  double init,final,inc,eb_no_db,es_no;

  char name1[80], name2[80];
  FILE *fp1,*fp2;

  if (argc != 9)
  {
    printf("Usage: %s n k d_min file_WD file_Pe init_SNR  final_SNR  inc_SNR\n",
    argv[0]);
    exit(1);
  }
  sscanf(argv[1], "%d", &n);
  sscanf(argv[2], "%d", &k);
  sscanf(argv[3], "%d", &dmin);
  sscanf(argv[4], "%s", name2);
  sscanf(argv[5], "%s", name1);
  sscanf(argv[6], "%lf", &init);
  sscanf(argv[7], "%lf", &final);
  sscanf(argv[8], "%lf", &inc);

  rate = (double) k / (double) n;
  t = floor(dmin-1)/2;

  fp1 = fopen(name1,"w");
  fp2 = fopen(name2,"r");

  n_max = 0;
  while(fscanf(fp2,"%d %lf\n",&d_hi[n_max],&n_d[n_max])!=EOF)
    {
      n_max++;
    }
  fclose(fp2);

  for (i=0;i<n_max;i++)
    d_h[i] = (double) d_hi[i];

  for (eb_no_db = init; eb_no_db<=final; eb_no_db+=inc)
    {
    es_no = pow(10.0,(eb_no_db/10.0));   /* snr per bit */
    es_no = es_no*rate;                  /* snr per symbol */

    p = Q(sqrt(2.0*es_no));

    Pe = 0.0;
    for (i=0; i<n_max; i++)
      {
      sum = 0.0;
      for (ki=0; ki<=t; ki++)
        sum += prob_sphere(n, ki, d_hi[i], p);
      Pe += ( n_d[i] * sum );
      }
    fprintf(fp1,"%lf %e\n", eb_no_db, Pe);
    }
}



// Probability of error patterns centered at a word of weight j, within
// radius k. For words of length n and probability of a bit error p.
//
double prob_sphere(int n, int k, int j, double p)
{
double sum;
double term1, term2, term3, term4;
double comb(double a, double b);
int r;
  
  sum = 0;
  for (r=0; r<=k; r++)
    {
    term1 = comb((double)     j, (double)(k-r));
    term2 = comb((double) (n-j),    (double) r);
    term3 = pow(      p,   (double)(j-k+2*r));
    term4 = pow((1.0-p), (double)(n-j+k-2*r));

// printf("n = %d, k = %d, j = %d, r = %d, j-k+2r = %d,  n-j+k-2r = %d\n", 
//        n, k, j, r, (j-k+2*r), (n-j+k-2*r));

    sum += term1 * term2 * term3 * term4;
    }
  return(sum);
}



double fact(double a)
{
double i,tot;
  
    tot = 1.0;
    for (i=a;i>0.001;i=i-1.0)
    {
        tot = tot * i;
    }
    return(tot);
}



double comb(double a,double b)
{
double z,tot;
 
    if (a<b)
       return(1.0);
    if (b > (a-b))
    {
       tot = 1.0/fact(a-b);
       for (z=a;z>b;z=z-1.0)
       {
          tot = tot * z;
       }
    }
    else
    {
       tot = 1.0/fact(b);
       for (z=a;z>a-b;z=z-1.0)
       {
          tot = tot * z;
       }
    }
    return(tot);
}



double Q(double a)
{
        double erfc( double x );
    return(0.5*erfc(a/sqrt(2.0)));
}