pe_313_hd_awgn.c

error correction code

// ------------------------------------------------------------------------
// File: pe_313_hd_awgn.c
//
// Probability of error of a repetition code with hard-decision decoding
// over an AWGN channel
//
// ------------------------------------------------------------------------
// 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 1024          // Maximum code length 

main(argc,argv)
int argc;
char **argv;
{
  int n;
  int t;
  double a[NMAX+1];        // Weight distribution
  double p, Pc;
  int i;
  double nd, id;

  double fact(double a);
  double comb(double a, double b);

  double init,final,inc,  eb_no_db, es_no, rate;

  double Q(double a);
  double erfc( double x );

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

  if (argc != 4)
  {
    printf("Usage: %s n t file_Pe\n",
    argv[0]);
    exit(1);
  }
  sscanf(argv[1], "%d", &n);
  sscanf(argv[2], "%d", &t);
  sscanf(argv[3], "%s", name1);

  fp1 = fopen(name1,"w");

  rate = 1/(double)n;
  init = -5.0;
  final = 17.0;
  inc = 0.5;

  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));
      
      Pc = 0.0;
      for (i=0; i<=t; i++)
        {
        nd = (double) n;
        id = (double) i;
        Pc += comb(nd,id) * pow(p,id)*pow((1.0-p),(nd-id));
        }
      fprintf(fp1,"%lf\t%e\n", eb_no_db, 1.0-Pc);
    }

}

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