awgn_simulator.c

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

  total_iterations += loop;
  return ret;
}

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) return 0.0;
  if (a > 0) return 1.0;
  if (a < 0) return -1.0;
}

double min(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(fabs(a),fabs(b));
}

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

  int i,j;
  int v;
  int *counter;
  double a,b;
  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(int*)*m); // check to variable message
  for (i=0;i<=m-1;i++) alpha[i] = malloc(sizeof(int)*rmax);
  beta  = malloc(sizeof(int*)*m); // variable to check message
  for (i=0;i<=m-1;i++) beta[i] = malloc(sizeof(int)*rmax);
  tmp_bit = malloc(sizeof(int)*n); // tentative decoding result

				// set up of table for check node operation
  R_tbl = malloc(sizeof(int*) * ASIZE);
  for (i = 0; i <= ASIZE-1; i++) {
    R_tbl[i] = malloc(sizeof(int) * ASIZE);
  }
  for (i = LB2; i <= UB2; i++) {
    for (j = LB2; j <= UB2; j++) {
      a = DELTA * i;
      b = DELTA * j;
      if (minsum == 1) 
	R_tbl[i+UB][j+UB] = quantize(minsumR(a,b));
      else
	R_tbl[i+UB][j+UB] = quantize(R(a,b));
    }
  }
  qLLR = malloc(sizeof(int)*n);	// quantized LLR
  forward = malloc(sizeof(int)*rmax);
  backward= malloc(sizeof(int)*rmax);
}

/**********************************************************************************/
// simulation status and result report functions
/**********************************************************************************/

void report_result2(FILE *out)
{
  fprintf(out,"%4.3f %7.6e %7.6e %d %d %d %d %5.2f %d %d\n",
	  snr,
	  (double)error_bits/total_bits,
	  (double)error_blocks/total_blocks,
	   error_bits,
	  total_bits,
	  error_blocks,
	  total_blocks,
	  (double)total_iterations/total_blocks,
	  n,
	  m
	  );
}

void report_result(FILE *out)
{
  fprintf(out,"%4.3f %7.6e %7.6e %d %d %d %d %5.2f var=%4.3f sdv=%4.3f mis=%d seed=%d mitr=%d n=%d m=%d r=%4.3f depth=%d DLETA=%4.3e range=%4.3e config=%s scond=%d serr=%d target=%4.3e log=%d punc=%d met=%d msum=%d norm=%4.3f pchk=%d\n",
	  snr,
	  (double)error_bits/total_bits,
	  (double)error_blocks/total_blocks,
	  error_bits,
	  total_bits,
	  error_blocks,
	  total_blocks,
	  (double)total_iterations/total_blocks,
	  channel_var,
	  sqrt(channel_var),
	  total_miscorrection,
	  seed,
	  max_itr,
	  n,
	  m,
	  code_rate,
	  quantization_level,
	  DELTA,
	  UB2*DELTA,
	  filename,
	  stop_condition,
	  stop_errors,
	  target_error_prob,
	  error_logging,
	  puncture,
	  met,
	  minsum,
	  norm_factor,
	  parity_check
	  );
}

/**********************************************************************************/
// simulation functions
/**********************************************************************************/

void make_received_word(double *r)
{
  int i;
  for (i = 0; i <= n-1; i++) {
    r[i] = 1.0-2.0*cword[i] + nrnd(channel_var);
  }
}

