diff --git a/src/bch.c b/src/bch.c new file mode 100644 index 0000000..415cdb5 --- /dev/null +++ b/src/bch.c @@ -0,0 +1,589 @@ +/* + * File: bch3.c + * Title: Encoder/decoder for binary BCH codes in C (Version 3.1) + * Author: Robert Morelos-Zaragoza + * Date: August 1994 + * Revised: June 13, 1997 + * + * =============== Encoder/Decoder for binary BCH codes in C ================= + * + * Version 1: Original program. The user provides the generator polynomial + * of the code (cumbersome!). + * Version 2: Computes the generator polynomial of the code. + * Version 3: No need to input the coefficients of a primitive polynomial of + * degree m, used to construct the Galois Field GF(2**m). The + * program now works for any binary BCH code of length such that: + * 2**(m-1) - 1 < length <= 2**m - 1 + * + * Note: You may have to change the size of the arrays to make it work. + * + * The encoding and decoding methods used in this program are based on the + * book "Error Control Coding: Fundamentals and Applications", by Lin and + * Costello, Prentice Hall, 1983. + * + * Thanks to Patrick Boyle (pboyle@era.com) for his observation that 'bch2.c' + * did not work for lengths other than 2**m-1 which led to this new version. + * Portions of this program are from 'rs.c', a Reed-Solomon encoder/decoder + * in C, written by Simon Rockliff (simon@augean.ua.oz.au) on 21/9/89. The + * previous version of the BCH encoder/decoder in C, 'bch2.c', was written by + * Robert Morelos-Zaragoza (robert@spectra.eng.hawaii.edu) on 5/19/92. + * + * NOTE: + * The author is not responsible for any malfunctioning of + * this program, nor for any damage caused by it. Please include the + * original program along with these comments in any redistribution. + * + * For more information, suggestions, or other ideas on implementing error + * correcting codes, please contact me at: + * + * Robert Morelos-Zaragoza + * 5120 Woodway, Suite 7036 + * Houston, Texas 77056 + * + * email: r.morelos-zaragoza@ieee.org + * + * COPYRIGHT NOTICE: This computer program is free for non-commercial purposes. + * You may implement this program for any non-commercial application. You may + * also implement this program for commercial purposes, provided that you + * obtain my written permission. Any modification of this program is covered + * by this copyright. + * + * == Copyright (c) 1994-7, Robert Morelos-Zaragoza. All rights reserved. == + * + * 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) + * + */ + +#include +#include +#include +#include +#include "bch.h" + +int bch_init(bch_t *bch, int m, int length, int t) { + + int p[21], n; + + if (bch == NULL) { + return -1; + } + + if (m < 2 || m > 20) { + return -2; + } + + bch->m = m; + bch->length = length; + bch->t = t; + + for (int i=1; in = n; + int ninf = (n + 1) / 2 - 1; + + if (length < ninf || length > n) { + return -3; + } + +/* + * 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; + + bch->alpha_to = malloc(n * sizeof(int)); + bch->index_of = malloc(n * sizeof(int)); + + mask = 1; + bch->alpha_to[m] = 0; + for (int i = 0; i < m; i++) { + bch->alpha_to[i] = mask; + bch->index_of[bch->alpha_to[i]] = i; + if (p[i] != 0) + bch->alpha_to[m] ^= mask; + mask <<= 1; + } + bch->index_of[bch->alpha_to[m]] = m; + mask >>= 1; + for (int i = m + 1; i < n; i++) { + if (bch->alpha_to[i - 1] >= mask) + bch->alpha_to[i] = bch->alpha_to[m] ^ ((bch->alpha_to[i - 1] ^ mask) << 1); + else + bch->alpha_to[i] = bch->alpha_to[i - 1] << 1; + bch->index_of[bch->alpha_to[i]] = i; + } + bch->index_of[0] = -1; + +/* + * 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 (bch->m > 9) { + printf("Computing cycle sets modulo %d\n", bch->n); + printf("(This may take some time)...\n"); + } + do { + /* Generate the jj-th cycle set */ + ii = 0; + do { + ii++; + cycle[jj][ii] = (cycle[jj][ii - 1] * 2) % bch->n; + size[jj]++; + aux = (cycle[jj][ii] * 2) % bch->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 < (bch->n - 1))); + if (!(test)) { + jj++; /* next cycle set index */ + cycle[jj][0] = ll; + size[jj] = 1; + } + } while (ll < (bch->n - 1)); + nocycles = jj; /* number of cycle sets modulo n */ + + int 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++; + } + + bch-> k = length - rdncy; + + if (bch->k<0) + { + printf("Parameters invalid!\n"); + return -4; + } + + printf("This is a (%d, %d, %d) binary BCH code\n", bch->length, bch->k, d); + + /* Compute the generator polynomial */ + bch->g = malloc(rdncy * sizeof(int)); + bch->g[0] = bch->alpha_to[zeros[1]]; + bch->g[1] = 1; /* g(x) = (X + zeros[1]) initially */ + for (ii = 2; ii <= rdncy; ii++) { + bch->g[ii] = 1; + for (jj = ii - 1; jj > 0; jj--) + if (bch->g[jj] != 0) + bch->g[jj] = bch->g[jj - 1] ^ bch->alpha_to[(bch->index_of[bch->g[jj]] + zeros[ii]) % bch->n]; + else + bch->g[jj] = bch->g[jj - 1]; + bch->g[0] = bch->alpha_to[(bch->index_of[bch->g[0]] + zeros[ii]) % bch->n]; + } + printf("Generator polynomial:\ng(x) = "); + for (ii = 0; ii <= rdncy; ii++) { + printf("%d", bch->g[ii]); + } + printf("\n"); + + return 0; +} + +void generate_bch(bch_t *bch, int *data, int *bb) { +/* + * 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 feedback; + + for (int i = 0; i < bch->length - bch->k; i++) + bb[i] = 0; + for (int i = bch->k - 1; i >= 0; i--) { + feedback = data[i] ^ bb[bch->length - bch->k - 1]; + if (feedback != 0) { + for (int j = bch->length - bch->k - 1; j > 0; j--) + if (bch->g[j] != 0) + bb[j] = bb[j - 1] ^ feedback; + else + bb[j] = bb[j - 1]; + bb[0] = bch->g[0] && feedback; + } else { + for (int j = bch->length - bch->k - 1; j > 0; j--) + bb[j] = bb[j - 1]; + bb[0] = 0; + } + } +} + + +int +apply_bch(bch_t *bch, int *recd) +/* + * 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]; + + t2 = 2 * bch->t; + + /* first form the syndromes */ + printf("S(x) = "); + for (i = 1; i <= t2; i++) { + s[i] = 0; + for (j = 0; j < bch->length; j++) + if (recd[j] != 0) + s[i] ^= bch->alpha_to[(i * j) % bch->n]; + if (s[i] != 0) + syn_error = 1; /* set error flag if non-zero syndrome */ +/* + * 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] = bch->index_of[s[i]]; + printf("%3d ", s[i]); + } + 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] = bch->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] = + bch->alpha_to[(d[u] + bch->n - d[q] + elp[q][i]) % bch->n]; + for (i = 0; i <= l[u]; i++) { + elp[u + 1][i] ^= elp[u][i]; + elp[u][i] = bch->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] = bch->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] ^= bch->alpha_to[(s[u + 1 - i] + + bch->index_of[elp[u + 1][i]]) % bch->n]; + /* put d[u+1] into index form */ + d[u + 1] = bch->index_of[d[u + 1]]; + } + } while ((u < t2) && (l[u + 1] <= bch->t)); + + u++; + if (l[u] <= bch->t) {/* Can correct errors */ + /* put elp into index form */ + for (i = 0; i <= l[u]; i++) + elp[u][i] = bch->index_of[elp[u][i]]; + + printf("sigma(x) = "); + for (i = 0; i <= l[u]; i++) + printf("%3d ", elp[u][i]); + printf("\n"); + printf("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 <= bch->n; i++) { + q = 1; + for (j = 1; j <= l[u]; j++) + if (reg[j] != -1) { + reg[j] = (reg[j] + j) % bch->n; + q ^= bch->alpha_to[reg[j]]; + } + if (!q) { /* store root and error + * location number indices */ + root[count] = i; + loc[count] = bch->n - i; + count++; + printf("%3d ", bch->n - i); + } + } + 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; + return l[u]; + } + else { /* elp has degree >t hence cannot solve */ + printf("Incomplete decoding: errors detected\n"); + return -1; + } + } else { + return -1; + } + } else { + return 0; // No errors + } +} + +/* LEFT justified in hex */ +void bytes_to_bits(int *bytes, int *bit_dest, int num_bits) { + for (int i = 0; i < num_bits; i++) { + int index = i / 8; + int bit_pos = 7 - (i % 8); + int bit_mask = 1 << bit_pos; + bit_dest[i] = (bytes[index] & bit_mask) != 0; + } +} + +void dump_bch(bch_t *bch) { + printf("m: %d length: %d t: %d n: %d k: %d\n", bch->m, bch->length, bch->t, bch->n, bch->k); +} + +#define TEST_APPLY + +int main() +{ + int test[][8] = { + { 0xb2, 0x17, 0xa2, 0xb9, 0x53, 0xdd, 0xc5, 0x52 }, /* perfect random test */ + { 0xf0, 0x5a, 0x6a, 0x6a, 0x01, 0x63, 0x33, 0xd0 }, /* g001-cut-lenthened_457.938M.wav */ + { 0xf0, 0x81, 0x52, 0x6b, 0x71, 0xa5, 0x63, 0x08 }, /* 1st in eotd_received_data */ +/* 3 errors */ { 0xf0, 0x85, 0x50, 0x6a, 0x01, 0xe5, 0x6e, 0x84 }, /* 2nd in eotd_received_data */ +/* fixed */ { 0xf0, 0x85, 0x50, 0x6a, 0x01, 0xe5, 0x06, 0x84 }, /* 2nd, but with the bits fixed */ + { 0xf0, 0x85, 0x59, 0x5a, 0x01, 0xe5, 0x6e, 0x84 }, /* 3rd */ + { 0xf1, 0x34, 0x50, 0x1a, 0x01, 0xe5, 0x66, 0xfe }, /* 4th */ + { 0xf0, 0xeb, 0x10, 0xea, 0x01, 0x6e, 0x54, 0x1c }, /* 5th */ + { 0xf0, 0xea, 0x5c, 0xea, 0x01, 0x6e, 0x55, 0x0e }, /* 6th */ + { 0xe0, 0x21, 0x10, 0x1a, 0x01, 0x32, 0xbc, 0xe4 }, /* Sun Mar 20 05:41:00 2022 */ + { 0xf0, 0x42, 0x50, 0x5b, 0xcf, 0xd5, 0x64, 0xe4 }, /* Sun Mar 20 12:58:43 2022 */ + { 0xf0, 0x8c, 0x10, 0xaa, 0x01, 0x73, 0x7b, 0x1a }, /* Sun Mar 20 13:35:48 2022 */ + { 0xf0, 0x8c, 0x10, 0xb1, 0xc0, 0xe0, 0x90, 0x64 }, /* Sun Mar 20 13:37:05 2022 */ + { 0xf0, 0x8c, 0x10, 0x6a, 0x01, 0x64, 0x7a, 0xe8 }, /* Sun Mar 20 13:37:48 2022 */ + { 0x50, 0x8c, 0x12, 0x6a, 0x01, 0x64, 0x7a, 0xe8 }, + }; + + int bits[63]; + bch_t bch; + + bch_init(&bch, 6, 63, 3); + + for (int count = 0; count < sizeof(test) / sizeof(*test); count++) { + bytes_to_bits(test[count], bits, 63); +#ifdef TEST_BYTES_TO_BITS + printf("ORIG pkt [%d]\n", count); + for (int i = 0; i < 8; i++) { + printf("%02x ", test[count][i]); + } + printf("\n"); + + printf("ORIG pkt[%d] bits\n", count); + for (int i = 0; i < 63; i++) { + printf("%d ", bits[i]); + } + printf("\n"); +#endif +#ifdef TEST_GENERATE + int bch_code[18]; + generate_bch(&bch, bits, bch_code); + printf("generated BCH\n"); + for (int i = 0; i < 18; i++) { + printf("%d ", bch_code[i]); + } + printf("\n"); +#endif +#ifdef TEST_APPLY + int recv[63]; + // backwards, for now + for (int i = 0; i < 45; i++) { + recv[i + 18] = bits[i]; + } + + for (int i = 0; i < 18; i++) { + recv[i] = bits[i + 45]; + } + + printf("rearranged packet: "); + for (int i = 0; i < 63; i++) { + printf("%d ", recv[i]); + } + printf("\n"); + + int corrected = apply_bch(&bch, recv); + + printf("corrected [%d] packet: ", corrected); + for (int i = 0; i < 63; i++) { + printf("%d ", recv[i]); + } + printf("\n"); +#endif + } +} diff --git a/src/bch.h b/src/bch.h new file mode 100644 index 0000000..1f90da5 --- /dev/null +++ b/src/bch.h @@ -0,0 +1,22 @@ +#ifndef __BCH_H +#define __BCH_H + +struct bch { + int m; // 2^m - 1 is max length, n + int length; // Actual packet size + int n; // 2^m - 1 + int k; // Length of data portion + int t; // Number of correctable bits + + int *g; // Calculated polynomial of length n - k + int *alpha_to; + int *index_of; +}; + +typedef struct bch bch_t; + +int bch_init(bch_t *bch, int m, int length, int t); + +void generate_bch(bch_t *bch, int *data, int *bb); + +#endif