awgn_simulator.c

此程序是c语言的有关LDPC译码仿真程序

/*
**********************************************************************************
  	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.4987
quantize_depth	10
max_iteration	100
random_seed	1
display		1
stop_condition 	1
stop_errors	100
target_prob	1e-5
error_logging	0
puncture	0
multi_edge_type	0
min-sum		0
norm_factor	0.0
parity_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_constraint
4 3                               # largest_row_weight largest_column_weight
4 4 4 4 4 4 4 4 4 4 4 4           # row_weight_profile
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3   # column_weight_profile
3 8 10 13                         # parity check matrix(caution:colum # started from 1!)
4 7 9 13
2 5 7 10
4 6 11 14
3 9 15 16
1 6 9 10
4 8 12 15
2 6 12 16
1 7 14 16
3 5 12 14
2 11 13 15
1 5 8 11                    

The above spmat format corresponds to the following matrix:

16
4
 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 0

If puncture = 1, append puncture pattern like
 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0
to 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 bit
int **alpha;		// check to variable message
int **beta;		// variable to check message
int n;			// code length
int m;			// number of constraints
int rmax;		// maximum number of row weight
int cmax;		// maximum number of column weight
int *row_weight;	// row weights
int *col_weight;	// column weights
int **row_list;		// row list of parity check matrix
int **col_list_r;	// column list 1
int **col_list_c;	// column list 1
int LB,UB,LB2,UB2;	// upper and lower bound for quantized messages
int ASIZE;		// number of quantization levels
double DELTA;		// quantization step size
int **R_tbl;		// table for check node operation
double *LLR;		// array for received LLR values
int *qLLR;		// quantized LLR values
int *forward;		// forward array for check node operation
int *backward;		// backward array for check node operation
int max_itr;		// maximum number of iterations in BP

//-------------------------------------for simulation_loop()
double code_rate;	// code rate
int total_iterations;	// total number of iterations
int seed;		// seed for random number generator
int quantization_level;	// quantization depth
char *filename;		// filename of configuration file
int total_blocks;	// total number of tested blocks
int total_bits;		// total number of tested bits
int error_blocks;	// number of error blocks
int error_bits;		// number of error bits
int total_miscorrection; // number of mis corrected blocks
int display_result;	 // display switch
int stop_condition;	 // stop condition
int stop_errors;	 // number of errors for stop condition
int error_logging;	 // error logging switch
double target_error_prob; // target error probability
int *punc_ptn;		 // puncture pattern
int puncture;		 // puncture switch
int met;		 // multi-edge type switch
int *vnode_type;	 // variable node type
int num_vnode_type;	 // number of variable node type
FILE *error_log;	 // error log file
int *vnode_errors;	 // number of variable node errors
int minsum;		 // min-sum algorithm switch
double norm_factor;	 // normalization factor
int parity_check;	 // parity check switch
int encoding;

//-------------------------------------for encoding
int e_n;                // code length
int e_m;		// number of parity symbols
int e_rmax;		// maximum number of row weight
int e_cmax;		// maximum number of column weight
int *e_row_weight;	// row weights
int *e_col_weight;	// column weights
int **e_row_list;	// parity check matrix

//-------------------------------------channel dependent parameters
int *cword;             // sent codeword
double *rword;		// received word at receiver
double channel_var;     // channel variance
double 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;
      }
    }
  }