rs_eedec_bm_c4p3.c

详细讲述纠错码的书籍

// ------------------------------------------------------------------
// File:    rs_eedec_bm_c4p3.c
// Author:  Robert Morelos-Zaragoza
// Date:    July 20, 2006
//
// 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 BERLEKAMP-MASSEY'S ALGORITHM
//
// Copyright 2002 (c) Robert H. Morelos-Zaragoza. All rights reserved
//
// MODIFIED FOR CHAPTER 4 PROBLEM 3 OF TEXTBOOK.
// ------------------------------------------------------------------

#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, 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[])
{
  m=4;
  length=15;
  red=6;
  init_zero=1;

  k = length - red;
  t = red/2;
  t2 = 2*t;

  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 */

  for (i=0; i<length; i++) 
     recd[i]=0;
  recd[3] = 11;
  recd[6] = 8;
  recd[12] = 3;

  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[] */

}





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;
  p[1] = 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);
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 elp[33][1024], d[1026], l[1026], u_lu[1026], 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 */


#ifdef PRINT_SYNDROME
printf("\ntau =    ");
for (i=0; i<=numera; i++)
  printf("%4d ", tau[i]);
printf("\nforney = ");
#endif


        // Compute FORNEY modified syndrome: