new_rs_erasures.c

一些纠错编码的c程序实现

	    if (discr_r == 0){
	       /* 3 lines below: B(x) <-- x*B(x) */
	       tmp_pol[0] = 0;
	       for (i=1;i < 2*tt+1;i++) tmp_pol[i] = b[i-1];
	       for (i=0;i < 2*tt+1;i++) b[i] = tmp_pol[i];
	    }
	    else{
	       /* 5 lines below: T(x) <-- lambda(x) - discr_r*x*b(x) */
	       T[0] = lambda[0];
	       for (i=1;i < 2*tt+1;i++){
		    tmp =  (b[i-1] == 0)? 0 : alpha_to[(index_of[discr_r]+index_of[b[i-1]])%nn];
		    T[i] = lambda[i]^tmp;
	       }

	       if (2*L <= r+no_eras-1){
		    L = r+no_eras-L;
		    /* 2 lines below: B(x) <-- inv(discr_r) * lambda(x) */
		    for (i=0;i < 2*tt+1;i++)
		        b[i] = (lambda[i] == 0)? 0 : alpha_to[(index_of[lambda[i]]-index_of[discr_r]+nn)%nn];
		    for (i=0;i < 2*tt+1;i++) lambda[i] = T[i]; 
	       }
	       else{
		    for (i=0;i < 2*tt+1;i++) lambda[i] = T[i]; 
	            /* 3 lines below: B(x) <-- x*B(x) */
	            tmp_pol[0] = 0;
	            for (i=1;i < 2*tt+1;i++) tmp_pol[i] = b[i-1];
	            for (i=0;i < 2*tt+1;i++) b[i] = tmp_pol[i];
	       }
	    }
          }
/* Put lambda(x) into index form */
    for (i=0; i < 2*tt+1; i++)
	lambda[i] = index_of[lambda[i]];
    
/* Compute deg(lambda(x)) */
    deg_lambda = 2*tt;
    while ((lambda[deg_lambda] == -1) && (deg_lambda > 0)) 
	--deg_lambda;
    if (deg_lambda <= 2*tt){
/* Find roots of the error+erasure locator polynomial. By Chien Search */
         for (i=1; i < 2*tt+1; i++) reg[i] = lambda[i] ;
         count = 0 ; /* Number of roots of lambda(x) */
         for (i=1; i <= nn; i++)
          {  q = 1 ;
             for (j=1; j <= deg_lambda; j++)
              if (reg[j] != -1)
                { reg[j] = (reg[j]+j)%nn ;
                  q ^= alpha_to[(reg[j])%nn] ;
                } ;
             if (!q)        /* store root (index-form) and error location number */
              { root[count] = i;
                loc[count] = nn-i;
                count++;
              };
          } ;
#ifndef NO_PRINT
	  printf("\n Final error positions:\t");
	  for (i=0;i < count;i++) printf("%d ",loc[i]);
	  printf("\n");
#endif 
	if (deg_lambda == count){ /* correctable error */
/* Compute err+eras evaluator poly omega(x) = s(x)*lambda(x) (modulo x**(nn-kk)). in poly-form */
	     for (i=0;i < 2*tt;i++){
		omega[i] = 0;
		for (j=0;(j < deg_lambda+1) && (j < i+1);j++){
		    if ((s[i+1-j] != -1) && (lambda[j] != -1))
		       	tmp = alpha_to[(s[i+1-j]+lambda[j])%nn];
		    else
			tmp = 0;
		    omega[i] ^= tmp;	
		}
	     }
	     omega[2*tt] = 0;
/* Compute lambda_pr(x) = formal derivative of lambda(x) in poly-form */
	for (i=0;i < tt;i++){
	   lambda_pr[2*i+1] = 0;
	   lambda_pr[2*i] = (lambda[2*i+1] == -1)? 0 : alpha_to[lambda[2*i+1]];
	}
	lambda_pr[2*tt] = 0;
/* Compute deg(omega(x)) */
    deg_omega = 2*tt;
    while ((omega[deg_omega] == 0) && (deg_omega > 0)) 
	--deg_omega;
/* Compute error values in poly-form. num1 = omega(inv(X(l))), 
  num2 = inv(X(l))**(b0-1) and den = lambda_pr(inv(X(l))) all in poly-form */
        for (j=0;j < count;j++){
	   pres_root = root[j];
	   pres_loc = loc[j];
	   num1 = 0;
	   for (i=0;i < deg_omega+1;i++){
	      if (omega[i] != 0) 
		tmp = alpha_to[(index_of[omega[i]]+i*pres_root)%nn];
	      else
		tmp = 0;
	      num1 ^= tmp;
	   }
	   num2 = alpha_to[(pres_root*(b0-1))%nn];
	   den = 0;
	   for (i=0;i < deg_lambda+1;i++){
	      if (lambda_pr[i] != 0)
		tmp = alpha_to[(index_of[lambda_pr[i]]+i*pres_root)%nn];
	      else
		tmp = 0;
	      den ^= tmp;
	   }
	   if (den == 0){
		printf("\n ERROR: denominator = 0\n");
	   }
	   err[pres_loc] = 0;
	   if (num1 != 0){
		err[pres_loc] = alpha_to[(index_of[num1]+index_of[num2]+(nn-index_of[den]))%nn];
	   }
	 }
/* Correct word by subtracting out error bytes. First convert recd[] into poly-form */
		for (i=0;i < nn;i++) recd[i] = (recd[i] == -1)? 0 : alpha_to[recd[i]];
		for (j=0;j < count;j++)
		    recd[loc[j]] ^= err[loc[j]];
		return(1);
	}
	else /* deg(lambda) unequal to number of roots => uncorrectable error detected */
	     return(2);
     }
     else /* deg(lambda) > 2*tt => uncorrectable error detected */
	 return(3);
   }
   else{  /* no non-zero syndromes => no errors: output received codeword */
       for (i=0; i < nn; i++){
          if (recd[i] != -1)        /* convert recd[] to polynomial form */
           recd[i] = alpha_to[(recd[i])%nn];
          else
           recd[i] = 0;
       }
       return(0);
   }
}


