bch_awgn.c

The intention is to show you how to incorporate the AWGN/Rayleigh fading models in the basic decodin

// ------------------------------------------------------------------------
// File:    bch_awgn.c
// Date:    January 12, 2000
//
// Simulation of hard-decision decoding of a binary BCH code with binary
// transmission over an AWGN channel. Decoding using the BM 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 <math.h>
#include <stdio.h>
#include <float.h>
#include <limits.h>
#include <stdlib.h>
#define MAX_RANDOM LONG_MAX    // Maximum value of random() 

int i;
int m, n, length, k, t, d;
int p[10];
int alpha_to[1024], index_of[1024], g[1024];
int recd[1024], data[1024], bb[1024];
int numerr, errpos[1024], decerror = 0;
FILE *fp2;

double rate;
float init_snr;
float final_snr;
float snr_increment;
double snr;
double num_sim;
double sim;
double ber;
double amp;
long seed;
int error;

char filename[40], name2[40];

void read_p(void);
void generate_gf(void);
void gen_poly(void);
void encode_bch(void);
void decode_bch(void);
void bpsk_awgn(void);

main(int argc, char *argv[])
{
  	// Command line processing
  	if (argc != 10)
    	{
	printf("\nSimulation of BCH codes over a binary input AWGN channel\n");
	printf("Copyright 1994-2000. Robert Morelos-Zaragoza\n\n");
      	printf("Usage: %s m length t init_snr final_snr snr_inc num_sim output_file seed\n", 
                      argv[0]);
	printf(" - m is the order of GF(2^m)\n");
	printf(" - length of the BCH code\n");
 	printf(" - t = [d-1/2]\n");
	printf(" - init_snr is the initial value of Eb/No (dB)\n");
	printf(" - final_snr is the final value of Eb/No (dB)\n");
	printf(" - snr_inc is the increment in Eb/No (dB)\n");
	printf(" - num_sim is the number of simulations per Eb/No value\n");
	printf(" - output_file is the name of a file with Eb/No and BER\n");
	printf(" - seed is the value used in srandom\n");
      	exit(0);
    	}

  	sscanf(argv[1],"%d", &m);
  	sscanf(argv[2],"%d", &length);
  	sscanf(argv[3],"%d", &t);
  	sscanf(argv[4],"%f", &init_snr);
  	sscanf(argv[5],"%f", &final_snr);
  	sscanf(argv[6],"%f", &snr_increment);
  	sscanf(argv[7],"%lf",&num_sim);
  	sscanf(argv[8],"%s", name2);
  	sscanf(argv[9],"%lf",&seed);

  	fp2 = fopen(name2,"w");

	read_p();        /* Read m */
	generate_gf();   /* Construct the Galois Field GF(2**m) */
	gen_poly();      /* Compute the generator polynomial of BCH code */
	srandom(seed);

  	rate = (float) k / (float) n;

	snr = init_snr;

	while ( snr < (final_snr+0.001) )
	{ 
  		amp = sqrt(2.0*rate*pow(10.0,(snr/10.0)));
  		ber = 0.0;
  		sim = 0.0;

		while (sim < num_sim)
		{

		/* Randomly generate DATA */
		for (i = 0; i < k; i++)
		data[i] = ( random() & 65536 ) >> 16;

		encode_bch();    /* encode data */
		for (i = 0; i < length - k; i++) recd[i] = bb[i];
		for (i = 0; i < k; i++)	recd[i + length - k] = data[i];

        	bpsk_awgn();
	
		decode_bch();	/* DECODE received codeword recv[] */

		// DECODING ERRORS? we compare only the data portion
        	decerror = 0;
		for (i = length - k; i < length; i++)
			if (data[i - length + k] != recd[i])
			decerror++;
        	ber += decerror;
        	sim += 1.0;
		}
    	printf("%f %8.0f %8.0f %13.8e\n", snr, ber, (k*sim), (ber/(sim*k))); 
    	fflush(stdout);
    	fprintf(fp2, "%f %13.8e\n", snr, (ber/(sim*k)) );
    	fflush(fp2);
    	snr += snr_increment;
	}
}


