rs_de.c

一份比较有价值的有关ＲＳ编码的文章 欢迎下载

/***************************************************************************
* File:     	RSDeCode.c
* Title:    	Decoder for RS codes in C
* Authors:  	lzg
* Date:     	2003.4
* Input:	mm,kk,b0,Generate the GF(2^mm) and conversion table
*		recd,received data that is Waiting for been decoded
*
* Return:	recd,data that has been decoded
*		return value 
*			0:no error,so correct decode codewords
*			1:correct decode codewords
*			2:cann't find root,so cann't correct decode codewords
*			3:error number > tt ,so dismissed for correcting decode codewords
* Note:		two ways for decodeing
*****************************************************************************/
/*===========================================================================
*Fuction:
*	decode_rs()	performs errors-only decoding of RS codes
*	
*
============================================================================*/

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

#ifndef Gen_GF
#define Gen_GF
//#include "Gen_GF.c"
#endif

#ifndef GFMM_KK_B0
#define GFMM_KK_B0
#define mm		8
#define nn		255
#define kk		243
#define tt		6
#define b0		120
#define nn_short	31
#define kk_short	19
#endif

/* Performs ERRORS-ONLY decoding of RS codes. If decoding is successful,
 writes the codeword into recd[] itself. Otherwise recd[] is unaltered. 
 If channel caused no more than "tt" errors, the tranmitted codeword will
 be recovered.
*/ 

int decode_rs(int *recd)
 {
   register int i,j,u,q ;
   int elp[nn-kk+2][nn-kk], d[nn-kk+2], l[nn-kk+2], u_lu[nn-kk+2], s[nn-kk+1] ;
   int count=0, syn_error=0, root[tt], loc[tt], z[tt+1], err[nn], reg[tt+1] ;

   int *Alpha_to,*Index_of,*gg;
//   FILE *fp_out;
   /*
   int tt;
   tt=(nn-kk)/2;
   */
   /*nn=pow(2,mm) - 1;*/
   Alpha_to = (int *) calloc((nn + 1),sizeof(int));
   Index_of = (int *) calloc((nn + 1),sizeof(int));
   gg = (int *) calloc((nn - kk),sizeof(int));

   /*generate GF(2^mm) and conversion table for RS encode*/
   gen_gf(Alpha_to,Index_of,gg);


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


   /********************************计算伴随式S************************************/
   /* 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]+(long)(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 errors, try and correct */
    {

/* compute the error location polynomial via the Berlekamp iterative algorithm,
   following the terminology of Lin and Costello :   d[u] is the 'mu'th
   discrepancy, where u = 'mu'+ 1 and 'mu' (the Greek letter!) is the step number
   ranging from -1 to 2*tt (see L&C),  l[u] is the
   degree of the elp at that step, and u_l[u] is the difference between the
   step number and the degree of the elp.
  
   The notation is the same as that in Lin and Costello's book; pages 155-156 and 175. */


/*********************************BM叠代法计算错误位置多项式elp******************************/
      /* initialise table entries */
      d[0] = 0 ;           	/* index form elp差值d*/
      d[1] = s[1] ;        	/* index form */
      elp[0][0] = 0 ;      	/* index form 错误位置多项式*/
      elp[1][0] = 1 ;      	/* polynomial form */
      for (i=1; i<nn-kk; i++)
        { 
          elp[0][i] = -1 ;   	/* index form */
          elp[1][i] = 0 ;   	/* polynomial form */
        }
      l[0] = 0 ;		/*elp的度数D*/
      l[1] = 0 ;
      u_lu[0] = -1 ;		/*i-D*/
      u_lu[1] = 0 ;
      u = 0 ;			/*叠代次数*/

      do
      {
        u++ ;
        if (d[u]==-1)		/*差值d为0*/
          { 
            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 */
          {
	/*************************查找满足i-D最大的q值********************************/ 
            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) ;
             } ;
	/********************查找满足i-D最大的q值结束********************************/ 


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


	/* form new elp(x) */
	/******************************计算elp*****************************************/
            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] ;
	/*******************************计算elp结束****************************************/


	/* form (u+1)th discrepancy */
	/********************************计算差值d*****************************************/
        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 */
          }
	/********************************计算差值d结束***********************************/

      } while ((u < nn-kk) && l[u+1] <= tt);
/****************************************叠代结束************************************************/
	
      u++ ;
 
      if (l[u] <= tt)         /* can correct error */
       {

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

	/*****************************钱搜索计算错误位置**************************************/
	/* 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++ ;
              };
          } ;
	/*****************************钱搜索计算错误位置结束***********************************/
      
         if (count == l[u])    /* no. roots = degree of elp hence <= tt errors */
          {
		/* 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 */
            } ;
            
  	/* 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]]]+(long)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 */
               }
            /***************************计算错误值结束******************************/
            }
	    free(Alpha_to);
   	    free(Index_of);
   	    free(gg);
	    return(1);
          }
          else		/* no. roots != degree of elp => >tt errors and cannot solve */
          {    
           	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 */
		}
		free(Alpha_to);
   	    	free(Index_of);
   	    	free(gg);
		return(2);
	  }
       }
      else{         /* elp has degree has degree >tt hence cannot solve */
       	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 */
	}
	free(Alpha_to);
   	free(Index_of);
   	free(gg);
	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]] ;
           else  
	 	recd[i] = 0 ;
 	}
 	free(Alpha_to);
   	free(Index_of);
   	free(gg);
	return(0);
   }
 }