rs_eedec_euc.c

The Reed-Somolon code is specified by the finite field, the length (length <= 2^m-1), the numbe

// ------------------------------------------------------------------
// File:    rs_eedec_euc.c
// Author:  Robert Morelos-Zaragoza
// Date:    March 20, 2002
//
// The Reed-Somolon code is specified by the finite field, the length 
// (length <= 2^m-1), the number of redundant symbols (length-k), and 
// the initial zero of the code, init_zero, such that the zeros are:
//    init_zero, init_zero+1, ..., init_zero+length-k-1
//
// Update 8/17/03: Added a line to check for zero in computation of
// erasure polynomial. Many thanks to Matteo Albanese, a graduate
// student at Politecnico di Milano, for pointing this out.
//
// ERRORS AND ERASURES CORRECTION WITH EXTENDED EUCLIDEAN ALGORITHM
//
// Copyright 2002 (c) Robert H. Morelos-Zaragoza. All rights reserved
// ------------------------------------------------------------------

#include <math.h>
#include <stdio.h>
#include <float.h>
#include <limits.h>
#include <stdlib.h>
#define MAX_RANDOM LONG_MAX    // Maximum value of random() 

int i;
int m, n, length, k, t, t2, t21, d, red;
int init_zero;
int p[10];
int alpha_to[1024], index_of[1024], g[1024];
int recd[1024], data[1024], b[1024];
int numerr, errpos[512], errval[512], decerror;
int biterror, error;
char filename[40], name2[40];
int numera;
int era[512], eraval[512];

void read_p(void);
void generate_gf(void);
void gen_poly(void);
void encode_rs(void);
void bpsk_awgn(void);
void decode_rs(void);
int weight(int word);

main(int argc, char *argv[])
{
  // Command line processing
  if (argc != 5)
  {
  printf("\nSimulation of RS codes \n");
  printf("Copyright 2002 (c) Robert Morelos-Zaragoza\n\n");
  printf("Usage: %s m length red init_zero \n", argv[0]);
  printf(" - m is the order of GF(2^m)\n");
  printf(" - length is the length of the RS code\n");
  printf(" - red = length - dimension\n");
  printf(" - init_zero = inital zero of the code\n");
  exit(0);
  }

  sscanf(argv[1],"%d", &m);
  sscanf(argv[2],"%d", &length);
  sscanf(argv[3],"%d", &red);
  sscanf(argv[4],"%d", &init_zero);

  k = length - red;
  t = red/2;
  t2 = 2*t;
  t21 = t2+1;

  read_p();        // Read m */
  generate_gf();   // Construct the Galois Field GF(2^m) */
  printf("--> This is an RS(%d,%d,%d) code over GF(2^%d), t = %d\n",
          length, k, red+1, m, t);
  printf("    with zeros: ");
  for (i=0;i<t2;i++) printf("%d ", init_zero+i);
  printf("\n\n");

  gen_poly();      // Compute the generator polynomial of RS code */

  // Randomly generate DATA */
  for (i = 0; i < k; i++)
    data[i] = (random()>>10) % n;

  encode_rs();    // encode data */

  // Codeword is c(X) = data(X)*X**(length-k)+ b(X)
  for (i = 0; i < length - k; i++) recd[i] = b[i];
  for (i = 0; i < k; i++) recd[i + length - k] = data[i];

  for (i = 0; i < length; i++) recd[i] = 0;

  // Introduce errors and erasures at will ...

  /* CHANNEL NOISE: Read errors and erasures */
  printf("Number of errors, e: "); scanf("%d", &numerr);
  if (numerr)
    for (i=0; i<numerr; i++)
      {
      printf("position: "); scanf("%d", &errpos[i]); // as power of alpha
      printf("value:    "); scanf("%d", &errval[i]); // as vector
      recd[errpos[i]] ^= errval[i];
      }
  printf("Number of erasures, b: "); scanf("%d", &numera);
  if (numera)
    for (i=0; i<numera; i++)
      {
      printf("position: "); scanf("%d", &era[i]); // as power of alpha
      printf("value:    "); scanf("%d", &eraval[i]); // as vector
      recd[era[i]] ^= eraval[i];
      }


  printf("\n\nRec =");
  for (i=0; i<length; i++) {
     printf("%4d ", index_of[recd[i]]);
     }
  printf("\n     ");
  for (i=0; i<length; i++) {
     printf("%4d ", recd[i]);
     }
  printf("\n");

  decode_rs();   // DECODE received codeword recv[] */

  // DECODING ERRORS? we compare only the data portion
  decerror = 0;
  biterror = 0;
  for (i = length-k; i < length; i++)
    if (data[i-length+k] != recd[i])
      {
      decerror++;
      biterror += weight(data[i-length+k]- recd[i]);
      }
}


