ldpc_awgn.c

用C编写的通信过程中编译码,主要用于纠错.望大家指教

/*
**********************************************************************************
  	FILENAME: ldpc_awgn.c
	version 1.0
	(C) jftian, 2003

        An LDPC code simulator based on spa algorithm for AWGN chanel
*********************************************************************************
        ChangeLog:
		1st,Jan., 2004 : Build

**********************************************************************************
  How to run
*********************************************************************************

    ldpc_awgn filename_of_configuration_file(.cfg) snr
    Example: ldpc_awgn H_15_7 1.0

/**********************************************************************************/
// Inlucde files and macros
/**********************************************************************************/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <time.h>

#define OFF 0
#define ON 1
#define BIT 0
#define BLOCK 1
#define MAX_STR_LEN 100
#define ran01 (rand()/((double)RAND_MAX+0.1))
//#define TEST1		// simulate for spa()
//#define TEST_tian	// simulate with no input
#define TEST2		// simulate with encoding
//#define debug

/**********************************************************************************/
// Global variables
/**********************************************************************************/

//-------------------------------------for spa()
int *tmp_bit;		// tentetative decision bit
double **alpha;		// check to variable message
double **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
double *LLR;		// array for received LLR values
double *forward;		// forward array for check node operation
double *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
char *file_name;		// filename of configuration file(no suffix)
char *char_snr;
char char_snr_add[MAX_STR_LEN];
char config_file[MAX_STR_LEN], log_file[MAX_STR_LEN];
char *suffix1;
char *suffix2;
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 usage(void)
{ fprintf (stderr, 
   "Usage:   ldpc_awgn config-file[*.cfg] SNR\n");
  exit(1);
}

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 <= e_m-1; i++) seq[i] = 0;
  for (i = e_m; i <= n-1; i++) {
    if (ran01> 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 = 0; j < e_m; j++) {
    parity = 0;
    for (i = 0; i < e_row_weight[j]; i++) 
      parity += word[e_m+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;break;}
  }
  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 * ran01 - 1;
      r2 = 2 * ran01 - 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);
  }
}



double tanh(double x)
{
	return (exp(2*x)-1)/(exp(2*x)+1);
}

double atanh(double x)
{
	return log((1+x)/(1-x))/2;
}

double fi_ldpc(double x)
{
	return -log(tanh(x/2.0));
}

double R(double a, double b)
				// function for check node operation
{
  return 2.0 * atanh(tanh(a/2.0)*tanh(b/2.0));
}

double sign(double a)
{
  if (a == 0.0) return 0.0;
  else if (a > 0.0) return 1.0;
  else return -1.0;
}


int float_spa(double *llr) 
{
  int i,j,k;
  double sum;
  int loop;
  double posterior;
  double a,b;
  int flag,parity,ret;
  double beta_tmp;
  int rwm1;
  double coef;
				// 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 = llr[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 = llr[i] + sum;
      else 
	posterior = sum;
#ifdef TEST1
      printf("%f ",posterior);
#endif
      if (posterior > 0) tmp_bit[i] = 0;
      if (posterior < 0) tmp_bit[i] = 1;
      if (posterior == 0) {
	if (ran01 < 0.5) tmp_bit[i] = 0;
	else tmp_bit[i] = 1;
      }
    }				
#ifdef TEST1                 
    printf("\n");
#endif 
				// check node operation (use spa )
    for (i = 0; i <= m-1; i++) { 
      rwm1 = row_weight[i]-1;
				// foward computation
      for (k = 0; k <= rwm1; k++) { 
	if (k == 0) {
	  a = fi_ldpc (fabs(beta[i][k]));
	  coef = sign(beta[i][k]);
	}
	else {
	  b = fi_ldpc(fabs(beta[i][k]));
	  a = a+b;     //R_tbl[a+UB][b+UB];
	  coef = coef*sign(beta[i][k]);
	}
	#ifdef debug
	  printf("%f %f %f\n",a,b,coef);
	  
#endif
	forward[k] = a;
      }
				// backward computation
      for (k = rwm1; k >= 0 ; k--) { 
	if (k == rwm1) {
	  a = fi_ldpc(fabs(beta[i][k]));
	}
	else {
	  b = fi_ldpc(fabs(beta[i][k]));
	  a = a+b;     //R_tbl[a+UB][b+UB];
	}
	backward[k] = a;
      }
				// update of alpha 
      for (k = 0; k <= rwm1; k++) { 
	if (k == 0) alpha[i][k] = fi_ldpc(backward[1])*coef/sign(beta[i][k]);
	else if (k == rwm1) alpha[i][k] = fi_ldpc(forward[rwm1-1])*coef/sign(beta[i][k]);
	else alpha[i][k] = fi_ldpc(forward[k-1]+backward[k+1])*coef/sign(beta[i][k]);  

#ifdef debug
	printf("%f\n",alpha[i][k]);
#endif
	           //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;
      }
    }
  }	
  total_iterations += loop;
  return ret;
}


