bch.c

Go to the documentation of this file.

/* MUX project: bch.c
 *
 * This file contains the functions for the doing the BCH coding of the 
 * header.  
 * Most of the BCH code is taken from :
 *
 *   File:    bch3.c
 *   Title:   Encoder/decoder for binary BCH codes in C (Version 3.1)
 *   Author:  Robert Morelos-Zaragoza (robert@spectra.eng.hawaii.edu)
 *   Date:    August 1994
 *
 * Brendan Dowling 8/31/95 (repitition code)
 * Max Luttrell 4/4/96 (BCH code)
 *
 */

#include <stdio.h>
#include <math.h>

#include "bytes.h"
#include "bits.h"

#include "bch.h"
#include "input.h"

/*  
 * DESCRIPTION OF VARIABLES
 *
 * m = order of the Galois field GF(2**m) 
 * n = 2**m - 1 = size of the multiplicative group of GF(2**m)
 * length = length of the BCH code
 * t = error correcting capability (max. no. of errors the code corrects)
 * d = 2*t + 1 = designed min. distance = no. of consecutive roots of g(x) + 1
 * k = n - deg(g(x)) = dimension (no. of information bits/codeword) of the code
 * p[] = coefficients of a primitive polynomial used to generate GF(2**m)
 * g[] = coefficients of the generator polynomial, g(x)
 * alpha_to [] = log table of GF(2**m) 
 * index_of[] = antilog table of GF(2**m)
 * data[] = information bits = coefficients of data polynomial, i(x)
 * bb[] = coefficients of redundancy polynomial x^(length-k) i(x) modulo g(x)
 * numerr = number of errors 
 * errpos[] = error positions 
 * recd[] = coefficients of the received polynomial 
 * decerror = number of decoding errors (in _message_ positions) 
 *
 */

static int             m, n, length, k, t, d;
static int             p[21];
static int             alpha_to[1048576], index_of[1048576], g[548576];
static int             recd[1048576], data[1048576], bb[548576];

extern int debug;

/* int             numerr, errpos[1024], decerror = 0; */


static void 
generate_gf()
/*
 * Generate field GF(2**m) from the irreducible polynomial p(X) with
 * coefficients in p[0]..p[m].
 *
 * Lookup tables:
 *   index->polynomial form: alpha_to[] contains j=alpha^i;
 *   polynomial form -> index 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;
}

static void 
gen_poly()
/*
 * Compute the generator polynomial of a binary BCH code. Fist generate the
 * cycle sets modulo 2**m - 1, cycle[][] =  (i, 2*i, 4*i, ..., 2^l*i). Then
 * determine those cycle sets that contain integers in the set of (d-1)
 * consecutive integers {1..(d-1)}. The generator polynomial is calculated
 * as the product of linear factors of the form (x+alpha^i), for every i in
 * the above cycle sets.
 */
{
        register int    ii, jj, ll, kaux;
        register int    test, aux, nocycles, root, noterms, rdncy;
        int             cycle[1024][21], size[1024], min[1024], zeros[1024];

        /* Generate cycle sets modulo n, n = 2**m - 1 */
        cycle[0][0] = 0;
        size[0] = 1;
        cycle[1][0] = 1;
        size[1] = 1;
        jj = 1;                 /* cycle set index */
        if (m > 9)  {
                printf("BCH: gen_poly() - Computing cycle sets modulo %d\n", n);
        }
        do {
                /* Generate the jj-th cycle set */
                ii = 0;
                do {
                        ii++;
                        cycle[jj][ii] = (cycle[jj][ii - 1] * 2) % n;
                        size[jj]++;
                        aux = (cycle[jj][ii] * 2) % n;
                } while (aux != cycle[jj][0]);
                /* Next cycle set representative */
                ll = 0;
                do {
                        ll++;
                        test = 0;
                        for (ii = 1; ((ii <= jj) && (!test)); ii++)     
                        /* Examine previous cycle sets */
                          for (kaux = 0; ((kaux < size[ii]) && (!test)); kaux++)
                             if (ll == cycle[ii][kaux])
                                test = 1;
                } while ((test) && (ll < (n - 1)));
                if (!(test)) {
                        jj++;   /* next cycle set index */
                        cycle[jj][0] = ll;
                        size[jj] = 1;
                }
        } while (ll < (n - 1));
        nocycles = jj;          /* number of cycle sets modulo n */

        d = 2 * t + 1;
        /* Search for roots 1, 2, ..., d-1 in cycle sets */
        kaux = 0;
        rdncy = 0;
        for (ii = 1; ii <= nocycles; ii++) {
                min[kaux] = 0;
                test = 0;
                for (jj = 0; ((jj < size[ii]) && (!test)); jj++)
                        for (root = 1; ((root < d) && (!test)); root++)
                                if (root == cycle[ii][jj])  {
                                        test = 1;
                                        min[kaux] = ii;
                                }
                if (min[kaux]) {
                        rdncy += size[min[kaux]];
                        kaux++;
                }
        }
        noterms = kaux;
        kaux = 1;
        for (ii = 0; ii < noterms; ii++)
                for (jj = 0; jj < size[min[ii]]; jj++) {
                        zeros[kaux] = cycle[min[ii]][jj];
                        kaux++;
                }
        k = length - rdncy;

        /* Compute the generator polynomial */
        g[0] = alpha_to[zeros[1]];
        g[1] = 1;               /* g(x) = (X + zeros[1]) initially */
        for (ii = 2; ii <= rdncy; ii++) {
          g[ii] = 1;
          for (jj = ii - 1; jj > 0; jj--)
            if (g[jj] != 0)
              g[jj] = g[jj - 1] ^ alpha_to[(index_of[g[jj]] + zeros[ii]) % n];
            else
              g[jj] = g[jj - 1];
          g[0] = alpha_to[(index_of[g[0]] + zeros[ii]) % n];
        }
        
        if (debug >= 5) {
          printf("BCH: gen_poly() - g(x) = ");
          for (ii = 0; ii <= rdncy; ii++) {
            printf("%d", g[ii]);
            if (ii && ((ii % 50) == 0))
              printf("\n");
          }
          printf("\n");
        }
}


