ldpc_awgn.c

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

  forward = malloc(sizeof(double)*rmax);
  backward= malloc(sizeof(double)*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,"-------------------------------------\n");
	fprintf(out,"SNR=%4.3f BER=%7.6e FER=%7.6e err_bits=%d tot_bits=%d err_blk=%d tot_blk%d  mis=%d\n",
		         snr,
				 (double)error_bits/total_bits,
				 (double)error_blocks/total_blocks,
				 error_bits,//4
				 total_bits,//5
				 error_blocks,//6
				 total_blocks,
				 total_miscorrection);
	fprintf(out,"var=%4.3f sdv=%4.3f seed=%d max_itr=%d n=%d m=%d r=%4.3f config=%s \n",
			  channel_var,//8
	          sqrt(channel_var),
			  seed,//11
	          max_itr,//12
	          n,//13
	          m,//14
	          code_rate,//15
	          config_file);
	fprintf(out,"scond=%d serr=%d target=%4.3e log=%d punc=%d met=%d msum=%d norm=%4.3f pchk=%d\n",
		stop_condition,
	    stop_errors,
	    target_error_prob,
	    error_logging,
	    puncture,
	    met,
	    minsum,
	    norm_factor,
	    parity_check);
/*
	fprintf(out,"%4.3f %7.6e %7.6e %d %d %d %d var=%4.3f sdv=%4.3f mis=%d seed=%d mitr=%d n=%d m=%d r=%4.3f config=%s scond=%d serr=%d target=%4.3e log=%d punc=%d met=%d msum=%d norm=%4.3f pchk=%d\n",
	  snr,//1
	  (double)error_bits/total_bits,//2
	  (double)error_blocks/total_blocks,//3
	  error_bits,//4
	  total_bits,//5
	  error_blocks,//6
	  total_blocks,//7
//	  (double)total_iterations/total_blocks,
	  channel_var,//8
	  sqrt(channel_var),//9
	  total_miscorrection,//10
	  seed,//11
	  max_itr,//12
	  n,//13
	  m,//14
	  code_rate,//15
	  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;
  int ttt;

// 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;
//	print_col();

    if (encoding == 1) {
      gen_info_seq(cword);
//		cword[3]=1;
//		cword[4] = 0;
//		cword[5] = 1;
	  encode_word(cword);
//	  for (i=0; i<n; i++)
//	  {
//	  printf("%d ",cword[i]);
//	  }
//	  printf("\n");
//	  re_order(cword);
//	  for (i=0; i<6; i++)
//	  {
//	  printf("%d ",cword[i]);
//	  }
	  ttt = parity_check_func(cword);
	  if (ttt!=0)
	  {
		  printf("encoding wrong!!\n");
		  exit(-1);
	  }
//	  printf("parity::%d\n",ttt);
//	  if (ttt!=0) 
//	  {
//		  printf("parity is wrong!\n");
//			  exit(1);
//	  }
    } 
    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);  // revised by tian 12-21
	ret = float_spa(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,"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
  if(seed>0)  srand(seed);		// set up random seed
  else 
  srand( (unsigned) time(NULL) );


  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(log_file,"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 (strcmp(str1,"generate_matrix")!=0) terminate_program("error in gen_matrix\n");
  fscanf(fp,"%s",str2);
  if ((encoder = fopen(str2,"r")) == NULL) {
    fprintf(stderr,"Can't open %s.\n",str1);
    exit(-1);
  }

/*
  fscanf(fp,"%s",str1);
  if (strcmp(str1,"order_file")!=0) terminate_program("error in order_file\n");
  fscanf(fp,"%s",str2);
  if ((file_order = fopen(str2,"r")) == NULL) {
    fprintf(stderr,"Can't open %s.\n",str2);
    exit(-1);
  }
*/  

  read_encfile(encoder);
//  read_order_file(file_order);
  
  //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;
  suffix1=".cfg";
  suffix2=".log";

#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}; 
  double r[] = {0.2,0.2,-0.9,0.6,0.5,-1.1,-0.4,-1.2};

  channel_var = 0.5; //0.794328;
  printf("basic test\n");
  fp = fopen("12.21.spmat","r");
  puncture = OFF;
  met = OFF;
  parity_check = OFF;
  max_itr = 10;

  init_dec(fp);
  LLR = malloc(sizeof(int)*n);
  for (i = 0;i < n; i++) LLR[i] = 2.0*r[i]/channel_var;
  #ifdef debug
	for (i=0;i<n;i++)
		printf("%f  ",LLR[i]);
	printf("\n");
  #endif
  float_spa(LLR);
  exit(0);
//-----------------------------------------------------------------
#endif

// user interface

// add by tjf 2003-12-21
  #ifdef TEST_tian
	printf("use the parameter defined by tian \n");
	file_name = "H_255_175";
	snr = 3.0;
  #else
  if (argc != 3) {
    usage();
  }

  if (!(file_name = argv[1]) )
  {  	 usage();
  }


  
  snr = atof(argv[2]);		// get signal to noise ratio
  char_snr = argv[2];

  
  #endif
  
  strcpy(config_file,file_name);
  strcat(config_file,suffix1);

  strcpy(log_file,file_name);
  strcpy(char_snr_add,char_snr);
  strcat(log_file,"_");
  strcat(log_file,char_snr_add);
  strcat(log_file,"dB");
  strcat(log_file,suffix2);
  
  if ((fp = fopen(config_file,"r")) == NULL) {
    fprintf(stderr,"Can't open %s.\n",config_file);
    exit(-1);
  }
  read_config_file(fp);		// reading configulation file
  code_rate = (double)e_n/n;
  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;
}