/***********************************************************************************  	FILENAME: awgn_simulator.c	version 1.5	(C) Tadashi Wadayama, 2003        An LDPC code simulator based on quantized-BP algorithm for AWGN chanel*********************************************************************************        ChangeLog:	Aug. 14, 2003 : start	Aug. 17, 2003 : min-sum algorithm	Aug. 18, 2003 : encoding	Aug. 23, 2003 : output format change, bug fix of str[], addition of  "puncture" 	Aug. 23, 2003 : range of "loop", tie break rule, check node update, interger-BP	Aug. 24, 2003 : comments, error handling in read_config_file()	Aug. 24, 2003 : TEST1,TEST2********************************************************************************** How to compile*********************************************************************************	gcc -O2 -o awgn_simulator awgn_simulator.c -lm	with -D TEST1: test for quantized_BP()	with -D TEST2: test for simulation with encoding**********************************************************************************  How to run*********************************************************************************    awgn_simulator filename_of_configuration_file snr    Example: awgn_simulator config 1.0**********************************************************************************  Configuation file format*********************************************************************************Example:--------------------------------------------------------matrix_file 	9974.5000 code_rate	0.4987quantize_depth	10max_iteration	100random_seed	1display		1stop_condition 	1stop_errors	100target_prob	1e-5error_logging	0puncture	0multi_edge_type	0min-sum		0norm_factor	0.0parity_check	1--------------------------------------------------------NOTE:	matrix_file 	: low density parity check matrix in spmat form	code_rate	: code rate of target code	quantize_depth	: quantization depth (recommendation: 8--13 bits )			  decoding performance depends on the depth	max_itreation	: maximum number of iterations in BP 	random_seed	: seed of radom number generator	display		: display mode of simulation status during simulation	                  0: no display, 1: full display, 2: simple display	stop_condition 	: conditions for  termination of simulation                          0: bit error, 1: block errors 			  if stop_condition >=2, simulation will terminates after			  stop_condition-loops.	stop_errors	: number of errors enough to stop simulation	target_prob	: target error probability you wish to obtain	error_logging	: error logging switch	puncture	: switch for puncture	                  0: no puncture, 1:puncture			  if puncture == 1, spmat file must contain			  puncture pattern after matrix definition	multi_edge_type	: switch for multi-edge type	                  usullaly, set to zero.	min-sum		: switch for BP algorithm                          0:sum-product, 1:min-sum	norm_factor	: normalization factor for min-sum algorithm	parity_check	: parity ckeck switch	                0: no parity check in BP, 1:perform parity check in BP**********************************************************************************  simulation status and report format (See also report_result2(), report_result())*********************************************************************************snr pb pB #ebits #bits #eblks #blks aveitr var sdv mis seed mitr n m r depth DLETA range config scond serr target log punc met msum norm pchk	snr   : signal to noise ratio Eb/N0	pb    : bit error probability after decoding 	pB    : block error probability after decoding 	#ebits: number of error bits	#bits : number of total bits generated	#eblks: number of error blocks	#blks : number of total blocks generated	aveitr: average number of iterations per block	var   : variance of Gaussian noise	sdv   : standard deviation of Gaussian noise	mis   : number of miscorrected blocks	seed  : seed of random number generator	mitr  : maximum number of iterations	n     : code length	m     : number of parity constraints	r     : code rate	depth : bit depth of quantization	DELTA : quantization step	range : range of quantized messsage	config: file name of configuration file	scond : 	        0->simulation stops when #ebits reaches to  serr	        1->simulation stops when #eblks reaches to  serr	serr  : number of errors enough to stop a simulation	target: target error probability	log   : logging for error blocks	punc  : puncture	met   : multi-edge type	msum  : min-sum algorithm	norm  : normalization factor	pchk  : parity check in BP algorithm**********************************************************************************"spmat form" for sparce parity check matrix*********************************************************************************Example (Do not include "#...." in a spmat file.)16 12	                          # code_length number_of_parity_check_constraint4 3                               # largest_row_weight largest_column_weight4 4 4 4 4 4 4 4 4 4 4 4           # row_weight_profile3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3   # column_weight_profile3 8 10 13                         # parity check matrix(caution:colum # started from 1!)4 7 9 132 5 7 104 6 11 143 9 15 161 6 9 104 8 12 152 6 12 161 7 14 163 5 12 142 11 13 151 5 8 11                    The above spmat format corresponds to the following matrix:164 0 0 1 0 0 0 0 1 0 1 0 0 1 0 0 0 0 0 0 1 0 0 1 0 1 0 0 0 1 0 0 0 0 1 0 0 1 0 1 0 0 1 0 0 0 0 0 0 0 0 0 1 0 1 0 0 0 0 1 0 0 1 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 1 1 0 0 0 0 1 0 0 1 1 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 1 0 0 1 0 0 1 0 0 0 1 0 0 0 0 0 1 0 0 0 1 1 0 0 0 0 0 1 0 0 0 0 0 0 1 0 1 0 0 1 0 1 0 0 0 0 0 0 1 0 1 0 0 0 1 0 0 0 0 0 0 0 0 1 0 1 0 1 0 1 0 0 0 1 0 0 1 0 0 1 0 0 0 0 0If puncture = 1, append puncture pattern like 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0to the bottom of a spmat file. In the above example, 1st, 2nd and 3rd bit will be punctured.*//**********************************************************************************/// Inlucde files and macros/**********************************************************************************/#include <stdio.h>#include <stdlib.h>#include <math.h>#define OFF 0#define ON 1#define BIT 0#define BLOCK 1#define MAX_STR_LEN 100/**********************************************************************************/// Global variables/**********************************************************************************///-------------------------------------for quantized_BP()int *tmp_bit;		// tentetative decision bitint **alpha;		// check to variable messageint **beta;		// variable to check messageint n;			// code lengthint m;			// number of constraintsint rmax;		// maximum number of row weightint cmax;		// maximum number of column weightint *row_weight;	// row weightsint *col_weight;	// column weightsint **row_list;		// row list of parity check matrixint **col_list_r;	// column list 1int **col_list_c;	// column list 1int LB,UB,LB2,UB2;	// upper and lower bound for quantized messagesint ASIZE;		// number of quantization levelsdouble DELTA;		// quantization step sizeint **R_tbl;		// table for check node operationdouble *LLR;		// array for received LLR valuesint *qLLR;		// quantized LLR valuesint *forward;		// forward array for check node operationint *backward;		// backward array for check node operationint max_itr;		// maximum number of iterations in BP//-------------------------------------for simulation_loop()double code_rate;	// code rateint total_iterations;	// total number of iterationsint seed;		// seed for random number generatorint quantization_level;	// quantization depthchar *filename;		// filename of configuration fileint total_blocks;	// total number of tested blocksint total_bits;		// total number of tested bitsint error_blocks;	// number of error blocksint error_bits;		// number of error bitsint total_miscorrection; // number of mis corrected blocksint display_result;	 // display switchint stop_condition;	 // stop conditionint stop_errors;	 // number of errors for stop conditionint error_logging;	 // error logging switchdouble target_error_prob; // target error probabilityint *punc_ptn;		 // puncture patternint puncture;		 // puncture switchint met;		 // multi-edge type switchint *vnode_type;	 // variable node typeint num_vnode_type;	 // number of variable node typeFILE *error_log;	 // error log fileint *vnode_errors;	 // number of variable node errorsint minsum;		 // min-sum algorithm switchdouble norm_factor;	 // normalization factorint parity_check;	 // parity check switchint encoding;//-------------------------------------for encodingint e_n;                // code lengthint e_m;		// number of parity symbolsint e_rmax;		// maximum number of row weightint e_cmax;		// maximum number of column weightint *e_row_weight;	// row weightsint *e_col_weight;	// column weightsint **e_row_list;	// parity check matrix//-------------------------------------channel dependent parametersint *cword;             // sent codeworddouble *rword;		// received word at receiverdouble channel_var;     // channel variancedouble snr;		// signal to noise ratio (Eb/N0)/**********************************************************************************/// encoding functions (used only for TEST2)/**********************************************************************************/void read_encfile(FILE *fp){  int i,j;  int v;  fscanf(fp,"%d %d\n",&e_n,&e_m);  fscanf(fp,"%d %d\n",&e_rmax,&e_cmax);  e_row_weight = malloc(sizeof(int)*e_m);  for (i = 0; i <= e_m-1; i++) fscanf(fp,"%d",&e_row_weight[i]);  e_col_weight = malloc(sizeof(int)*e_n);  for (i = 0; i <= e_n-1; i++) fscanf(fp,"%d",&e_col_weight[i]);  e_row_list = malloc(sizeof(int*)*e_m);  for (i=0;i<=e_m-1;i++) e_row_list[i] = malloc(sizeof(int)*e_rmax);  for (j = 0; j <= e_m-1; j++) {    for (i = 0; i <= e_row_weight[j]-1; i++) {      fscanf(fp,"%d",&v);      v--;      e_row_list[j][i] = v;    }  }}void print_encfile(){  int i,j;  printf("encoder file ----------------------\n");  printf("%d %d\n",e_n,e_m);  printf("%d %d\n",e_rmax,e_cmax);  for (i=0;i<=e_m-1;i++) {    for (j=0;j<=e_row_weight[i]-1;j++)  printf("%d ",e_row_list[i][j]+1);    printf("\n");  }}void gen_info_seq(int *seq){  int i;  for (i = 0; i <= m-1; i++) seq[i] = 0;  for (i = m; i <= n-1; i++) {    if (drand48() > 0.5) seq[i] = 1;    else seq[i] = 0;  }  return;}void print_seq(int *seq){  int i;  for (i = 0; i <= n-1; i++) printf("%d ",seq[i]);  printf("\n");}void encode_word(int* word){  int i,j;  int parity;  for (j = m-1; j >= 0; j--) {    parity = 0;    for (i = 1; i <= e_row_weight[j]-1; i++)       parity += word[e_row_list[j][i]];    word[j] = parity % 2;  }}int parity_check_func(int *seq){  int i,j;  int flag,parity;  flag = 0;  for (i=0; i <= m-1; i++) {    parity = 0;    for (j=0; j <= row_weight[i]-1; j++)       parity = (parity + seq[row_list[i][j]]) % 2;    if (parity == 1) flag = -1;  }  return flag;}void print_spmat(){  int i,j;  printf("----------------------\n");  printf("%d %d\n",n,m);  printf("%d %d\n",rmax,cmax);  for (i=0;i<=m-1;i++) {    for (j=0;j<=row_weight[i]-1;j++) printf("%d ",row_list[i][j]);    printf("\n");  }  for (i=0;i<=n-1;i++) printf("%d",punc_ptn[i]);  printf("\n");}/**********************************************************************************/// functions for quantized BP algorithm/**********************************************************************************/double nrnd(double var)				// gaussian random number generator{  static int sw = 0;  static double r1, r2, s;  if (sw == 0) {    sw = 1;    do {      r1 = 2 * drand48() - 1;      r2 = 2 * drand48() - 1;      s = r1 * r1 + r2 * r2;    } while (s > 1 || s == 0);    s = sqrt(-2 * log(s) / s);    return r1 * s * sqrt(var);  } else {    sw = 0;    return r2 * s * sqrt(var);  }}int quantize(double v){  int idx;  if (v >= DELTA/2.0) idx = (int)floor(v/DELTA + 0.5);  if (v <= -DELTA/2.0) idx = (int)ceil(v/DELTA - 0.5);  if ((v < DELTA/2.0) && (v > -DELTA/2.0)) idx = 0;  if (v > UB2*DELTA) idx = UB2;  if (v < LB2*DELTA) idx = LB2;  return idx;}int quantized_BP(double *llr) {  int i,j,k;  int sum;  int loop;  int posterior;  int a,b;  int flag,parity,ret;  int beta_tmp;  int rwm1;				// initialization of alpha  for (i=0; i <= m-1; i++) {    for (j=0; j <= row_weight[i]-1; j++) {      alpha[i][j] = 0;    }  }  for (i = 0; i <= n-1; i++) qLLR[i] = quantize(llr[i]);  ret = -1;  for (loop = 0; loop <= max_itr-1; loop++) { // start of iteration     				// variable node operation    for (i = 0; i <= n-1; i++) {      sum = 0;      for (k = 0; k <= col_weight[i]-1; k++) 	sum += alpha[col_list_r[i][k]][col_list_c[i][k]];      for (j = 0; j <= col_weight[i]-1; j++) {		if (punc_ptn[i] == 0) 	  beta_tmp = qLLR[i] + sum-alpha[col_list_r[i][j]][col_list_c[i][j]];	else 	  beta_tmp = sum-alpha[col_list_r[i][j]][col_list_c[i][j]];	if ((minsum == ON) && (norm_factor > 0.0)) 	  beta_tmp = quantize((double)beta_tmp * DELTA/norm_factor);	if (beta_tmp < LB2) beta_tmp = LB2;	if (beta_tmp > UB2) beta_tmp = UB2;	beta[col_list_r[i][j]][col_list_c[i][j]] = beta_tmp;      }			      if (punc_ptn[i] == 0) 	posterior = qLLR[i] + sum;      else 	posterior = sum;#ifdef TEST1      printf("%f ",posterior * DELTA);#endif      if (posterior > 0) tmp_bit[i] = 0;      if (posterior < 0) tmp_bit[i] = 1;      if (posterior == 0) {	if (drand48() < 0.5) tmp_bit[i] = 0;	else tmp_bit[i] = 1;      }    }				#ifdef TEST1    printf("\n");#endif 				// check node operation (using table lookup)    for (i = 0; i <= m-1; i++) {       rwm1 = row_weight[i]-1;				// foward computation      for (k = 0; k <= rwm1; k++) { 	if (k == 0) {	  a = beta[i][k];	}	else {	  b = beta[i][k];	  a = R_tbl[a+UB][b+UB];	}	forward[k] = a;      }				// backward computation      for (k = rwm1; k >= 0 ; k--) { 	if (k == rwm1) {	  a = beta[i][k];	}	else {	  b = beta[i][k];	  a = R_tbl[a+UB][b+UB];	}	backward[k] = a;      }				// update of alpha       for (k = 0; k <= rwm1; k++) { 	if (k == 0) alpha[i][k] = backward[1];	else if (k == rwm1) alpha[i][k] = forward[rwm1-1];	else alpha[i][k] = R_tbl[forward[k-1]+UB][backward[k+1]+UB];      }    }							/* parity check */    if (parity_check == ON) {      flag = 0;      for (i=0; i <= m-1; i++) {	parity = 0;	for (j=0; j <= row_weight[i]-1; j++) 	  parity = (parity + tmp_bit[row_list[i][j]]);	if (parity % 2 == 1) flag = 1;      }      if (flag == 0) {	ret = 0;	break;      }    }  }