static void 
mencode_bch()
/*
 * Compute redundacy bb[], the coefficients of b(x). The redundancy
 * polynomial b(x) is the remainder after dividing x^(length-k)*data(x)
 * by the generator polynomial g(x).
 */
{
        register int    i, j;
        register int    feedback;

        for (i = 0; i < length - k; i++)
                bb[i] = 0;
        for (i = k - 1; i >= 0; i--) {
                feedback = data[i] ^ bb[length - k - 1];
                if (feedback != 0) {
                        for (j = length - k - 1; j > 0; j--)
                                if (g[j] != 0)
                                        bb[j] = bb[j - 1] ^ feedback;
                                else
                                        bb[j] = bb[j - 1];
                        bb[0] = g[0] && feedback;
                } else {
                        for (j = length - k - 1; j > 0; j--)
                                bb[j] = bb[j - 1];
                        bb[0] = 0;
                }
        }
}


static int mdecode_bch(int correction)
/*
 * Simon Rockliff's implementation of Berlekamp's algorithm.
 *
 * Assume we have received bits in recd[i], i=0..(n-1).
 *
 * Compute the 2*t syndromes by substituting alpha^i into rec(X) and
 * evaluating, storing the syndromes in s[i], i=1..2t (leave s[0] zero) .
 * Then we use the Berlekamp algorithm to find the error location polynomial
 * elp[i].
 *
 * If the degree of the elp is >t, then we cannot correct all the errors, and
 * we have detected an uncorrectable error pattern. We output the information
 * bits uncorrected.
 *
 * If the degree of elp is <=t, we substitute alpha^i , i=1..n into the elp
 * to get the roots, hence the inverse roots, the error location numbers.
 * This step is usually called "Chien's search".
 *
 * If the number of errors located is not equal the degree of the elp, then
 * the decoder assumes that there are more than t errors and cannot correct
 * them, only detect them. We output the information bits uncorrected.
 */
{
        register int    i, j, u, q, t2, count = 0, syn_error = 0;
        int             elp[1026][1024], d[1026], l[1026], u_lu[1026], s[1025];
        int             root[200], loc[200], err[1024], reg[201];
        int             returnval = 1;

        t2 = 2 * t;

        /* first form the syndromes */
        if (debug >= 3) {
          printf("BCH: mdecode_bch() - S(x) = "); 
        }

        for (i = 1; i <= t2; i++) {
          s[i] = 0;
          for (j = 0; j < length; j++)
            if (recd[j] != 0)
              s[i] ^= alpha_to[(i * j) % n];
          if (s[i] != 0) {
                        syn_error = 1; /* set error flag if non-zero syndrome */
                        /* if we are only detecting, bail out */
                        if (!correction) 
                                return 0;
          }
          /*
           * Note:    If the code is used only for ERROR DETECTION, then
           *          exit program here indicating the presence of errors.
           */
          /* convert syndrome from polynomial form to index form  */
          s[i] = index_of[s[i]];
          if (debug >= 3) {
            printf("%3d ", s[i]);
          }
        }
        if (debug >= 3) {
          printf("\n");
        }
        
        if (syn_error) {        /* if there are errors, try to correct them */
                /*
                 * Compute the error location polynomial via the Berlekamp
                 * iterative algorithm. Following the terminology of Lin and
                 * Costello's book :   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*t (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. 
                 */
                /* 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 < t2; 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 {
                        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 < t2; 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] + n - d[q] + elp[q][i]) % n];
                                for (i = 0; i <= l[u]; i++) {
                                        elp[u + 1][i] ^= elp[u][i];
                                        elp[u][i] = index_of[elp[u][i]];
                                }
                        }
                        u_lu[u + 1] = u - l[u + 1];
 
                        /* form (u+1)th discrepancy */
                        if (u < t2) {   
                        /* 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]]) % n];
                          /* put d[u+1] into index form */
                          d[u + 1] = index_of[d[u + 1]];        
                        }
                } while ((u < t2) && (l[u + 1] <= t));
 
                u++;
                if (l[u] <= t) {/* Can correct errors */
                        /* put elp into index form */
                        for (i = 0; i <= l[u]; i++)
                                elp[u][i] = index_of[elp[u][i]];

                        if (debug >= 5) {
                          printf("BCH: mencode_bch() - sigma(x) = ");
                          for (i = 0; i <= l[u]; i++)
                                printf("%3d ", elp[u][i]);
                          printf("\n");
                        } 
                        
                        if (debug >= 3) {
                          printf("BCH: mencode_bch() - Roots: ");
                        }
                        /* Chien search: 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 <= n; i++) {
                                q = 1;
                                for (j = 1; j <= l[u]; j++)
                                        if (reg[j] != -1) {
                                                reg[j] = (reg[j] + j) % n;
                                                q ^= alpha_to[reg[j]];
                                        }
                                if (!q) {       /* store root and error
                                                 * location number indices */
                                        root[count] = i;
                                        loc[count] = n - i;
                                        count++;
                                        if (debug >=3 )
                                          printf("%3d ", n - i);
                                }
                        }
                        if (debug >= 3)
                          printf("\n");
                        if (count == l[u]) {    
                        /* no. roots = degree of elp hence <= t errors */
                                for (i = 0; i < l[u]; i++)
                                        recd[loc[i]] ^= 1;
                        }
                        else{   /* elp has degree >t hence cannot solve */
                                if (debug >=3)
                                  printf("mdecode_bch() - Incomplete decoding: errors detected\n");
                        }
                }
        }
        return returnval;  
}


