decode_rs.c

这个文件包括wimax中所有的编解码源代码

/*****************************************************************************/
/*   FIle Name : decode_rs.c                                                 */
/*   Description : Decoder of Reed-Solomon FEC                               */
/*                 The RS encoder may NOT be compliant with 802.16.2004 spec.*/
/*                 So this decoder is not be sophisticated enough.           */
/*   author : miffie                                                         */
/*   Date   : sep/20/05                                                      */
/*   Copyright (c) 2005 miffie   All rights reserved.                        */
/*****************************************************************************/

struct binaryset decode_rs( struct binaryset datain, char RR )
/* assume we have received bits grouped into mm-bit symbols in datain.data[i],
   i=0..(nn-1),  and datain.data[i] is index form (ie as powers of alpha).
   We first compute the 2*tt syndromes by substituting alpha**i into rec(X) and
   evaluating, storing the syndromes in s[i], i=1..2tt (leave s[0] zero) .
   Then we use the Berlekamp iteration to find the error location polynomial
   elp[i].   If the degree of the elp is >tt, we cannot correct all the errors
   and hence just put out the information symbols uncorrected. If the degree of
   elp is <=tt, we substitute alpha**i , i=1..n into the elp to get the roots,
   hence the inverse roots, the error location numbers. If the number of errors
   located does not equal the degree of the elp, we have more than tt errors
   and cannot correct them.  Otherwise, we then solve for the error value at
   the error location and correct the error.  The procedure is that found in
   Lin and Costello. For the cases where the number of errors is known to be too
   large to correct, the information symbols as received are output (the
   advantage of systematic encoding is that hopefully some of the information
   symbols will be okay and that if we are in luck, the errors are in the
   parity part of the transmitted codeword).  Of course, these insoluble cases
   can be returned as error flags to the calling routine if desired.   */
 { //decode_rs
   register int i,j,u,q ;
   static int NN = 255 ;
   char tt ;
   char unfixable ;
   short KK, pad ;
   int elp[RR+2][RR], d[RR+2], l[RR+2], u_lu[RR+2], s[RR+1] ;
   int count=0, syn_error=0, root[RR], loc[RR], z[RR+1], err[NN], reg[RR+1] ;
   //int count=0, syn_error=0, root[tt], loc[tt], z[tt+1], err[NN], reg[tt+1] ;
   int alpha_to [NN+1], index_of [NN+1] ;
   int received[NN+1] ;



// Initializations
 printf("...start decode_rs ( %d %d)\n", datain.size , RR ) ;
  /* generate the Galois Field GF(2**mm) */
  generate_gf(&alpha_to[0], &index_of[0]) ;
  //tt is the number of bytes to be corrected RR= 2*tt 
  tt = RR/2 ;

  //pad is the number of padding byte  as prefix
  KK = NN-RR ;
  pad = KK-(datain.size-RR) ;
 printf("KK pad ( %d %d)\n", KK , pad ) ;
  

  /* copy datain.data[i] into sram */
  for (i=0; i<RR; i++)
     received[i] = datain.data[i]&0xff ; //parity
  //add pads
  if(pad>0) for (i=RR; i<RR+pad; i++) received[i] = 0 ; //pads
  //
  for (i=RR+pad; i<NN; i++) 
     received[i] = datain.data[i-pad]&0xff ; //data

  //for (i=0; i<NN; i++) 
  //   printf("received(%d)=%x\n", i, received[i] ) ;

  unfixable = 0 ;

/* first form the syndromes */
   for (i=1; i<=RR; i++)
    { s[i] = 0 ;
      for (j=0; j<NN; j++) {
          s[i] ^= gf_multiply( (received[j]&0xff), alpha_to[(i*j)%NN] ) ;
      }
/* 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]&0xff] ;
    } ;

   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.
*/
      printf("syn_error  ") ;
/* initialise table entries */
      d[0] = 0 ;           /* index form */
      d[1] = s[1] ;        /* index form */
      elp[0][0] = 0 ;      /* index form */
      elp[1][0] = 1 ;      /* polynomial form */
      for (i=1; i<RR; i++)
        { elp[0][i] = -1 ;   /* index form */
          elp[1][i] = 0 ;   /* polynomial form */
        }
      l[0] = 0 ;
      l[1] = 0 ;
      u_lu[0] = -1 ;
      u_lu[1] = 0 ;
      u = 0 ;

      do
      {//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 ;
            //printf("l[%d]=%d\n", u+1, l[u+1] );

/* form new elp(x) */
            for (i=0; i<RR; 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<RR)    /* 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 */
            //printf("d[%d]=%x\n", u+1, d[u+1] ) ;
          }
      } while ((u<RR) && (l[u+1]<=tt)) ; //do

      u++ ;
      if (l[u]<=tt)         /* can correct error */
       {
         //printf("can correct\n") ;
/* 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++ ;
              };
          } ;
         //printf("count=%d l[%d]=%d\n", count, u, l[u]) ;
         if (count==l[u])    /* no. roots = degree of elp hence <= tt errors */
          {
           printf("counti=%d\n", count) ;
/* 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] ;
            } ;

  /* evaluate errors at locations given by error location numbers loc[i] */
           for (i=0; i<NN; i++)
              err[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++) {
                   err[loc[i]] ^= gf_multiply(z[j], alpha_to[(j*root[i])%NN]) ;
              }
              if (err[loc[i]]!=0)
               { err[loc[i]] = index_of[err[loc[i]]&0xff] ;
                 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] ;
                 //printf("loc(%x)=%x %x\n", i, loc[i], err[loc[i]]) ;
                 received[loc[i]] ^= err[loc[i]] ;  /*received[i] must be in polynomial form */
               }
            }
          }
         else { unfixable=1 ;}    /* no. roots != degree of elp => >tt errors and cannot solve */
       }
     else { unfixable=1 ;}    /* elp has degree has degree >tt hence cannot solve */
    }
   else { printf(" no errors found\n" ) ; }    /* no non-zero syndromes => no errors: output received codeword */

    //copy result to datain set
    for (i=0; i<(datain.size-RR); i++)
       datain.data[i]=received[i+RR+pad] ;
       
    //If it's unfixable , set datain.size same as input. 
    if (unfixable==0) datain.size -=  RR ;
    else printf(" unfixable error occured in rs_decoder\n" ) ;
    printf("\n") ;
    return( datain ) ;
 } //decode_rs