void simulation_loop()
{

  int i,j;
  int sum;
  int ret;

// set up of initial condition

  cword = malloc(sizeof(int)*n);
  rword = malloc(sizeof(double)*n);
  LLR = malloc(sizeof(double)*n);

  total_blocks = 0;
  total_bits = 0;
  error_blocks = 0;
  error_bits = 0;
  total_iterations = 0;
  total_miscorrection = 0;

  while(1) {
    total_blocks++;
    total_bits += n;

    if (encoding == 1) {
      gen_info_seq(cword);
      encode_word(cword);
    } 
    else for (i = 0; i <= n-1;i++) cword[i] = 0;
    //for (i =0; i <= n-1; i++) printf("%d",cword[i]);printf("\n");
    //printf("parity = %d\n",parity_check_func(cword));
    make_received_word(rword);
    for (i = 0; i <= n-1; i++) LLR[i] = 2.0 * rword[i]/channel_var;
    ret = quantized_BP(LLR);
    sum = 0;
    for(i = 0; i <= n-1;i++) {
      if (encoding == 1) tmp_bit[i] ^= cword[i];
      if (punc_ptn[i] == 0) sum += tmp_bit[i];
      if (met == ON) 
	vnode_errors[vnode_type[i]] += tmp_bit[i];
    }
    error_bits += sum;
    if (sum != 0) {
      error_blocks++;
      if (error_logging==1) {
	for (j=0; j<=n-1;j++) fprintf(error_log,"%d ",tmp_bit[j]);
	fprintf(error_log,"(error_weight=%d)\n",sum);
      }
    }
    if ((ret == 0)&&(sum != 0)) total_miscorrection++;
    if (display_result == 1) report_result(stderr); // full report
    if (display_result == 2) report_result2(stderr); // brief report

    if ((error_bits >= stop_errors)&&(stop_condition==BIT)&&(total_blocks>=100)) break;
    if ((error_blocks >= stop_errors)&&(stop_condition==BLOCK)) break;
    if ((total_bits >= stop_errors/target_error_prob)&&(stop_condition==BIT)) {
      fprintf(stdout, "stop by target_error_prob conditon\n");
      break;
    }
    if ((total_blocks >= stop_errors/target_error_prob)&&(stop_condition==BLOCK)) {
      fprintf(stdout, "stop by target_error_prob conditon\n");
      break;
    }
    if ((total_bits >= 10/target_error_prob)&&(error_bits==0)&&(stop_condition==BIT)) {
      fprintf(stdout, "stop by early stop condition\n");
      break;
    }
    if ((total_blocks >= 10/target_error_prob)&&(error_blocks==0)&&(stop_condition==BLOCK)) {
      fprintf(stdout, "stop by early stop condition\n");
      break;
    }
    if ((stop_condition >=2)&&(total_blocks >stop_condition)) break;
  }
}

void terminate_program(char *s)
{
    fprintf(stderr,"%s\n",s);
    exit(-1);
}