void 
read_p()
/*
 *	Read m, the degree of a primitive polynomial p(x) used to compute the
 *	Galois field GF(2**m). Get precomputed coefficients p[] of p(x). Read
 *	the code length.
 */
{
	int			i, ninf;

	printf("bch_awgn: Simulation of binary BCH codes over a binary input AWGN channel\n");
	printf("Copyright (c) 1994-2000. Robert Morelos-Zaragoza. All rights reserved.\n");
	printf("This program is free. Read the copyright notice.\n");
	for (i=1; i<m; i++)
		p[i] = 0;
	p[0] = p[m] = 1;
	if (m == 2)			p[1] = 1;
	else if (m == 3)	p[1] = 1;
	else if (m == 4)	p[1] = 1;
	else if (m == 5)	p[2] = 1;
	else if (m == 6)	p[1] = 1;
	else if (m == 7)	p[1] = 1;
	else if (m == 8)	p[4] = p[5] = p[6] = 1;
	else if (m == 9)	p[4] = 1;
	else if (m == 10)	p[3] = 1;
	else if (m == 11)	p[2] = 1;
	else if (m == 12)	p[3] = p[4] = p[7] = 1;
	else if (m == 13)	p[1] = p[3] = p[4] = 1;
	else if (m == 14)	p[1] = p[11] = p[12] = 1;
	else if (m == 15)	p[1] = 1;
	else if (m == 16)	p[2] = p[3] = p[5] = 1;
	else if (m == 17)	p[3] = 1;
	else if (m == 18)	p[7] = 1;
	else if (m == 19)	p[1] = p[5] = p[6] = 1;
	else if (m == 20)	p[3] = 1;
	printf("Primitive polynomial of the GF, p(x) = ");
    	n = 1;
	for (i = 0; i <= m; i++) {
        	n *= 2;
		printf("%1d", p[i]);
        }
	printf("\n");
	n = n / 2 - 1;

}


void 
generate_gf()
/*
 * Generate field GF(2**m) from the irreducible polynomial p(X) with
 * coefficients in p[0]..p[m].
 *
 * Lookup tables:
 *   index->polynomial form: alpha_to[] contains j=alpha^i;
 *   polynomial form -> index form:	index_of[j=alpha^i] = i
 *
 * alpha=2 is the primitive element of GF(2**m) 
 */
{
	register int    i, mask;

	mask = 1;
	alpha_to[m] = 0;
	for (i = 0; i < m; i++) {
		alpha_to[i] = mask;
		index_of[alpha_to[i]] = i;
		if (p[i] != 0)
			alpha_to[m] ^= mask;
		mask <<= 1;
	}
	index_of[alpha_to[m]] = m;
	mask >>= 1;
	for (i = m + 1; i < n; i++) {
		if (alpha_to[i - 1] >= mask)
		  alpha_to[i] = alpha_to[m] ^ ((alpha_to[i - 1] ^ mask) << 1);
		else
		  alpha_to[i] = alpha_to[i - 1] << 1;
		index_of[alpha_to[i]] = i;
	}
	index_of[0] = -1;
}


void 
gen_poly()
/*
 * Compute the generator polynomial of a binary BCH code. Fist generate the
 * cycle sets modulo 2**m - 1, cycle[][] =  (i, 2*i, 4*i, ..., 2^l*i). Then
 * determine those cycle sets that contain integers in the set of (d-1)
 * consecutive integers {1..(d-1)}. The generator polynomial is calculated
 * as the product of linear factors of the form (x+alpha^i), for every i in
 * the above cycle sets.
 */
{
	register int	ii, jj, ll, kaux;
	register int	test, aux, nocycles, root, noterms, rdncy;
	int             cycle[1024][21], size[1024], min[1024], zeros[1024];

	/* Generate cycle sets modulo n, n = 2**m - 1 */
	cycle[0][0] = 0;
	size[0] = 1;
	cycle[1][0] = 1;
	size[1] = 1;
	jj = 1;			/* cycle set index */
	if (m > 9)  {
		printf("Computing cycle sets modulo %d\n", n);
		printf("(This may take some time)...\n");
	}
	do {
		/* Generate the jj-th cycle set */
		ii = 0;
		do {
			ii++;
			cycle[jj][ii] = (cycle[jj][ii - 1] * 2) % n;
			size[jj]++;
			aux = (cycle[jj][ii] * 2) % n;
		} while (aux != cycle[jj][0]);
		/* Next cycle set representative */
		ll = 0;
		do {
			ll++;
			test = 0;
			for (ii = 1; ((ii <= jj) && (!test)); ii++)	
			/* Examine previous cycle sets */
			  for (kaux = 0; ((kaux < size[ii]) && (!test)); kaux++)
			     if (ll == cycle[ii][kaux])
			        test = 1;
		} while ((test) && (ll < (n - 1)));
		if (!(test)) {
			jj++;	/* next cycle set index */
			cycle[jj][0] = ll;
			size[jj] = 1;
		}
	} while (ll < (n - 1));
	nocycles = jj;		/* number of cycle sets modulo n */

	// printf("Enter the error correcting capability, t: ");
	// scanf("%d", &t);

	d = 2 * t + 1;

	/* Search for roots 1, 2, ..., d-1 in cycle sets */
	kaux = 0;
	rdncy = 0;