new_rs_erasures.c

一些纠错编码算法的源代码

      u_lu[0] = -1 ;
      u_lu[1] = 0 ;
      u = 0 ;

      do
      {
        u++ ;
        if (d[u]==-1)
          { l[u+1] = l[u] ;
            for (i=0; i<=l[u]; i++)
             {  elp[u+1][i] = elp[u][i] ;
                elp[u][i] = index_of[elp[u][i]] ;
             }
          }
        else
/* search for words with greatest u_lu[q] for which d[q]!=0 */
          { q = u-1 ;
            while ((d[q]==-1) && (q>0)) q-- ;
/* have found first non-zero d[q]  */
            if (q>0)
             { j=q ;
               do
               { j-- ;
                 if ((d[j]!=-1) && (u_lu[q]<u_lu[j]))
                   q = j ;
               }while (j>0) ;
             } ;

/* have now found q such that d[u]!=0 and u_lu[q] is maximum */
/* store degree of new elp polynomial */
            if (l[u]>l[q]+u-q)  l[u+1] = l[u] ;
            else  l[u+1] = l[q]+u-q ;

/* form new elp(x) */
            for (i=0; i<nn-kk; i++)    elp[u+1][i] = 0 ;
            for (i=0; i<=l[q]; i++)
              if (elp[q][i]!=-1)
                elp[u+1][i+u-q] = alpha_to[(d[u]+nn-d[q]+elp[q][i])%nn] ;
            for (i=0; i<=l[u]; i++)
              { elp[u+1][i] ^= elp[u][i] ;
                elp[u][i] = index_of[elp[u][i]] ;  /*convert old elp value to index*/
              }
          }
        u_lu[u+1] = u-l[u+1] ;

/* form (u+1)th discrepancy */
        if (u<nn-kk)    /* no discrepancy computed on last iteration */
          {
            if (s[u+1]!=-1)
                   d[u+1] = alpha_to[s[u+1]] ;
            else
              d[u+1] = 0 ;
            for (i=1; i<=l[u+1]; i++)
              if ((s[u+1-i]!=-1) && (elp[u+1][i]!=0))
                d[u+1] ^= alpha_to[(s[u+1-i]+index_of[elp[u+1][i]])%nn] ;
            d[u+1] = index_of[d[u+1]] ;    /* put d[u+1] into index form */
          }
      } while ((u < nn-kk) && (l[u+1] <= tt)) ;


      u++ ;
      if (l[u] <= tt)         /* can correct error */
       {
#ifdef DECODER_DEBUG
    printf("---------------------------------------------------------------\n");
      printf("The Final Error Locator Polynomial Lambda(x) in polynomial-form is: \n");
      for (i=0;i < l[u];i++){
	  printf("%#x x^%d + ",elp[u][i],i);
          if (i && (i % 7 == 0)) printf("\n"); /* 8 coefficients per line */
      }
      printf("%#x x^%d",elp[u][i],i);
      printf("\n");
#endif

/* put elp into index form */
         for (i=0; i<=l[u]; i++)   elp[u][i] = index_of[elp[u][i]] ;

#ifdef DECODER_DEBUG
      printf("The Final Error Locator Polynomial Lambda(x) in index-form is: \n");
      for (i=0;i < l[u];i++){
          printf("%d x^%d + ",elp[u][i],i);
          if (i && (i % 7 == 0)) printf("\n"); /* 8 coefficients per line */
      }  
      printf("%d x^%d",elp[u][i],i);
      printf("\n");
    printf("---------------------------------------------------------------\n");
#endif


/* find roots of the error location polynomial */
         for (i=1; i <= l[u]; i++)
           reg[i] = elp[u][i] ;
         count = 0 ;
         for (i=1; i <= nn; i++)
          {  q = 1 ;
             for (j=1; j<=l[u]; j++)
              if (reg[j]!=-1)
                { reg[j] = (reg[j]+j)%nn ;
                  q ^= alpha_to[reg[j]] ;
                } ;
             if (!q)        /* store root and error location number indices */
              { root[count] = i;
                loc[count] = nn-i ;
                count++ ;
              };
          } ;
#ifndef NO_PRINT
	  printf("Computed Error Locations\n");
	  for (i=0;i < l[u];i++) printf("%d ",loc[i]);
	  printf("\n");
#endif
         if (count == l[u])    /* no. roots = degree of elp hence <= tt errors */
          {
#ifdef DECODER_DEBUG
 	printf("\n No of roots of Lambda(x) found by Chien search equals with \n");
 	printf("the degree of Lambda(x). Hence channel presumably caused <= %d errors\n",tt);
	printf("The roots found are (in index-form) : \n");
	for (i=0;i < count;i++)
	    printf("%d ",root[i]);
        printf("\n");
#endif 
/* form polynomial z(x) */
           for (i=1; i <= l[u]; i++)        /* Z[0] = 1 always - do not need */
            { if ((s[i]!=-1) && (elp[u][i]!=-1))
                 z[i] = alpha_to[s[i]] ^ alpha_to[elp[u][i]] ;
              else if ((s[i]!=-1) && (elp[u][i]==-1))
                      z[i] = alpha_to[s[i]] ;
                   else if ((s[i]==-1) && (elp[u][i]!=-1))
                          z[i] = alpha_to[elp[u][i]] ;
                        else
                          z[i] = 0 ;
              for (j=1; j<i; j++)
                if ((s[j]!=-1) && (elp[u][i-j]!=-1))
                   z[i] ^= alpha_to[(elp[u][i-j] + s[j])%nn] ;
              z[i] = index_of[z[i]] ;         /* put into index form */
            } ;

#ifdef DECODER_DEBUG
    printf("---------------------------------------------------------------\n");
      printf(" \n The Final Error Evaluator Polynomial Omega(x) in index-form is: \n");
      printf("0 x^0 + ");
      for (i=1;i <= l[u];i++){
          printf("%d x^%d + ",z[i],i);
          if (i && (i % 7 == 0)) printf("\n"); /* 8 coefficients per line */
      }
      printf("%d x^%d",z[i],i);
      printf("\n");
#endif

#ifdef DECODER_DEBUG
      printf(" \n The Final Error Evaluator Polynomial Omega(x) in polynomial-form is: \n");
      printf("%#x x^%d + ",1,0);
      for (i=1;i <= l[u];i++){
          printf("%#x x^%d + ",z[i],i);
          if (i && (i % 7 == 0)) printf("\n"); /* 8 coefficients per line */
      }  
      printf("%#x x^%d",z[i],i);
      printf("\n");
    printf("---------------------------------------------------------------\n");
#endif

  /* evaluate errors at locations given by error location numbers loc[i] */
           for (i=0; i<nn; i++)
             { err[i] = 0 ;
               if (recd[i]!=-1)        /* convert recd[] to polynomial form */
                 recd[i] = alpha_to[recd[i]] ;
               else  recd[i] = 0 ;
             }
           for (i=0; i < l[u]; i++)    /* compute numerator of error term first */
            { err[loc[i]] = 1;       /* accounts for z[0] */
              for (j=1; j<=l[u]; j++){
                if (z[j]!=-1)
                  err[loc[i]] ^= alpha_to[(z[j]+j*root[i])%nn] ;
	      } /* z(x) evaluated at X(l)**(-1) */
              if (err[loc[i]] != 0) /* term X(l)**(1-b0) */
		 err[loc[i]] = alpha_to[(index_of[err[loc[i]]]+root[i]*(b0+nn-1))%nn];
              if (err[loc[i]]!=0)
               {
		 err[loc[i]] = index_of[err[loc[i]]] ;
                 q = 0 ;     /* form denominator of error term */
                 for (j=0; j<l[u]; j++)
                   if (j!=i)
                     q += index_of[1^alpha_to[(loc[j]+root[i])%nn]] ;
                 q = q % nn ;
                 err[loc[i]] = alpha_to[(err[loc[i]]-q+nn)%nn] ;
                 recd[loc[i]] ^= err[loc[i]] ;  /*recd[i] must be in polynomial form */
               }
            }
#ifdef DECODER_DEBUG
    		printf("---------------------------------------------------------------\n");
		printf("The computed Error Value Bytes are (in polynomial-form) :\n");
		printf("Location (Error) \n");
		for (i=0;i < l[u];i++){
		    printf("%d (%#x) ",loc[i],err[loc[i]]);
		}
		printf("\n");
    		printf("---------------------------------------------------------------\n");
#endif

	    return(1);
          }
          else{    /* no. roots != degree of elp => >tt errors and cannot solve */
#ifdef DECODER_DEBUG
		printf("No of roots of Lambda(x) found by Chien search\n");
		printf("does not equal the degree of the Lambda(x).\n");
		printf("Hence channel has definitely caused >= %d errors\n",tt);
#endif
           	for (i=0; i<nn; i++){        /* could return error flag if desired */
               	   if (recd[i]!=-1)        /* convert recd[] to polynomial form */
                 	recd[i] = alpha_to[recd[i]] ;
                   else  
			recd[i] = 0 ;     /* just output received word as is */
		}
		return(2);
	  }
       }
      else{         /* elp has degree has degree >tt hence cannot solve */
#ifdef DECODER_DEBUG
		printf("Degree of the Lambda(x) > %d. \n",tt);
		printf("Hence channel has definitely caused >= %d errors\n",tt);
#endif
       	for (i=0; i<nn; i++){       /* could return error flag if desired */
            if (recd[i]!=-1)        /* convert recd[] to polynomial form */
                recd[i] = alpha_to[recd[i]] ;
            else  
		recd[i] = 0 ;     /* just output received word as is */
	}
	return(3);
      }