void read_config_file(FILE *fp)
{
// set up for simulation parameters

  FILE *hfile,*encoder;
  char str1[MAX_STR_LEN],str2[MAX_STR_LEN];
  int i;

  fscanf(fp,"%s",str1);
  if (strcmp(str1,"matrix_file")!=0) terminate_program("error in matrix_file\n");
  fscanf(fp,"%s",str2);		// parity check matrix file (spmat form)
  if ((hfile = fopen(str2,"r")) == NULL) {
    fprintf(stderr,"Can't open %s.\n",str2);
    exit(-1);
  }
  fscanf(fp,"%s",str1);
  if (strcmp(str1,"code_rate")!=0) terminate_program("error in code_rate\n");
  fscanf(fp,"%lf",&code_rate);	// code rate of the target LDPC code
  if ((code_rate<0) || (code_rate>1)) terminate_program("invalid code rate\n");

  fscanf(fp,"%s",str1);
  if (strcmp(str1,"quantize_depth")!=0) terminate_program("error in quantize_depth\n");
  fscanf(fp,"%d",&quantization_level); // quantization depth
  if ((quantization_level<2) || (quantization_level>20)) 
    terminate_program("invalid quantization_level\n");
  ASIZE = (1 << quantization_level)-1; // number of quantization levels
  DELTA = 102.4/(ASIZE+1);	// quatization step

  LB = -(ASIZE-1)/2;		
  UB = (ASIZE-1)/2;
  LB2 = -(ASIZE-1)/4;		// lower bound of quantized value
  UB2 = (ASIZE-1)/4;		// upper bound of quantized value

  fscanf(fp,"%s",str1);
  if (strcmp(str1,"max_iteration")!=0) terminate_program("error in max_iteration\n");
  fscanf(fp,"%d",&max_itr);	// maximum number of iterations in BP
  if (max_itr < 0) terminate_program("invalid max_itr\n");

  fscanf(fp,"%s",str1);
  if (strcmp(str1,"random_seed")!=0) terminate_program("error in random_seed\n");
  fscanf(fp,"%d",&seed);	// seed of random number generator
  srand48(seed);		// set up random seed

  fscanf(fp,"%s",str1);
  if (strcmp(str1,"display")!=0) terminate_program("error in display\n");
  fscanf(fp,"%d",&display_result); // flag for display mode of simulation status
				   // 0:no display, 1:full display, 2:simple display
  fscanf(fp,"%s",str1);
  if (strcmp(str1,"stop_condition")!=0) terminate_program("error in stop_condition\n");
  fscanf(fp,"%d",&stop_condition); // stop condition of simulation
				   // 0:bit errors, 1:block errors, 
				   // other positive number: number of loops
  if (stop_condition < 0) terminate_program("invalid stop_condition\n");

  fscanf(fp,"%s",str1);
  if (strcmp(str1,"stop_errors")!=0) terminate_program("error in stop_errors\n");
  fscanf(fp,"%d",&stop_errors);	// When the number of errors reaches stop_errors,
				// simulation loop is terminated.
  if (stop_errors < 0) 
    terminate_program("invalid stop_errors\n");

  fscanf(fp,"%s",str1);
  if (strcmp(str1,"target_prob")!=0) terminate_program("error in target_prob\n");
  fscanf(fp,"%lf",&target_error_prob); // target error probability
  if (target_error_prob < 0.0) 
    terminate_program("invalid target_error_prob\n");

  fscanf(fp,"%s",str1);
  if (strcmp(str1,"error_logging")!=0) terminate_program("error in error_loging\n");
  fscanf(fp,"%d",&error_logging); // error logging switch 0 or 1
  if ((error_logging!=0)&&(error_logging!=1)) 
    terminate_program("invalid error_logging\n");
  if(error_logging == 1) {
    if ((error_log = fopen("error_log","w")) == NULL) {
      fprintf(stderr,"Can't open %s.\n","error_log");
      exit(-1);
    }
    fprintf(stderr,"error logging starts.\n");
  }

  fscanf(fp,"%s",str1);
  if (strcmp(str1,"puncture")!=0) terminate_program("error in puncture\n");
  fscanf(fp,"%d",&puncture);	// Are there puncture bits? 0 or 1
  if ((puncture!=0)&&(puncture!=1)) 
    terminate_program("invalid puncture\n");

  fscanf(fp,"%s",str1);
  if (strcmp(str1,"multi_edge_type")!=0) terminate_program("error in multi_edge_type\n");
  fscanf(fp,"%d",&met);		// Is it a multi-edge type LDPC code? 0 or 1
  if ((met!=0)&&(met!=1)) 
    terminate_program("invalid met\n");

  fscanf(fp,"%s",str1);
  if (strcmp(str1,"min-sum")!=0) terminate_program("error in min-sum\n");
  fscanf(fp,"%d",&minsum);	// switch of min-sum algorithm 0 or 1
				// if minsum == 1, then simulation with min-sum algorithm
				// otherwise simulation with sum-product algortihm
  if ((minsum!=0)&&(minsum!=1)) 
    terminate_program("invalid minsum\n");

  fscanf(fp,"%s",str1);
  if (strcmp(str1,"norm_factor")!=0) terminate_program("error in norm_factor\n");
  fscanf(fp,"%lf",&norm_factor); // normalization factor for min-sum algorithm
				 // if normalization == 0, then no normalization
				 // See papers of normalized-BP by Prof.Fossorier.
  if (norm_factor < 0.0) 
    terminate_program("invalid norm_factor\n");

  fscanf(fp,"%s",str1);
  if (strcmp(str1,"parity_check")!=0) terminate_program("error in parity_check\n");
  fscanf(fp,"%d",&parity_check); // switch of parity check in BP algortihm, 0 or 1
  if ((parity_check!=0)&&(parity_check!=1)) 
    terminate_program("invalid parity_check\n");

  #ifdef TEST2
  encoding = 1;
  printf("encoding is ON\n");
  fscanf(fp,"%s",str1);
  if ((encoder = fopen(str1,"r")) == NULL) {
    fprintf(stderr,"Can't open %s.\n",str1);
    exit(-1);
  }
  read_encfile(encoder);
  //print_encfile(encoder);
  #endif

  init_dec(hfile);		// initialization of decoder
}