/*************************************************************/
/* This is the stuff that we can work with.                  */
/*************************************************************/



static void bch_init(bch_type type)
{
  static bch_type old_type=-1;
  int i,j;

  /*
   * let's only initialize stuff if the parameters have changed.  
   * that should speed this up a lot.
   */

  if (old_type != type) {
    /*                             
     * initial stuff.  set p[] = 0's.      
     */
    for (i=0;i<21;i++)
      p[i] = 0;

    /*
     * setup code specific params.
     */
    switch(type) {
    case BCH15_5:  
      /*   (15,5,7) BCH Code */      
      k=5;      t=3;      m = 4;     length = 15;  
      break;
    case BCH15_7:  
      /*   (15,7,5) BCH code */      
      k=7;      t=2;      m = 4;     length = 15;  
      break;
    case BCH15_11: 
      /*   (15,11,3) BCH code */     
      k=11;     t=1;      m = 4;     length = 15;  
      break;
    case BCH31_6:  
      /*   (31,6,15) BCH code */     
      k=6;      t=7;      m = 5;     length = 31;  
      break;
    case BCH31_11: 
      /*   (31,11,11) BCH code */    
      k=11;     t=5;      m = 5;     length = 31;  
      break;
    case BCH63_7:  
      /*   (63,7,31) BCH code */     
      k=7;      t=15;     m = 6;     length = 63;  
      break;
    case BCH63_10: 
      /*   (63,10,37) BCH code */    
      k=10;     t=13;     m = 6;     length = 63;  
      break;
    }

    /*
     * next set n = 2^m-1.  
     */
    n=1;
    for (i=0;i<m;i++)
      n*=2;
    n-=1;

    /*
     * Set length specific params.
     */
    switch(type) {
    case BCH15_5: 
    case BCH15_7: 
    case BCH15_11:      
      /* p(x) = 11001 */
      p[0]=1;      p[1]=1;      p[4]=1;
      break;
    case BCH31_6:
    case BCH31_11:
      /* p(x) = 101001 */
      p[0]=1;      p[2]=1;      p[5]=1;
      break;
    case BCH63_7:
    case BCH63_10:
      /* p(x) = 1100001 */
      p[0]=1;      p[1]=1;      p[6]=1;
      break;
    }
    generate_gf();
    gen_poly();
  }
}