void read_p()
//      Read m, the degree of a primitive polynomial p(x) used to compute the
//      Galois field GF(2**m). Get precomputed coefficients p[] of p(x). Read
//      the code length.
{
  int i, ninf;

  printf("\nSimulation of RS codes \n");
  printf("Copyright 2002 (c)Robert Morelos-Zaragoza. All rights reserved.\n\n");
  for (i=1; i<m; i++)
    p[i] = 0;
  p[0] = p[m] = 1;
  if (m == 2)             p[1] = 1;
  else if (m == 3)        p[1] = 1;
  else if (m == 4)        p[3] = 1;
  // else if (m == 4)        p[1] = 1;  // Commented out to match example p. 68
  else if (m == 5)        p[2] = 1;
  else if (m == 6)        p[1] = 1;
  else if (m == 7)        p[1] = 1;
  else if (m == 8)        p[4] = p[5] = p[6] = 1;
  else if (m == 9)        p[4] = 1;
  else if (m == 10)       p[3] = 1;
  else if (m == 11)       p[2] = 1;
  else if (m == 12)       p[3] = p[4] = p[7] = 1;
  else if (m == 13)       p[1] = p[3] = p[4] = 1;
  else if (m == 14)       p[1] = p[11] = p[12] = 1;
  else if (m == 15)       p[1] = 1;
  else if (m == 16)       p[2] = p[3] = p[5] = 1;
  else if (m == 17)       p[3] = 1;
  else if (m == 18)       p[7] = 1;
  else if (m == 19)       p[1] = p[5] = p[6] = 1;
  else if (m == 20)       p[3] = 1;
  printf("Primitive polynomial of GF(2^%d), (LSB first)   p(x) = ",m);
  n = 1;
  for (i = 0; i <= m; i++) 
    {
      n *= 2;
      printf("%1d", p[i]);
    }
  printf("\n");
  n = n / 2 - 1;
}



void generate_gf()
// generate GF(2^m) from the irreducible polynomial p(X) in p[0]..p[m]
//
// lookup tables:  log->vector form           alpha_to[] contains j=alpha**i;
//                 vector form -> log form    index_of[j=alpha**i] = i
// alpha=2 is the primitive element of GF(2^m)
{
register int i, mask; 
  mask = 1; 
  alpha_to[m] = 0; 
  for (i=0; i<m; i++)
   { 
     alpha_to[i] = mask; 
     index_of[alpha_to[i]] = i; 
     if (p[i]!=0)
       alpha_to[m] ^= mask; 
     mask <<= 1; 
   }
  index_of[alpha_to[m]] = m; 
  mask >>= 1; 
  for (i=m+1; i<n; i++)
   { 
     if (alpha_to[i-1] >= mask)
        alpha_to[i] = alpha_to[m] ^ ((alpha_to[i-1]^mask)<<1); 
     else alpha_to[i] = alpha_to[i-1]<<1; 
     index_of[alpha_to[i]] = i; 
   }
  index_of[0] = -1; 

#ifdef PRINT_GF
printf("Table of GF(%d):\n",n+1);
printf("----------------------\n");
printf("   i\tvector \tlog\n");
printf("----------------------\n");
for (i=0; i<=n; i++)
printf("%4d\t%4x\t%4d\n", i, alpha_to[i], index_of[i]);
#endif
 }


void gen_poly()
// Compute the generator polynomial of the t-error correcting, length
// n=(2^m -1) Reed-Solomon code from the product of (X+alpha^i), for
// i = init_zero, init_zero + 1, ..., init_zero+length-k-1
{
register int i,j; 

   g[0] = alpha_to[init_zero];  //  <--- vector form of alpha^init_zero
   g[1] = 1;     // g(x) = (X+alpha^init_zero)
   for (i=2; i<=length-k; i++)
    { 
      g[i] = 1;
      for (j=i-1; j>0; j--)
        if (g[j] != 0)  
          g[j] = g[j-1]^ alpha_to[(index_of[g[j]]+i+init_zero-1)%n]; 
        else 
          g[j] = g[j-1]; 
      g[0] = alpha_to[(index_of[g[0]]+i+init_zero-1)%n];
    }
   // convert g[] to log form for quicker encoding 
   for (i=0; i<=length-k; i++)  g[i] = index_of[g[i]]; 
#ifdef PRINT_POLY
printf("Generator polynomial (independent term first):\ng(x) = ");
for (i=0; i<=length-k; i++) printf("%5d", g[i]);
printf("\n");
#endif
}


void encode_rs()
// Compute the 2t parity symbols in b[0]..b[2*t-1]
// data[] is input and b[] is output in polynomial form.
// Encoding is done by using a feedback shift register with connections
// specified by the elements of g[].
{
   register int i,j; 
   int feedback; 

   for (i=0; i<length-k; i++)   
     b[i] = 0; 
   for (i=k-1; i>=0; i--)
    {
    feedback = index_of[data[i]^b[length-k-1]]; 
      if (feedback != -1)
        { 
        for (j=length-k-1; j>0; j--)
          if (g[j] != -1)
            b[j] = b[j-1]^alpha_to[(g[j]+feedback)%n]; 
          else
            b[j] = b[j-1]; 
          b[0] = alpha_to[(g[0]+feedback)%n]; 
        }
       else
        { 
        for (j=length-k-1; j>0; j--)
          b[j] = b[j-1]; 
        b[0] = 0; 
        }
    }
}



void decode_rs()
{
   register int i,j,u,q;
   int qt[513], r[129][513];
   int b[12][513];
   int degr[129], degb[129];
   int elp[129][1024], l[1], s[1025], forney[1025];
   int count=0, syn_error=0, tau[512], root[512], loc[512], z[513]; 
   int err[1024], reg[513], aux[513], omega[1025], phi[1025], phiprime[1025];
   int degphi, ell, temp;

   // Compute the syndromes

#ifdef PRINT_SYNDROME
   printf("\ns =         0 ");
#endif

   for (i=1; i<=t2; i++)
    { 
      s[i] = 0; 
      for (j=0; j<length; j++)
        if (recd[j]!=0)
          s[i] ^= alpha_to[(index_of[recd[j]]+(i+init_zero-1)*j)%n];
      // convert syndrome from vector form to log form  */
      if (s[i]!=0)  
        syn_error=1;         // set flag if non-zero syndrome => error
      //
      // Note:    If the code is used only for ERROR DETECTION, then
      //          exit program here indicating the presence of errors.
      //
      s[i] = index_of[s[i]]; 

#ifdef PRINT_SYNDROME
   printf("%4d ", s[i]);
#endif

    }

   if (syn_error)       // if syndromes are nonzero then try to correct
    {

     s[0] = 0;  // S(x) = 1 + s_1x + ...

      // TO HANDLE ERASURES

      if (numera)
      // if erasures are present, compute the erasure locator polynomial, tau(x)
        {
          for (i=0; i<=t2; i++) 
            { tau[i] = 0; aux[i] = 0;}
          
          aux[1] = alpha_to[era[0]];
          aux[0] = 1;       // (X + era[0]) 

          if (numera>1)

            for (i=1; i<numera; i++)
              {
              p[1] = era[i];
              p[0] = 0;
              for (j=0; j<2; j++)
                for (ell=0; ell<=i; ell++)
                  // Line below added 8/17/2003
                  if ((p[j] !=-1) && (aux[ell]!=0))
                    tau[j+ell] ^= alpha_to[(p[j]+index_of[aux[ell]])%n];
              if (i != (numera-1))
              for (ell=0; ell<=(i+1); ell++)
                {
                aux[ell] = tau[ell];
                tau[ell] = 0;
                }
              }

          else {
            tau[0] = aux[0]; tau[1] = aux[1];
          }

        // Put in index (log) form
        for (i=0; i<=numera; i++)
          tau[i] = index_of[tau[i]]; /* tau in log form */