/**********************************************************************************/
// main 
/**********************************************************************************/

int main(int argc,char **argv)
{
  FILE *fp;			// pointer for configuration file
  int i;

#ifdef TEST1
/*
Test for quantized_BP()
-----------------------------------------------------------------
6.3.spmat
------------------
6 3
3 2
3 2 3
1 1 2 2 1 1
1 2 3 
3 4
4 5 6
1 0 0 0 0 1
------------------
log posterior probability ratio (by double-precision BP)

    4.0554    0.5884    0.2407    1.0694    1.6060    2.2915
    3.9537    0.2780    1.7111    1.7111    1.5387    2.2338
    4.3974    1.6925    1.7111    1.7111    2.0840    2.7033
    4.3974    1.6925    1.7111    1.7111    2.0840    2.7033
    4.3974    1.6925    1.7111    1.7111    2.0840    2.7033

*/
  double r[] = {1.620803,0.264281,-0.031637,-0.127654,0.746347,1.003543}; 
  int qlevel;

  channel_var = 0.794328;
  printf("basic test\n");
  fp = fopen("6.3.spmat","r");
  qlevel = 10;
  ASIZE = (1 << qlevel)-1; 
  DELTA = 102.4/(ASIZE+1);
  LB = -(ASIZE-1)/2;		
  UB = (ASIZE-1)/2;
  LB2 = -(ASIZE-1)/4;		
  UB2 = (ASIZE-1)/4;		
  puncture = OFF;
  met = OFF;
  parity_check = OFF;
  max_itr = 10;

  init_dec(fp);
  LLR = malloc(sizeof(int)*6);
  for (i = 0;i <= 5; i++) LLR[i] = 2.0*r[i]/channel_var;
  quantized_BP(LLR);
  exit(0);
//-----------------------------------------------------------------
#endif

// user interface

  if (argc != 3) {
    printf("usage: awgn_simulator config_file snr\n");
				// snr: Eb/N0
    exit(-1);
  }
  if ((fp = fopen(argv[1],"r")) == NULL) {
    fprintf(stderr,"Can't open %s.\n",argv[1]);
    exit(-1);
  }
  filename =argv[1];		// filename of configuration file
  snr = atof(argv[2]);		// get signal to noise ratio
  read_config_file(fp);		// reading configulation file
  channel_var = 0.5 * (1.0/pow(10.0,snr/10.0))/code_rate;

// simulation and report

  simulation_loop();		// main simulation loop
  report_result(stdout);	// print out the simulation result

  if (error_logging==ON) report_result(error_log);
				// for multi-edge type
  if (met == ON) {
    printf("# number of bit errors for each variable node type\n");
    printf("# ");
    for (i=0;i<=num_vnode_type-1;i++) printf("%d ",vnode_errors[i]);
    printf("\n");
    printf("# bit error probability for each variable node type\n");
    printf("# ");
    for (i=0;i<=num_vnode_type-1;i++) printf("%16.12e ",(double)vnode_errors[i]/total_bits);
    printf("\n");
  }
  return 0;
}