The need for reliable, high-speed transmission of video, audio, and data applications has never been greater. One of the keys to reliable communications is error recovery--and Bose-Chaudhuri-Hocquenghem (BCH) and Reed-Solomon (RS) codes are two powerful approaches to error-control coding. RS, in fact, is specified in virtually every proposed digital-TV transmission standard that tolerates less than one uncorrected error event per hour. At a multiplexed transmission rate of about 30 Mbits/sec, this corresponds to a bit-error rate (BER) in the order of 10. With transmission channels typically having a BER of 10e-3 to 10e-4, powerful error-correction techniques are obviously required.
In this article, I'll examine BCH and RS codes, using them to encode data in such a way that data errors can be recovered. In doing so, I'll present GALOIS.EXE, a C program for constructing Galois fields and BCH/RS generators and for encoding/ decoding of arbitrary files. The program also includes a feature to "scramble" an encoded file with a specified number of symbol errors per block to simulate error correctability. The complete package, including programmer notes, is available electronically; see the top of the first page.
Polynomials and Galois Fields
As Figure 1 illustrates, cyclic-linear block codes, like those discussed in this article, convert blocks of L symbols at the encoder input into blocks of N symbols (N>L). They're called "cyclic" because each shift of valid N-length code word yields another code word. For these codes, the input/output relation can easily be expressed using polynomials. The "generator polynomial" g(D) completely characterizes a particular code.
Figure 1
Finite fields (or Galois fields)have p elements denoted as GF(p). The field's elements are denoted by 0,1,_,p--1. On each pair of elements, both addition and multiplication operations are defined. If the fields have a prime number of arguments, addition and multiplication operations are performed the normal way, with the result being modulo p. One reason for the presence of the implied modulo operation is to keep the result of adding/multiplying two field elements contained in the same GF.
However, there exist fields with a non-prime number of arguments. In such cases, the number of GF elements is a power of a prime number, so the field can be expressed as GF(pn) with p prime. Here I'll refer to GF(p) as the "symbol field" and GF(pn) as the "locator field." You can express an element of GF(pn) as a row of n GF(p) elements and hence as a polynomial of a degree no larger than n-1, with each coefficient value belonging to GF(p).
How do you construct GF(pn)? The addition of two polynomials of degree n--1 yields another polynomial of degree n--1 or less. This is not the case for multiplication. The resulting polynomial may not be a field element. However, if you reduce the polynomial by performing the modulo operation using a degree-n polynomial p(D), the result will be of degree n--1 or less. To construct a valid field, you choose p(D) to be a "primitive" polynomial. For GF(32), such a polynomial is D2+D+2. Table 1 lists the representation of the elements of this Galois field. Note that you can represent a GF(pn) element ![[alpha]](http://twimgs.com/ddj/ddj/images/ddj9412d/alpha14.gif)
i (i=0..pni modulo g(D), g(D) being the primitive polynomial for primitive element . Since the primitive polynomial is of degree n and is irreducible, you can find a polynomial for generating GF(pn) by taking an irreducible polynomial of this degree and calculating Di mod g(D) for i=0,_,pn--2. If this gives you pn-1 (all elements of GF(pn) except 0), g(D) is primitive.
Table 1
All operations are performed on coefficient values belonging to the already-constructed symbol field. Listing One is the code for searching for such a primitive polynomial (there's always at least one for each GF). Afterwards, addition and multiplication operations on each pair of GF(pn) elements are defined by means of the polynomial representation of these operands: by adding or multiplying the corresponding polynomials modulo g(D). The result is defined to be the polynomial representation of the GF element corresponding to the sum or product of both operands. GALOIS.EXE stores this value into a table defining the field's operations.
Listing One
/* Finds a primitive polynomial of GF(p^n) */
void FindPrimitivePolynomial(void) {
int k,d;
int *pi1; /* array of (n+1) integers, each indicating */
/* a number between 0 and p-1 -> generation */
/* of all polynomials of degree n with */
/* coeff. in GF(p) */
int *pi2; /* idem for polynomials of degree <=n-1 */
bool rem0;
d=pg_l->Exp; /* degree of primitive polynomial */
... (allocate and init to zero pp_prim, pp_aux, pi1 & pi2) ...
pi1[d]=1; /* force it to be a degree-q polynomial */
pp_prim->Deg=d;
do{ /* loop for generation of all degree-d polynomials */
for(k=d;k>=0;k--)
pp_prim->c[k]=pi1[k];
if(d!=1){ /* all degree 1 polynomials are irreducible */
memset((void*)pi2,0,sizeof(int)*d);
pi2[1]=1; /* degree 1 polynomial */
do{ /* loop for generation of all degree 0..q-1 pol. */
for(k=d-1;k>=0;k--)
pp_aux->c[k]=pi2[k];
DEGREE(pp_aux);
rem0=DivPoly(pp_prim,pp_aux,NULL,NULL,po_s);
if(rem0) break; /* not irreducible */
}
while(update(pi2,po_s->MaxNo,d-1));
}
else
rem0=0;
if(!rem0) /* found an irreducible polynomial */
if(is_primitive())
break;
}
while(update(pi1,po_s->MaxNo,d));
... (free pp_aux, pi1, pi2) ...
}
/* Gets increment for (d+1) integers, each having maximum value p-1 */
static bool update(int *pi,int p,int q) {
int i,j,max;
max=p-1;
j=0;
while(pi[j]!=max && j<q) j++;
i=j;
while(pi[i]==max && pi[i+1]==max && i<q) i++;
if(i==q && j==0)
return FALSE;
if(j)
pi[0]++;
else {
pi[i+1]++;
for(j=0;j<=i;j++)
pi[j]=0;
}
return TRUE;
}
/* Checks for an irreducible polynomial to be primitive. If it is found to be
primitive, the function returns TRUE and <pg> contains all GF-elements. */
static bool is_primitive(void) {
int p,i,k,l;
pPOLYNOM pp_rem=NULL;
bool is_prim=TRUE;
p=pp_prim->Deg;
...(allocate pp_aux, pp_rem) ...
...(init pg_l->pe[0] to 0
pg_l->pe[1] to 1
remaining entries to 0) ...
for(i=2;i<pg_l->MaxNo && is_prim;i++) { /* `^1,..,`^(p^q-2) */
SET_POLY(pp_aux,i-1,1); /* D^i <-> `^i */
DivPoly(pp_aux,pp_prim,NULL,pp_rem,po_s);
for(k=p-1;k>=0;k--)
pg_l->pe[i].Alpha[k] = (byte)pp_rem->c[k];
for(l=0;l<i;l++) {
for(k=0;k<p;k++) {
if(pg_l->pe[l].Alpha[k] != pg_l->pe[i].Alpha[k])
break;
}
if(k==p) { /* then pe[l] == pe[i] -> not primitive */
is_prim=FALSE;
break;
}
}
}
...(free pp_aux, pp_rem) ...
return(is_prim);
}
Listing Two provides the syntax used in the Galois program to store these and additional one-dimensional tables for additive and multiplicative inverses. Table 2 shows addition and multiplication tables for GF(22), constructed out of GF(2) using D2+D+1 as a primitive polynomial over GF(2). Note that operations in GF(4) are not ordinary addition and multiplication modulo 4 (4 is not prime!). You could apply this procedure, for instance, to repeatedly generate GF(42) out of the previously created GF(4) by using D2+D+2 over GF(4). Table 3 shows the corresponding elements of GF(16).
Listing Two
/* constructs table of operations <po> of locator field <pg> given symbol field
** operations <op>. Returns <po> Since 2D operations commute, only 1/2 of the
** NxN tables are necessary. Proper indexing into 2D-arrays is done using
** ADD/MULT macros. */
#define I(i,j) ((i>j) ? j : i)
#define J(i,j) ((i>j) ? (i-j) : (j-i))
#define MULT(i,j) Mult[I(i,j)][J(i,j)]
#define ADD(i,j) Add[I(i,j)][J(i,j)]
pOPTABLE ConstructOpTable(pOPTABLE po,pGALOIS pg,pOPTABLE op) {
unsigned i,j;
for(i=0;i<po->MaxNo;i++) {
po->InvAdd[i] = (byte) i_add_inv(i,pg,op);
po->InvMult[i] = (byte) i_mult_inv(i,pg,op);
for(j=i;j<po->MaxNo;j++) {
po->Add[i][j-i] = (byte) i_add(i,j,pg,op);
po->Mult[i][j-i]= (byte) i_mult(i,j,pg,op);
}
}
return po;
}
/* Note: All references to GF-elements denote indexes to the <pg>-structure
** containing these elements. The convention is:
** index = 0 -> GF-element = 0
** 1 -> 1
** 2 -> `
** 3 -> `^2
** ...
** p^q - 1 -> `^(p^q-2)
*/
/* `^k = `^i * `^j */
int i_mult(int i,int j,pGALOIS pg,pOPTABLE po) {
int res;
int mod=pg->MaxNo;
if(po==NULL) /* symbol field created -> modulo */
return( (i*j)%pg->Base);
if(i==0 || j==0) /* 0 */
return 0;
res=i+j-1;
while(res>=mod)
res -= (mod-1);
return res;
}
/* `^k = `^i + `^j */
int i_add(int i,int j,pGALOIS pg,pOPTABLE po) {
byte res[8];
int k,l;
if(po==NULL) /* symbol field created -> modulo */
return( (i+j)%pg->Base);
for(k=pg->Exp - 1;k>=0;k--)
res[k]=po->ADD(pg->pe[i].Alpha[k],pg->pe[j].Alpha[k]);
for(l=0;;l++) {
for(k=0;k<pg->Exp;k++) {
if( res[k] != pg->pe[l].Alpha[k] )
break;
}
if(k==pg->Exp)
break;
}
return l;
}
/* `^k = (`^i)^n */
int i_pow(int i,int n,pGALOIS pg,pOPTABLE po) {
int res=1;
while (n--)
res=i_mult(res,i,pg,po);
return res;
}
/* `^k = - (`^i) */
int i_add_inv(int i,pGALOIS pg,pOPTABLE po) {
int j=0;
while( 0 != i_add(i,j,pg,po) ) j++;
return j;
}
/* `^k = (`^i)^(-1) */
int i_mult_inv(int i,pGALOIS pg,pOPTABLE po) {
int j=0;
if(!i) return 0; /* no multiplicative inverse for 0 */
while( 1 != i_mult(i,j,pg,po) ) j++;
return j;
}
Table 2
Table 3
Generating BCH/RS codes
BCH codes are designed for powerful error recovery--if you need more robustness than Hamming provides, BCH is the way to go. The first step in generating a specific BCH or RS code is the choice and construction of both GF(p) and GF(pn) fields. GF(p) contains the symbols and defines operations on these symbols. In the case of GF(2), a symbol is just a bit; if p is equal to 16, the code deals with 16-ary symbols, each one represented by four bits. The block length N of the code is dependent on your choice of n in GF(pn) since N=pn-1 for BCH and RS codes. So if you choose p=16 and n=2, the coder will generate blocks of 255 symbols, each symbol consisting of four bits. Just one parameter remains to be determined: the number of symbol errors (not bit errors) per block N that can be recovered; this is denoted by t. Up to four bit errors cause only one symbol error to occur if these four erroneous bits all occur within the same symbol boundary.
Coding theory proves that RS codes maximize error correctability for a set of possible code words having a specific minimum distance (the minimum Hamming distance between two valid N-length code words). The family of BCH codes approximates this theoretical limit. The generators for BCH and RS codes are given in Figure 2. RS codes actually constitute a subset of BCH codes: If both symbol and locator field are the same (say, GF(p) and GF(p1), respectively), the code generated is RS. In this case, all locator-field elements also belong to the symbol field. The minimal polynomial of a locator-field element ai then simply reduces to (D-i), yielding the RS canonical form. Note also that the number of recoverable symbol errors per block t always equals one-half the number of parity symbols (N-L) added to the data block, the latter indicated by the generator's degree. Listing Three shows how GALOIS.EXE constructs the minimal polynomial for GF element i and how to get to the generator polynomial for a specific BCH or RS code.
Figure 2
Listing Three
/* Gets the minimal polynomial for each field element (except for 0 and 1),
** given locator field's elements and operations. */
void FindMinimalPolynomials(void) {
int i,j,k,order;
pPOLYNOM pp_min;
...(allocate pp_aux and init it to D) ...
for(i=2;i<pg_l->MaxNo;i++) {
/* We first get the order of the element since this */
/* will determine the degree of the minimal polynomial */
for(order=1;
i_pow(i,(int)pow((double)pg_l->Base,(double)order),pg_l,po_l) != i;
order++) ;
...(allocate pp_min to this degree and init to 1)...
/* Now we'll actually start constructing */
/* it : m`i(D)=(D-`i)(D-`i^2)...(D-`i^q) */
/* where q equals the order of `i. */
for(j=0;j<order;j++) {
pp_aux->c[1]=1;
pp_aux->c[0]=
po_l->InvAdd[
i_pow(i,(int)pow((double)pg_l->Base,(double)j),pg_l,po_l)];
pp_min=MultPoly(pp_min,pp_aux,pp_min,po_l);
}
/* coeff. of <pp_min> n SYMBOL field, but */
/* operations are done in the locator */
/* field, so remap these values */
for(k=order;k>=0;k--)
pp_min->c[k]= (int)pg_l->pe[ pp_min->c[k] ].Alpha[0];
DEGREE(pp_min);
}
...(free pp_aux)...
}
/* Generates BCH/RS generator polynomial in symbol and locator field set up
** previously. */
extern int t; /* error correctability */
extern pPOLYNOM pp_gen; /* generator polynomial */
extern bool ComparePoly(pPOLYNOM p1,pPOLYNOM p2);
/* returns true if p1==p2 */
void GenerateBCH(void) {
bool common_factor;
int deg,i,j;
...(init pp_gen to 1) ...
deg=1+(t<<1); /* 2t+1 */
for(i=2;i<=deg;i++) {
common_factor=FALSE;
for(j=2;j<i;j++) {
if(ComparePoly(pg_l->pe[j].pp,pg_l->pe[i].pp)) {
common_factor=TRUE;
break;
}
}
if(!common_factor)
pp_gen=MultPoly(pp_gen,pg_l->pe[i].pp,pp_gen,po_s);
}
DEGREE(pp_gen);
}
END.
Decoding BCH/RS Codes
How can a receiver know about the occurrence of one or more errors during transmission or storage of the encoded block, given only the generator polynomial g(D) and the received code word y, represented by its polynomial y(D)? Assuming that code word x(D) is transmitted, you can write y(D)=x(D)+e(D), introducing the error polynomial e(D). Remember, x(D) equals u(D).g(D) and g(i)=0 (i=1..2t) for BCH/RS codes by construction of the generator. Hence, calculating y() when receiving the block yields a nonzero result only when errors have occurred. When their nonzero status points to an erroneous condition, these 2it values y() are called "syndromes" and are denoted S(i).
i
Now you know, on a per block basis, whether or not errors have occurred, but you don't know how many symbols are affected, nor the error positions. You can find these positions by constructing yet another polynomial: the "error locator" ![[sigma]](http://twimgs.com/ddj/ddj/images/ddj9412d/sigma12.gif)
(D). The inverses of its roots will indicate the positions in error. Figure 3 defines (D).
Figure 3
The key to decoding BCH/RS codes is the construction of this polynomial using the Massey-Berlekamp algorithm, which is generally more applicable for synthesizing a linear-shift feedback register (LSFR) of minimum length and generating a predefined output sequence. Figure 4 shows the algorithm. Table 4 illustrates the procedure in terms of LSFR generation for decoding a triple-error-correcting binary (15,5) BCH-encoded block with three positions in error. Note that a LSFR generating the 2t syndromes is constructed in this way. Each time the register fails to generate a syndrome, its length is incremented. When the algorithm ends, the LSFR will generate the 2t syndromes correctly if the number of errors in the block is less than or equal to t. As you see from the example, the LSFR is readily described by the polynomial (D). Solve (D) for its roots, and invert these to obtain the error locators.
Figure 4
Table 4
To decode binary BCH codes where the symbol field is GF(2), you only have to know the position of symbol errors. In nonbinary codes, however, you also need to know the error magnitudes. You can get around by modifying the error locator (to incorporate error values, not only error positions) and by defining an extra error evaluator polynomial. GALOIS.EXE incorporates the code for doing this.
Conclusion
While these error-correcting codes have been known for some time, their application in high-speed digital transmission has only recently come into vogue because of the complex decoding involved. Chips performing RS encoding/decoding using 256-ary symbols are now commercially available at bit rates of up to 40 Mbits/sec and higher. There's little doubt that they would be a useful building block for the construction of error-resilient, safe file systems.
References
Anderson, J.B. and S. Mohan. Source and Channel Coding: An Algorithmic Approach. Boston, MA: Kluwer Academic Publishers, 1991.
J.L. Massey. "Shift Register Synthesis and BCH Decoding." IEEE Transactions on Information Theory (January 1969).
Bart is project manager for digital video communications in the Broadcast & Cable division of Barco nv., Kortrijk, Belgium. He can be reached at [email protected].