#ifdef DECODER_DEBUG
    printf("+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n");
#endif
   }
   else{       /* no non-zero syndromes => no errors: output received codeword */
#ifdef DECODER_DEBUG
    printf("Received vector is a codeword. Channel has presumably caused no error\n");
    printf("+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n");
#endif
       	for (i=0; i < nn; i++){
       	   if (recd[i] != -1)        /* convert recd[] to polynomial form */
         	recd[i] = alpha_to[recd[i]] ;
           else  
	 	recd[i] = 0 ;
 	}
	return(0);
   }
 }

/* Performs ERRORS+ERASURES decoding of RS codes. If decoding is successful, writes
the codeword into recd[] itself. Otherwise recd[] is unaltered except in the
erased positions where it is set to zero. 
  First "no_eras" erasures are declared by the calling program. Then, the 
maximum # of errors correctable is t_after_eras = floor((2*tt-no_eras)/2). If the number
of channel errors is not greater than "t_after_eras" the transmitted codeword
 will be recovered. Details of algorithm can be found
in R. Blahut's "Theory ... of Error-Correcting Codes". */

int eras_dec_rs(eras_pos,no_eras)
int eras_pos[2*tt],no_eras;
{
  register int i,j,r,u,q,tmp,tmp1,tmp2,num1,num2,den,pres_root,pres_loc;
  int phi[2*tt+1],tmp_pol[2*tt+1]; /* The erasure locator in polynomial form */
  int U,syn_err,discr_r,deg_phi,deg_lambda,L,deg_omega,t_after_eras;
  int lambda[2*tt+1],s[2*tt+1],lambda_pr[2*tt+1];/* Err+Eras Locator poly and syndrome poly */ 
  int b[2*tt+1],T[2*tt+1],omega[2*tt+1];
  int syn_error=0, root[2*tt], z[tt+1], err[nn], reg[2*tt+1] ;
  int loc[2*tt],count = 0;

/* Maximum # ch errs correctable after "no_eras" erasures */
   t_after_eras = (int)floor((2.0*tt-no_eras)/2.0);

/* Compute erasure locator polynomial phi[x] */
  for (i=0;i < nn-kk+1;i++) tmp_pol[i] = 0;
  for (i=0;i < nn-kk+1;i++) phi[i] = 0; 
  if (no_eras > 0){
      phi[0] = 1; /* index form */
      phi[1] = alpha_to[eras_pos[0]];
      for (i=1;i < no_eras;i++){
          U = eras_pos[i];
	  for (j=1;j < i+2;j++){
	     tmp1 = index_of[phi[j-1]];
	     tmp_pol[j] = (tmp1 == -1)?  0 : alpha_to[(U+tmp1)%nn];
	  }
	  for (j=1;j < i+2;j++)
	      phi[j] = phi[j]^tmp_pol[j];
      }
  /* put phi[x] in index form */
     for (i=0;i < nn-kk+1;i++) phi[i] = index_of[phi[i]];
#ifdef ERASURE_DEBUG
/* find roots of the erasure location polynomial */
         for (i=1; i <= no_eras; i++) reg[i] = phi[i] ;
         count = 0 ;
         for (i=1; i <= nn; i++)
          {  q = 1 ;
             for (j=1; j <= no_eras; j++)
              if (reg[j]!=-1)
                { reg[j] = (reg[j]+j)%nn ;
                  q ^= alpha_to[(reg[j])%nn] ;
                } ;
             if (!q)        /* store root and error location number indices */
              { root[count] = i;
                loc[count] = nn-i ;
                count++ ;
              };
          } ;
	  if (count != no_eras){ printf("\n phi(x) is WRONG\n"); exit(1);}
#ifndef NO_PRINT 
          printf("\n Erasure positions as determined by roots of Eras Loc Poly:\n");
          for (i=0;i < count;i++) printf("%d ",loc[i]);
	  printf("\n");
#endif
#endif
  }

/* recd[] is in polynomial form, convert to index form */
   for (i=0;i < nn;i++) recd[i] = index_of[recd[i]];

/* first form the syndromes; i.e., evaluate recd(x) at roots of g(x) namely
 @**(b0+i), i = 0, ... ,(2*tt-1) */
   for (i=1; i <= nn-kk; i++)
    { s[i] = 0 ;
      for (j=0; j < nn; j++)
        if (recd[j] != -1)
          s[i] ^= alpha_to[(recd[j]+(b0+i-1)*j)%nn] ;      /* recd[j] in index form */
/* convert syndrome from polynomial form to index form  */
      if (s[i] != 0)  syn_error = 1 ;   /* set flag if non-zero syndrome => error */
      s[i] = index_of[s[i]] ;
    };


   if (syn_error){ /* if syndrome is zero, modified recd[] is a codeword */
/* Begin Berlekamp-Massey algorithm to determine error+erasure locator polynomial */
   	r = no_eras;
   	deg_phi = no_eras;
   	L = no_eras;
   	if (no_eras > 0){
        	/* Initialize lambda(x) and b(x) (in poly-form) to phi(x) */
        	for (i=0;i < deg_phi+1;i++) lambda[i] = (phi[i] == -1)? 0 : alpha_to[phi[i]];
        	for (i=deg_phi+1;i < 2*tt+1;i++) lambda[i] = 0;
        	deg_lambda = deg_phi;
        	for (i=0;i < 2*tt+1;i++) b[i] = lambda[i];
   	}
   	else{
		lambda[0] = 1;
		for (i=1;i < 2*tt+1;i++) lambda[i] = 0;
		for (i=0;i < 2*tt+1;i++) b[i] = lambda[i];
   	}
	while (++r <= 2*tt){ /* r is the step number */
	        /* Compute discrepancy at the r-th step in poly-form */
	        discr_r = 0;
	        for (i=0;i < 2*tt+1;i++){
		    if ((lambda[i] != 0) && (s[r-i] != -1)){
		       tmp = alpha_to[(index_of[lambda[i]]+s[r-i])%nn];
		       discr_r ^= tmp;
		    }
	        }