void encode_bch(unsigned long bits, byte *block, bch_type type)
/*
 * external interface to BCH encoding.
 */
{
  int i,j;

  if (debug >= 3) {
    printf("BCH: encode_bch() type = %s\n",  bch_type_str(type));
  }

  if (type == BCHNONE) {
    block[0] = bits & 0xFF;
    return;
  }

  /* 
   * init codec if necessary
   */
  bch_init(type);
 
  /* 
   * put input into encoder.
   */

  for (i=0;i<k;i++) 
    data[i] = getbitint(bits,(16-i)); /* data stores the field in REVERSE */

  /*
   * Run the encoder.
   */
  mencode_bch();   

  /* 
   * Put the redundancy + data into recd[] 
   */
  for (i = 0; i < length - k; i++)
    recd[i] = bb[i];
  for (i = 0; i < k; i++)
    recd[i + length - k] = data[i];
  
  /*
   * recd[] now contains the code.  put this into Block and return 
   */
    for (i=0;i<=(length/8);i++)
      for(j=1;j<=8;j++)
        setbit(&(block[i]),j,recd[(8*i)+j-1]);

}


int decode_bch(unsigned long *bits, byte *block, bch_type type, int correction)
     /*
      *
      * For detection mode,
      * Returns 1 if decoding ok.
      * Returns 0 if decoding not ok.
      *
      * For correction mode,
      * always returns 1
      *
      * correction is 0 if bch codec used for detection, 1 if for correction
      * at this point, only correction mode works properly
      */
{
  int i, j;
  unsigned int r=0;
  int returnval=0;


  if ((correction!=0)&&(correction!=1))
    error("decode_bch","incorrect correction mode");
  
  /*
   * begin the decoding process 
   */
  if (debug >= 3) {
    printf("BCH: decode_bch() - bch_code = %s\n", bch_type_str(type));
  }

  if (type == BCHNONE) {
    *bits = block[0];
    return 1;
  }

  /* 
   * init codec if necessary
   */
  bch_init(type);
 

  /* 
   * copy input to decoder to recd[]
   */
  
  for (i=0;i<m-1;i++) 
    for (j=1;j<8;j++)
      recd[(8*i)+j-1] = getbit(block[i],j);
  
  /*
   * run the BCH decoding algorithm
   */
  returnval = mdecode_bch(correction);  

  /*
   * put output of bch coder into *bits;
   */

  j=1;
  r = 0;
  for (i=length-k;i<length;i++) {
    r += j*recd[i];
    j*=2;
  }
  *bits = r;
 
  if (debug >= 3) {
    printf("BCH: decode_bch() bits = %d\n", *bits);
  }

  return (returnval);
}




bch_type strtobch_type(char *str)
{
  int i;
  i=strlen(str)-1;
  while (i > 0 && (str[i]==' ' || str[i]=='\t'))
    str[i]=0;

  if (equals(str, "BCH15_5"))
    return BCH15_5;
  else if (equals(str, "BCH15_7"))
    return BCH15_7;
  else if (equals(str, "BCH15_11"))
    return BCH15_11;
  else if (equals(str, "BCH31_6"))
    return BCH31_6;
  else if (equals(str, "BCH31_11"))
    return BCH31_11;
  else if (equals(str, "BCH63_7"))
    return BCH63_7;
  else if (equals(str, "BCH63_10"))
    return BCH63_10;
  else if (equals(str, "GOLAY23_12"))
    return GOLAY23_12;
  else
    return BCHNONE;
}

char *bch_type_str(bch_type type)
{
  switch(type) {
  case BCHNONE:  return "None";
  case BCH15_5:  return "(15,5)";
  case BCH15_7:  return "(15,7)";
  case BCH15_11: return "(15,11)";
  case BCH31_6:  return "(31,6)";
  case BCH31_11: return "(31,11)";
  case BCH63_7:  return "(63,7)";
  case BCH63_10: return "(63,10)";
  }
}


int code_length(bch_type b)
{
  switch(b)
    {
    case BCHNONE:
      return 7;
    case BCH15_5: 
    case BCH15_7: 
    case BCH15_11:      
      return 15;
    case BCH31_6:
    case BCH31_11:
      return 31;
    case BCH63_7:
    case BCH63_10:
      return 63;
    case GOLAY23_12:
      return 23;
    }
  return 0;
}


int info_length(bch_type b)
{
  switch(b)
    {
    case BCHNONE: return 8;
    case BCH15_5: return 5;
    case BCH15_7: return 7;
    case BCH15_11: return 11;
    case BCH31_6: return 6;
    case BCH31_11: return 11;
    case BCH63_7: return 7;
    case BCH63_10: return 10;
    }
  return 0;
}

---