|
| 1 | +/* |
| 2 | + * BCHCode.c |
| 3 | + * |
| 4 | + * Copyright (C) 2015 Craig Shelley (craig@microtron.org.uk) |
| 5 | + * |
| 6 | + * BCH Encoder/Decoder - Adapted from GNURadio for use with Multimon |
| 7 | + * |
| 8 | + * This program is free software; you can redistribute it and/or modify |
| 9 | + * it under the terms of the GNU General Public License as published by |
| 10 | + * the Free Software Foundation; either version 2 of the License, or |
| 11 | + * (at your option) any later version. |
| 12 | + * |
| 13 | + * This program is distributed in the hope that it will be useful, |
| 14 | + * but WITHOUT ANY WARRANTY; without even the implied warranty of |
| 15 | + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
| 16 | + * GNU General Public License for more details. |
| 17 | + * |
| 18 | + * You should have received a copy of the GNU General Public License |
| 19 | + * along with this program; if not, write to the Free Software |
| 20 | + * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. |
| 21 | + */ |
| 22 | + |
| 23 | +#include <math.h> |
| 24 | +#include <stdlib.h> |
| 25 | +#include "BCHCode.h" |
| 26 | + |
| 27 | +struct BCHCode { |
| 28 | + int * p; // coefficients of primitive polynomial used to generate GF(2**5) |
| 29 | + int m; // order of the field GF(2**5) = 5 |
| 30 | + int n; // 2**5 - 1 = 31 |
| 31 | + int k; // n - deg(g(x)) = 21 = dimension |
| 32 | + int t; // 2 = error correcting capability |
| 33 | + int * alpha_to; // log table of GF(2**5) |
| 34 | + int * index_of; // antilog table of GF(2**5) |
| 35 | + int * g; // coefficients of generator polynomial, g(x) [n - k + 1]=[11] |
| 36 | + int * bb; // coefficients of redundancy polynomial ( x**(10) i(x) ) modulo g(x) |
| 37 | +}; |
| 38 | + |
| 39 | + |
| 40 | + |
| 41 | + |
| 42 | + |
| 43 | + |
| 44 | +static void generate_gf(struct BCHCode * BCHCode_data) { |
| 45 | + if (BCHCode_data==NULL) return; |
| 46 | + /* |
| 47 | + * generate GF(2**m) from the irreducible polynomial p(X) in p[0]..p[m] |
| 48 | + * lookup tables: index->polynomial form alpha_to[] contains j=alpha**i; |
| 49 | + * polynomial form -> index form index_of[j=alpha**i] = i alpha=2 is the |
| 50 | + * primitive element of GF(2**m) |
| 51 | + */ |
| 52 | + |
| 53 | + register int i, mask; |
| 54 | + mask = 1; |
| 55 | + BCHCode_data->alpha_to[BCHCode_data->m] = 0; |
| 56 | + for (i = 0; i < BCHCode_data->m; i++) { |
| 57 | + BCHCode_data->alpha_to[i] = mask; |
| 58 | + BCHCode_data->index_of[BCHCode_data->alpha_to[i]] = i; |
| 59 | + if (BCHCode_data->p[i] != 0) |
| 60 | + BCHCode_data->alpha_to[BCHCode_data->m] ^= mask; |
| 61 | + mask <<= 1; |
| 62 | + } |
| 63 | + BCHCode_data->index_of[BCHCode_data->alpha_to[BCHCode_data->m]] = BCHCode_data->m; |
| 64 | + mask >>= 1; |
| 65 | + for (i = BCHCode_data->m + 1; i < BCHCode_data->n; i++) { |
| 66 | + if (BCHCode_data->alpha_to[i - 1] >= mask) |
| 67 | + BCHCode_data->alpha_to[i] = BCHCode_data->alpha_to[BCHCode_data->m] ^ ((BCHCode_data->alpha_to[i - 1] ^ mask) << 1); |
| 68 | + else |
| 69 | + BCHCode_data->alpha_to[i] = BCHCode_data->alpha_to[i - 1] << 1; |
| 70 | + BCHCode_data->index_of[BCHCode_data->alpha_to[i]] = i; |
| 71 | + } |
| 72 | + BCHCode_data->index_of[0] = -1; |
| 73 | +} |
| 74 | + |
| 75 | + |
| 76 | +static void gen_poly(struct BCHCode * BCHCode_data) { |
| 77 | + if (BCHCode_data==NULL) return; |
| 78 | + /* |
| 79 | + * Compute generator polynomial of BCH code of length = 31, redundancy = 10 |
| 80 | + * (OK, this is not very efficient, but we only do it once, right? :) |
| 81 | + */ |
| 82 | + |
| 83 | + register int ii, jj, ll, kaux; |
| 84 | + int test, aux, nocycles, root, noterms, rdncy; |
| 85 | + int cycle[15][6], size[15], min[11], zeros[11]; |
| 86 | + /* Generate cycle sets modulo 31 */ |
| 87 | + cycle[0][0] = 0; size[0] = 1; |
| 88 | + cycle[1][0] = 1; size[1] = 1; |
| 89 | + jj = 1; /* cycle set index */ |
| 90 | + do { |
| 91 | + /* Generate the jj-th cycle set */ |
| 92 | + ii = 0; |
| 93 | + do { |
| 94 | + ii++; |
| 95 | + cycle[jj][ii] = (cycle[jj][ii - 1] * 2) % BCHCode_data->n; |
| 96 | + size[jj]++; |
| 97 | + aux = (cycle[jj][ii] * 2) % BCHCode_data->n; |
| 98 | + } while (aux != cycle[jj][0]); |
| 99 | + /* Next cycle set representative */ |
| 100 | + ll = 0; |
| 101 | + do { |
| 102 | + ll++; |
| 103 | + test = 0; |
| 104 | + for (ii = 1; ((ii <= jj) && (!test)); ii++) { |
| 105 | + /* Examine previous cycle sets */ |
| 106 | + for (kaux = 0; ((kaux < size[ii]) && (!test)); kaux++) { |
| 107 | + if (ll == cycle[ii][kaux]) { |
| 108 | + test = 1; |
| 109 | + } |
| 110 | + } |
| 111 | + } |
| 112 | + } while ((test) && (ll < (BCHCode_data->n - 1))); |
| 113 | + if (!(test)) { |
| 114 | + jj++; /* next cycle set index */ |
| 115 | + cycle[jj][0] = ll; |
| 116 | + size[jj] = 1; |
| 117 | + } |
| 118 | + } while (ll < (BCHCode_data->n - 1)); |
| 119 | + nocycles = jj; /* number of cycle sets modulo BCHCode_data->n */ |
| 120 | + /* Search for roots 1, 2, ..., BCHCode_data->d-1 in cycle sets */ |
| 121 | + kaux = 0; |
| 122 | + rdncy = 0; |
| 123 | + for (ii = 1; ii <= nocycles; ii++) { |
| 124 | + min[kaux] = 0; |
| 125 | + for (jj = 0; jj < size[ii]; jj++) { |
| 126 | + for (root = 1; root < (2*BCHCode_data->t + 1); root++) { |
| 127 | + if (root == cycle[ii][jj]) { |
| 128 | + min[kaux] = ii; |
| 129 | + } |
| 130 | + } |
| 131 | + } |
| 132 | + if (min[kaux]) { |
| 133 | + rdncy += size[min[kaux]]; |
| 134 | + kaux++; |
| 135 | + } |
| 136 | + } |
| 137 | + noterms = kaux; |
| 138 | + kaux = 1; |
| 139 | + for (ii = 0; ii < noterms; ii++) { |
| 140 | + for (jj = 0; jj < size[min[ii]]; jj++) { |
| 141 | + zeros[kaux] = cycle[min[ii]][jj]; |
| 142 | + kaux++; |
| 143 | + } |
| 144 | + } |
| 145 | + //printf("This is a (%d, %d, %d) binary BCH code\n", BCHCode_data->n, BCHCode_data->k, BCHCode_data->d); |
| 146 | + /* Compute generator polynomial */ |
| 147 | + BCHCode_data->g[0] = BCHCode_data->alpha_to[zeros[1]]; |
| 148 | + BCHCode_data->g[1] = 1; /* g(x) = (X + zeros[1]) initially */ |
| 149 | + for (ii = 2; ii <= rdncy; ii++) { |
| 150 | + BCHCode_data->g[ii] = 1; |
| 151 | + for (jj = ii - 1; jj > 0; jj--) { |
| 152 | + if (BCHCode_data->g[jj] != 0) |
| 153 | + BCHCode_data->g[jj] = BCHCode_data->g[jj - 1] ^ BCHCode_data->alpha_to[(BCHCode_data->index_of[BCHCode_data->g[jj]] + zeros[ii]) % BCHCode_data->n]; |
| 154 | + else |
| 155 | + BCHCode_data->g[jj] = BCHCode_data->g[jj - 1]; |
| 156 | + } |
| 157 | + BCHCode_data->g[0] = BCHCode_data->alpha_to[(BCHCode_data->index_of[BCHCode_data->g[0]] + zeros[ii]) % BCHCode_data->n]; |
| 158 | + } |
| 159 | + //printf("g(x) = "); |
| 160 | + //for (ii = 0; ii <= rdncy; ii++) { |
| 161 | + // printf("%d", BCHCode_data->g[ii]); |
| 162 | + // if (ii && ((ii % 70) == 0)) { |
| 163 | + // printf("\n"); |
| 164 | + // } |
| 165 | + //} |
| 166 | + //printf("\n"); |
| 167 | +} |
| 168 | + |
| 169 | + |
| 170 | +void BCHCode_Encode(struct BCHCode * BCHCode_data, int data[]) { |
| 171 | + if (BCHCode_data==NULL) return; |
| 172 | + /* |
| 173 | + * Calculate redundant bits bb[], codeword is c(X) = data(X)*X**(n-k)+ bb(X) |
| 174 | + */ |
| 175 | + |
| 176 | + register int i, j; |
| 177 | + register int feedback; |
| 178 | + for (i = 0; i < BCHCode_data->n - BCHCode_data->k; i++) { |
| 179 | + BCHCode_data->bb[i] = 0; |
| 180 | + } |
| 181 | + for (i = BCHCode_data->k - 1; i >= 0; i--) { |
| 182 | + feedback = data[i] ^ BCHCode_data->bb[BCHCode_data->n - BCHCode_data->k - 1]; |
| 183 | + if (feedback != 0) { |
| 184 | + for (j = BCHCode_data->n - BCHCode_data->k - 1; j > 0; j--) { |
| 185 | + if (BCHCode_data->g[j] != 0) { |
| 186 | + BCHCode_data->bb[j] = BCHCode_data->bb[j - 1] ^ feedback; |
| 187 | + } else { |
| 188 | + BCHCode_data->bb[j] = BCHCode_data->bb[j - 1]; |
| 189 | + } |
| 190 | + } |
| 191 | + BCHCode_data->bb[0] = BCHCode_data->g[0] && feedback; |
| 192 | + } else { |
| 193 | + for (j = BCHCode_data->n - BCHCode_data->k - 1; j > 0; j--) { |
| 194 | + BCHCode_data->bb[j] = BCHCode_data->bb[j - 1]; |
| 195 | + } |
| 196 | + BCHCode_data->bb[0] = 0; |
| 197 | + }; |
| 198 | + }; |
| 199 | +}; |
| 200 | + |
| 201 | + |
| 202 | +int BCHCode_Decode(struct BCHCode * BCHCode_data, int recd[]) { |
| 203 | + if (BCHCode_data==NULL) return -1; |
| 204 | + /* |
| 205 | + * We do not need the Berlekamp algorithm to decode. |
| 206 | + * We solve before hand two equations in two variables. |
| 207 | + */ |
| 208 | + |
| 209 | + register int i, j, q; |
| 210 | + int elp[3], s[5], s3; |
| 211 | + int count = 0, syn_error = 0; |
| 212 | + int loc[3], reg[3]; |
| 213 | + int aux; |
| 214 | + int retval=0; |
| 215 | + /* first form the syndromes */ |
| 216 | + // printf("s[] = ("); |
| 217 | + for (i = 1; i <= 4; i++) { |
| 218 | + s[i] = 0; |
| 219 | + for (j = 0; j < BCHCode_data->n; j++) { |
| 220 | + if (recd[j] != 0) { |
| 221 | + s[i] ^= BCHCode_data->alpha_to[(i * j) % BCHCode_data->n]; |
| 222 | + } |
| 223 | + } |
| 224 | + if (s[i] != 0) { |
| 225 | + syn_error = 1; /* set flag if non-zero syndrome */ |
| 226 | + } |
| 227 | + /* NOTE: If only error detection is needed, |
| 228 | + * then exit the program here... |
| 229 | + */ |
| 230 | + /* convert syndrome from polynomial form to index form */ |
| 231 | + s[i] = BCHCode_data->index_of[s[i]]; |
| 232 | + //printf("%3d ", s[i]); |
| 233 | + }; |
| 234 | + //printf(")\n"); |
| 235 | + if (syn_error) { /* If there are errors, try to correct them */ |
| 236 | + if (s[1] != -1) { |
| 237 | + s3 = (s[1] * 3) % BCHCode_data->n; |
| 238 | + if ( s[3] == s3 ) { /* Was it a single error ? */ |
| 239 | + //printf("One error at %d\n", s[1]); |
| 240 | + recd[s[1]] ^= 1; /* Yes: Correct it */ |
| 241 | + } else { |
| 242 | + /* Assume two errors occurred and solve |
| 243 | + * for the coefficients of sigma(x), the |
| 244 | + * error locator polynomail |
| 245 | + */ |
| 246 | + if (s[3] != -1) { |
| 247 | + aux = BCHCode_data->alpha_to[s3] ^ BCHCode_data->alpha_to[s[3]]; |
| 248 | + } else { |
| 249 | + aux = BCHCode_data->alpha_to[s3]; |
| 250 | + } |
| 251 | + elp[0] = 0; |
| 252 | + elp[1] = (s[2] - BCHCode_data->index_of[aux] + BCHCode_data->n) % BCHCode_data->n; |
| 253 | + elp[2] = (s[1] - BCHCode_data->index_of[aux] + BCHCode_data->n) % BCHCode_data->n; |
| 254 | + //printf("sigma(x) = "); |
| 255 | + //for (i = 0; i <= 2; i++) { |
| 256 | + // printf("%3d ", elp[i]); |
| 257 | + //} |
| 258 | + //printf("\n"); |
| 259 | + //printf("Roots: "); |
| 260 | + /* find roots of the error location polynomial */ |
| 261 | + for (i = 1; i <= 2; i++) { |
| 262 | + reg[i] = elp[i]; |
| 263 | + } |
| 264 | + count = 0; |
| 265 | + for (i = 1; i <= BCHCode_data->n; i++) { /* Chien search */ |
| 266 | + q = 1; |
| 267 | + for (j = 1; j <= 2; j++) { |
| 268 | + if (reg[j] != -1) { |
| 269 | + reg[j] = (reg[j] + j) % BCHCode_data->n; |
| 270 | + q ^= BCHCode_data->alpha_to[reg[j]]; |
| 271 | + } |
| 272 | + } |
| 273 | + if (!q) { /* store error location number indices */ |
| 274 | + loc[count] = i % BCHCode_data->n; |
| 275 | + count++; |
| 276 | + //printf("%3d ", (i%n)); |
| 277 | + } |
| 278 | + } |
| 279 | + //printf("\n"); |
| 280 | + if (count == 2) { |
| 281 | + /* no. roots = degree of elp hence 2 errors */ |
| 282 | + for (i = 0; i < 2; i++) |
| 283 | + recd[loc[i]] ^= 1; |
| 284 | + } else { /* Cannot solve: Error detection */ |
| 285 | + retval=1; |
| 286 | + //for (i = 0; i < 31; i++) { |
| 287 | + // recd[i] = 0; |
| 288 | + //} |
| 289 | + //printf("incomplete decoding\n"); |
| 290 | + } |
| 291 | + } |
| 292 | + } else if (s[2] != -1) {/* Error detection */ |
| 293 | + retval=1; |
| 294 | + //for (i = 0; i < 31; i++) recd[i] = 0; |
| 295 | + //printf("incomplete decoding\n"); |
| 296 | + } |
| 297 | + } |
| 298 | + |
| 299 | + return retval; |
| 300 | +} |
| 301 | + |
| 302 | +/* |
| 303 | + * Example usage BCH(31,21,5) |
| 304 | + * |
| 305 | + * p[] = coefficients of primitive polynomial used to generate GF(2**5) |
| 306 | + * m = order of the field GF(2**5) = 5 |
| 307 | + * n = 2**5 - 1 = 31 |
| 308 | + * t = 2 = error correcting capability |
| 309 | + * d = 2*BCHCode_data->t + 1 = 5 = designed minimum distance |
| 310 | + * k = n - deg(g(x)) = 21 = dimension |
| 311 | + * g[] = coefficients of generator polynomial, g(x) [n - k + 1]=[11] |
| 312 | + * alpha_to [] = log table of GF(2**5) |
| 313 | + * index_of[] = antilog table of GF(2**5) |
| 314 | + * data[] = coefficients of data polynomial, i(x) |
| 315 | + * bb[] = coefficients of redundancy polynomial ( x**(10) i(x) ) modulo g(x) |
| 316 | + */ |
| 317 | +struct BCHCode * BCHCode_New(int p[], int m, int n, int k, int t) { |
| 318 | + struct BCHCode * BCHCode_data=NULL; |
| 319 | + |
| 320 | + BCHCode_data=(struct BCHCode *) malloc(sizeof (struct BCHCode)); |
| 321 | + |
| 322 | + if (BCHCode_data!=NULL) { |
| 323 | + BCHCode_data->alpha_to=(int *) malloc(sizeof(int) * (n+1)); |
| 324 | + BCHCode_data->index_of=(int *) malloc(sizeof(int) * (n+1)); |
| 325 | + BCHCode_data->p=(int *) malloc(sizeof(int) * (m+1)); |
| 326 | + BCHCode_data->g=(int *) malloc(sizeof(int) * (n-k+1)); |
| 327 | + BCHCode_data->bb=(int *) malloc(sizeof(int) * (n-k+1)); |
| 328 | + |
| 329 | + if ( |
| 330 | + BCHCode_data->alpha_to == NULL || |
| 331 | + BCHCode_data->index_of == NULL || |
| 332 | + BCHCode_data->p == NULL || |
| 333 | + BCHCode_data->g == NULL || |
| 334 | + BCHCode_data->bb == NULL |
| 335 | + ) { |
| 336 | + BCHCode_Delete(BCHCode_data); |
| 337 | + BCHCode_data=NULL; |
| 338 | + } |
| 339 | + } |
| 340 | + |
| 341 | + if (BCHCode_data!=NULL) { |
| 342 | + int i; |
| 343 | + for (i=0; i<(m+1); i++) { |
| 344 | + BCHCode_data->p[i]=p[i]; |
| 345 | + } |
| 346 | + BCHCode_data->m=m; |
| 347 | + BCHCode_data->n=n; |
| 348 | + BCHCode_data->k=k; |
| 349 | + BCHCode_data->t=t; |
| 350 | + |
| 351 | + generate_gf(BCHCode_data); /* generate the Galois Field GF(2**m) */ |
| 352 | + gen_poly(BCHCode_data); /* Compute the generator polynomial of BCH code */ |
| 353 | + } |
| 354 | + |
| 355 | + return BCHCode_data; |
| 356 | +} |
| 357 | + |
| 358 | +void BCHCode_Delete(struct BCHCode * BCHCode_data) { |
| 359 | + if (BCHCode_data==NULL) return; |
| 360 | + |
| 361 | + if (BCHCode_data->alpha_to != NULL) free(BCHCode_data->alpha_to); |
| 362 | + if (BCHCode_data->index_of != NULL) free(BCHCode_data->index_of); |
| 363 | + if (BCHCode_data->p != NULL) free(BCHCode_data->p); |
| 364 | + if (BCHCode_data->g != NULL) free(BCHCode_data->g); |
| 365 | + if (BCHCode_data->bb != NULL) free(BCHCode_data->bb); |
| 366 | + |
| 367 | + free(BCHCode_data); |
| 368 | +} |
0 commit comments