/* 
    Generates the Galois field and then the generator polynomial.
  Must be recompiled for each different generator polynomial after 
  setting parameters at the top of this file. 
*/

main()
{
  int i,j,no_cws,nn_short,kk_short,byte,error_byte,no_ch_errs,iter,error_flag;
  int received[1000],cw[1000],decode_flag,no_decoder_errors;
  long seed;
  unsigned short int tmp;
  char cw_file[30],file_in[30];
  FILE *fp_out,*fp;
  /**** Storage for random seed ****/
  unsigned short int store[3];
  long int nrand48();
  double x,erand48(),prob_symb_err;
  void srand48();

 /* Generate the Galois field GF(2**mm) */
 printf("Generating the Galois field GF(%d)\n\n",n); 
 generate_gf();

#ifdef DEBUG
  printf("Look-up tables for GF(2**%2d)\n",mm) ;
  printf("i\t\t alpha_to[i]\t index_of[i]\n") ;
  for (i=0; i < n; i++)
   printf("%3d \t\t %#x \t\t %d\n",i,alpha_to[i],index_of[i]) ;
  printf("\n The polynomial-form representation of @**i, 0 <= i < %d \n",n);
  printf("is obtained as follows: It is given by the binary representation of \n");
  printf("the integer alpha_to[i] with LS-Bit corresponding to @**0\n");
  printf("and the next bit to @**1 and so on.\n");
  printf("\n The index of a Galois field element x  whose polynomial-form \n");
  printf("representation is the binary representation of integer \n");
  printf("j, 0 <= j < %d is given by the integer index_of[j]\n",n);
  printf("Therefore, @**j = x \n");
  printf("\n EXAMPLES IN GF(256) \n");
  printf("The polynomial-form of @^8 is the binary representation of the integer alpha_to[8] viz. 0x1d. i.e., @**8 = @**4 + @**3 + @**2 + @**0 \n");
  printf("The index of x whose polynomial-form is 0x72 (integer 155) is index_of[155] = 217, so @**217 = x \n");
	 
#endif
 
 
  printf("Enter b0 = index of the first root \n");
  scanf("%d",&b0);

  gen_poly();

  printf("\n g(x) of the %d-error correcting RS (%d,%d) code: \n",tt,nn,kk);
  printf("Coefficient is the exponent of @ = primitive element\n");
  for (i=0;i <= nn-kk ;i++){
      printf("%d x^%d ",gg[i],i);
      if (i < nn-kk) printf(" + ");
      if (i && (i % 7 == 0)) printf("\n"); /* 8 coefficients per line */
  }
  printf("\n");

  printf("\n Coefficient is the representation in basis {@^%d,...,@,1}\n",mm-1);
  for (i=0;i <= nn-kk ;i++){
      printf("%#x x^%d ",alpha_to[gg[i]],i);
      if (i < nn-kk) printf(" + ");
      if (i && (i % 7 == 0)) printf("\n"); /* 8 coefficients per line */
  }
  printf("\n\n");

#ifdef ENCODE
  printf("Enter length of the shortened code\n");
  scanf("%d",&nn_short);
  if ((nn_short < 2*tt) || (nn_short > nn)){
     printf("Invalid entry %d for shortened length\n",nn_short);
     exit();
  }
  kk_short = kk - (nn-nn_short); /* compute dimension of shortened code */
  printf("The (%d,%d) %d-error correcting RS code\n",nn_short,kk_short,tt);
  printf("Enter the number of codewords desired \n");
  scanf("%d",&no_cws);
  printf("Enter the filename for storing cws\n");
  scanf("%s",cw_file);
  printf("Enter 3 positive integers as seed \n");
  for (i=0;i < 3;i++){ scanf("%hu",&tmp);store[i] = tmp;}
  
  if ((fp_out = fopen(cw_file,"w")) == NULL){
 	printf("Could not open %s\n",cw_file);
	exit();
  }
  
  /**** BEGIN: Encoding random information vectors ****/
     /* Pad with zeros rightmost (kk-kk_short) positions */
  for (j=kk_short;j < kk;j++) data[j] = 0;
  for (i=0;i < no_cws;i++){
      for (j=0;j < (int)kk_short;j++) data[j] = (int) (nrand48(store) % 256);

      encode_rs();

      for (j=0;j < (int)nn_short-(int)kk_short;j++) /* parity bytes first */
	 fprintf(fp_out,"%02x ",bb[j]);
      for (j=0;j < (int)kk_short;j++) /* info bytes last */
	 fprintf(fp_out,"%02x ",data[j]);
      fprintf(fp_out,"\n");
  }
  fclose(fp_out);
#endif

#ifdef DECODE
  printf("Enter length of the shortened code\n");
  scanf("%d",&nn_short);
  if ((nn_short < 2*tt) || (nn_short > nn)){
     printf("Invalid entry %d for shortened length\n",nn_short);
     exit();
  }
  kk_short = kk - (nn-nn_short); /* compute dimension of shortened code */
  printf("The (%d,%d) %d-error correcting RS code\n",nn_short,kk_short,tt);
  printf("Enter a random positive integer\n");
  scanf("%ld",&seed);
  srand48(seed);

  printf("Enter filename from which to read codewords\n"); scanf("%s",file_in);
  printf("Enter no of codewords \n"); scanf("%d",&no_cws);

  if ((fp = fopen(file_in,"r")) == NULL){             
        printf("Could not open %s\n",file_in);
        exit();
  }

  prob_symb_err = (double)tt/ (8.0 * (double)nn_short);
  /* CHANGE above probablity of symbol error if desired */
  /* printf("Probability of symbol error = %6e\n",prob_symb_err);*/

  no_decoder_errors = 0;
  iter = -1;
  while (++iter < no_cws){
     	/* read in codewords from file */
  	for (i=0;i < nn_short;i++){
      	    fscanf(fp,"%02x ",&byte);
      	    cw[i] = byte;
  	}
  	/* printf("The Transmitted Codeword\n");
  	for (i=0;i < nn_short;i++) printf("%02x ",cw[i]); printf("\n");*/

	
#ifndef NO_PRINT
	 printf("\n\n\n\n Transmitting codeword %d \n",iter);
	 printf("Channel caused following errors Location (Error): \n");
#endif
	 no_ch_errs = 0;
	 for (i=0;i < nn_short;i++){
	      x = erand48(store);
	      if (x < prob_symb_err){
		 error_byte = (int) (lrand48() % 256);
		 received[i] = cw[i] ^ error_byte;
		 ++no_ch_errs;
#ifndef NO_PRINT
		 printf("%d (%#x) ",i,error_byte);
#endif
	      }
	      else
		 received[i] = cw[i];
	  }
#ifndef NO_PRINT
	  printf("\n");
	  printf("Channel caused %d errors\n",no_ch_errs);
	  /* for (i=0;i < nn_short;i++) printf("%02x ",received[i]); printf("\n");*/
#endif

	  /* Pad with zeros and decode as if in (255,kk) tt-error correcting code */
          for (i=0;i < nn-kk;i++)  recd[i] = received[i]; /* parity bytes */
	  for (i=nn-kk;i < nn-kk+kk_short;i++) recd[i] = received[i]; /* info bytes */ 
	  for (i=nn-kk+kk_short;i < nn;i++) recd[i] = 0; /* zero padding */
	
	  decode_flag = decode_rs();
#ifndef NO_PRINT
	  printf("decode_rs() returned %d\n",decode_flag);
          error_flag = 0;
          for (i=0; i < kk_short; i++){ 
	      if (recd[i+nn-kk] != cw[i+nn-kk]){
		 error_flag = 1;
		 break;
	      }
	  }
#endif

	  if (decode_flag == 2){
             if (no_ch_errs <= max_errs){
                 printf("%d ch errs  <=  max # correctable errs but \n",no_ch_errs,max_errs);
                 printf("\n DECODER ERROR: deg(lambda) unequal to #roots \n");
                 exit(2); /* DECODER ERROR condition */
             }
	  }
	  else if (decode_flag == 3){
	     if (no_ch_errs <= max_errs){
                printf(" %d ch errs  <= max # correctable errs but \n",no_ch_errs,max_errs);
       		printf("\n deg(lambda) > 2*tt \n");
             	exit(3);/* DECODER ERROR condition */
             }
	  }
	  else{
	     if ((no_ch_errs <= max_errs) && (error_flag == 1)){
		printf("DECODER FAILED TO CORRECT %d ERRORS\n",no_ch_errs);
		++no_decoder_errors;
#ifndef NO_PRINT
		printf("The Transmitted Codeword\n");
		for (i=0;i < nn_short;i++) printf("%02x ",cw[i]); printf("\n");
#endif 
	     }
	  }

   } /* closing iteration loop */
   if (no_decoder_errors > 0)
      printf("The number of decoder errors were %d\n",no_decoder_errors);
   else
      printf("Decoder corrected all occurrences of %d or less errors\n",tt);
#endif

}