double min_d(double a, double b)
{
  if (a > b) return b;
  else return a;
}

double minsumR(double a, double b)
				// function for check node operation of min-sum algotithm
{
  return sign(a) * sign(b) * min_d(fabs(a),fabs(b));
}

void init_dec(FILE *fp)
{
/**********************************************************************************/
// initialization of spa decoder
/**********************************************************************************/

  int i,j;
  int v;
  int *counter;
  int tmp;

// read parity check matrix 

  fscanf(fp,"%d %d\n",&n,&m);	// n:code length, m:number of parity checks
  fscanf(fp,"%d %d\n",&rmax,&cmax); // rmax:maximum number of row weight
				    // cmax:maximum number of column weight
  row_weight = malloc(sizeof(int)*m); // array for row weight profile
  for (i = 0; i <= m-1; i++) fscanf(fp,"%d",&row_weight[i]);
  col_weight = malloc(sizeof(int)*n); // array for column weight profile
  for (i = 0; i <= n-1; i++) fscanf(fp,"%d",&col_weight[i]);
  counter = malloc(sizeof(int)*n); 
  for (i=0; i <= n-1; i++) counter[i] = 0;  

  row_list = malloc(sizeof(int*)*m); // parity check matrix in row form
  for (i=0;i<=m-1;i++) row_list[i] = malloc(sizeof(int)*rmax);
  col_list_r = malloc(sizeof(int*)*n); 
  for (i=0;i<=n-1;i++) col_list_r[i] = malloc(sizeof(int)*cmax);
  col_list_c = malloc(sizeof(int*)*n); 
  for (i=0;i<=n-1;i++) col_list_c[i] = malloc(sizeof(int)*cmax);

// set up for parity check matrix

  for (j = 0; j <= m-1; j++) {
    for (i = 0; i <= row_weight[j]-1; i++) {
      fscanf(fp,"%d",&v);
      v--;
      row_list[j][i] = v;
      col_list_r[v][counter[v]] = j;
      col_list_c[v][counter[v]] = i;
      counter[v]++;
    }
  }
  punc_ptn = malloc(sizeof(int) * n); // puncture pattern (0:no puncture, 1:puncture)
  if (puncture == ON)
    for (i = 0; i<=n-1; i++) fscanf(fp,"%d",&punc_ptn[i]);
  else 
    for (i = 0; i<=n-1; i++) punc_ptn[i] = 0;
  if (met == ON) {
    vnode_type = malloc(sizeof(int)*n);
    fscanf(fp,"%d",&tmp);
    for (i = 0; i<=m-1; i++) fscanf(fp,"%d",&tmp);
    fscanf(fp,"%d",&num_vnode_type);    
    for (i = 0; i<=n-1; i++) fscanf(fp,"%d",&vnode_type[i]);
    vnode_errors = malloc(sizeof(int)*num_vnode_type);
    for (i = 0; i <= num_vnode_type-1; i++) vnode_errors[i] = 0;
  }

  // set up for BP related global variables

  alpha = malloc(sizeof(double*)*m); // check to variable message
  for (i=0;i<=m-1;i++) alpha[i] = malloc(sizeof(double)*rmax);
  beta  = malloc(sizeof(double*)*m); // variable to check message
  for (i=0;i<=m-1;i++) beta[i] = malloc(sizeof(double)*rmax);
  tmp_bit = malloc(sizeof(int)*n); // tentative decoding result

				// set up of table for check node operation