Skip to content
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.

Commit 26ac27b

Browse files
dtillerdtiller
committedMar 20, 2022
Added beginnings of BCH processing.
1 parent 3b17d45 commit 26ac27b

File tree

2 files changed

+1297
-0
lines changed

2 files changed

+1297
-0
lines changed
 

‎src/bch3.c

Lines changed: 677 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,677 @@
1+
/*
2+
* File: bch3.c
3+
* Title: Encoder/decoder for binary BCH codes in C (Version 3.1)
4+
* Author: Robert Morelos-Zaragoza
5+
* Date: August 1994
6+
* Revised: June 13, 1997
7+
*
8+
* =============== Encoder/Decoder for binary BCH codes in C =================
9+
*
10+
* Version 1: Original program. The user provides the generator polynomial
11+
* of the code (cumbersome!).
12+
* Version 2: Computes the generator polynomial of the code.
13+
* Version 3: No need to input the coefficients of a primitive polynomial of
14+
* degree m, used to construct the Galois Field GF(2**m). The
15+
* program now works for any binary BCH code of length such that:
16+
* 2**(m-1) - 1 < length <= 2**m - 1
17+
*
18+
* Note: You may have to change the size of the arrays to make it work.
19+
*
20+
* The encoding and decoding methods used in this program are based on the
21+
* book "Error Control Coding: Fundamentals and Applications", by Lin and
22+
* Costello, Prentice Hall, 1983.
23+
*
24+
* Thanks to Patrick Boyle (pboyle@era.com) for his observation that 'bch2.c'
25+
* did not work for lengths other than 2**m-1 which led to this new version.
26+
* Portions of this program are from 'rs.c', a Reed-Solomon encoder/decoder
27+
* in C, written by Simon Rockliff (simon@augean.ua.oz.au) on 21/9/89. The
28+
* previous version of the BCH encoder/decoder in C, 'bch2.c', was written by
29+
* Robert Morelos-Zaragoza (robert@spectra.eng.hawaii.edu) on 5/19/92.
30+
*
31+
* NOTE:
32+
* The author is not responsible for any malfunctioning of
33+
* this program, nor for any damage caused by it. Please include the
34+
* original program along with these comments in any redistribution.
35+
*
36+
* For more information, suggestions, or other ideas on implementing error
37+
* correcting codes, please contact me at:
38+
*
39+
* Robert Morelos-Zaragoza
40+
* 5120 Woodway, Suite 7036
41+
* Houston, Texas 77056
42+
*
43+
* email: r.morelos-zaragoza@ieee.org
44+
*
45+
* COPYRIGHT NOTICE: This computer program is free for non-commercial purposes.
46+
* You may implement this program for any non-commercial application. You may
47+
* also implement this program for commercial purposes, provided that you
48+
* obtain my written permission. Any modification of this program is covered
49+
* by this copyright.
50+
*
51+
* == Copyright (c) 1994-7, Robert Morelos-Zaragoza. All rights reserved. ==
52+
*
53+
* m = order of the Galois field GF(2**m)
54+
* n = 2**m - 1 = size of the multiplicative group of GF(2**m)
55+
* length = length of the BCH code
56+
* t = error correcting capability (max. no. of errors the code corrects)
57+
* d = 2*t + 1 = designed min. distance = no. of consecutive roots of g(x) + 1
58+
* k = n - deg(g(x)) = dimension (no. of information bits/codeword) of the code
59+
* p[] = coefficients of a primitive polynomial used to generate GF(2**m)
60+
* g[] = coefficients of the generator polynomial, g(x)
61+
* alpha_to [] = log table of GF(2**m)
62+
* index_of[] = antilog table of GF(2**m)
63+
* data[] = information bits = coefficients of data polynomial, i(x)
64+
* bb[] = coefficients of redundancy polynomial x^(length-k) i(x) modulo g(x)
65+
* numerr = number of errors
66+
* errpos[] = error positions
67+
* recd[] = coefficients of the received polynomial
68+
* decerror = number of decoding errors (in _message_ positions)
69+
*
70+
*/
71+
72+
#include <math.h>
73+
#include <stdio.h>
74+
#include <stdlib.h>
75+
#include <string.h>
76+
77+
int m, n, length, k, t, d;
78+
int p[21];
79+
int alpha_to[1048576], index_of[1048576], g[548576];
80+
int recd[1048576], data[1048576], bb[548576];
81+
int seed;
82+
int numerr, errpos[1024], decerror = 0;
83+
84+
int real_data[][45] = {
85+
{ 1,1,1,1,0,0,0,0,
86+
0,1,0,1,1,0,1,0,
87+
0,1,1,0,1,0,1,0,
88+
0,1,1,0,1,0,1,0,
89+
0,0,0,0,0,0,0,1,
90+
0,1,1,0,0},
91+
92+
// 0=f0 1=81 2=52 3=6b 4=71 5=a5 6=63 7=08
93+
{ 1,1,1,1,0,0,0,0,
94+
1,0,0,0,0,0,0,1,
95+
0,1,0,1,0,0,1,0,
96+
0,1,1,0,1,0,1,1,
97+
0,1,1,1,0,0,0,1,
98+
1,0,1,0,0},
99+
100+
// 0=f0 1=85 2=50 3=6a 4=01 5=e5 6=6e 7=84
101+
{ 1,1,1,1,0,0,0,0,
102+
1,0,0,0,0,1,0,1,
103+
0,1,0,1,0,0,0,0,
104+
0,1,1,0,1,0,1,0,
105+
0,0,0,0,0,0,0,1,
106+
1,1,1,0,0},
107+
108+
// 0=f0 1=85 2=59 3=5a 4=01 5=e5 6=6e 7=84
109+
{ 1,1,1,1,0,0,0,0,
110+
1,0,0,0,0,1,0,1,
111+
0,1,0,1,1,0,1,0,
112+
0,0,0,0,0,0,0,1,
113+
1,1,1,0,0},
114+
115+
// 0=f1 1=34 2=50 3=1a 4=01 5=e5 6=66 7=fe
116+
{ 1,1,1,1,0,0,0,1,
117+
0,0,1,1,0,1,0,0,
118+
0,1,0,1,0,0,0,0,
119+
0,0,0,1,1,0,1,0,
120+
0,0,0,0,0,0,0,1,
121+
1,1,1,0,0},
122+
// 0=f0 1=eb 2=10 3=ea 4=01 5=6e 6=54 7=1c
123+
{ 1,1,1,1,0,0,0,0,
124+
1,1,1,0,1,0,1,1,
125+
0,0,0,1,0,0,0,0,
126+
1,1,1,0,1,0,1,0,
127+
0,0,0,0,0,0,0,1,
128+
0,1,1,0,1},
129+
130+
// 0=f0 1=ea 2=5c 3=ea 4=01 5=6e 6=55 7=0e
131+
{ 1,1,1,1,0,0,0,0,
132+
1,1,1,0,1,0,1,0,
133+
0,1,0,1,1,1,0,0,
134+
1,1,1,0,1,0,1,0,
135+
0,0,0,0,0,0,0,1,
136+
0,1,1,0,1},
137+
};
138+
139+
int expected[][18] = {
140+
{ 0,1,1, 0,0,1,1, 0,0,1,1, 1,1,0,1, 0,0,0 },
141+
{ 1,0,1, 0,1,1,0, 0,0,1,1, 0,0,0,0, 1,0,0 },
142+
//orig { 1,0,1, 0,1,1,0, 1,1,1,0, 1,0,0,0, 0,1,0 },
143+
{ 1,0,1, 0,0,0,0, 0,1,1,0, 1,0,0,0, 0,1,0 }, // CORRECTED
144+
{ 1,0,1, 0,1,1,0, 1,1,1,0, 1,0,0,0, 0,1,0 },
145+
{ 1,0,1, 0,1,1,0, 0,1,1,0, 1,1,1,1, 1,1,1 },
146+
{ 1,1,0, 0,1,0,1, 0,1,0,0, 0,0,0,1, 1,1,0 },
147+
{ 1,1,0, 0,1,0,1, 0,1,0,1, 0,0,0,0, 1,1,1 },
148+
};
149+
150+
void
151+
read_p()
152+
/*
153+
* Read m, the degree of a primitive polynomial p(x) used to compute the
154+
* Galois field GF(2**m). Get precomputed coefficients p[] of p(x). Read
155+
* the code length.
156+
*/
157+
{
158+
int i, ninf;
159+
160+
printf("bch3: An encoder/decoder for binary BCH codes\n");
161+
printf("Copyright (c) 1994-7. Robert Morelos-Zaragoza.\n");
162+
printf("This program is free, please read first the copyright notice.\n");
163+
printf("\nFirst, enter a value of m such that the code length is\n");
164+
printf("2**(m-1) - 1 < length <= 2**m - 1\n\n");
165+
do {
166+
printf("Enter m (between 2 and 20): ");
167+
scanf("%d", &m);
168+
} while ( !(m>1) || !(m<21) );
169+
for (i=1; i<m; i++)
170+
p[i] = 0;
171+
p[0] = p[m] = 1;
172+
if (m == 2) p[1] = 1;
173+
else if (m == 3) p[1] = 1;
174+
else if (m == 4) p[1] = 1;
175+
else if (m == 5) p[2] = 1;
176+
else if (m == 6) p[1] = 1;
177+
else if (m == 7) p[1] = 1;
178+
else if (m == 8) p[4] = p[5] = p[6] = 1;
179+
else if (m == 9) p[4] = 1;
180+
else if (m == 10) p[3] = 1;
181+
else if (m == 11) p[2] = 1;
182+
else if (m == 12) p[3] = p[4] = p[7] = 1;
183+
else if (m == 13) p[1] = p[3] = p[4] = 1;
184+
else if (m == 14) p[1] = p[11] = p[12] = 1;
185+
else if (m == 15) p[1] = 1;
186+
else if (m == 16) p[2] = p[3] = p[5] = 1;
187+
else if (m == 17) p[3] = 1;
188+
else if (m == 18) p[7] = 1;
189+
else if (m == 19) p[1] = p[5] = p[6] = 1;
190+
else if (m == 20) p[3] = 1;
191+
printf("p(x) = ");
192+
n = 1;
193+
for (i = 0; i <= m; i++) {
194+
n *= 2;
195+
printf("%1d", p[i]);
196+
}
197+
printf("\n");
198+
n = n / 2 - 1;
199+
ninf = (n + 1) / 2 - 1;
200+
do {
201+
printf("Enter code length (%d < length <= %d): ", ninf, n);
202+
scanf("%d", &length);
203+
} while ( !((length <= n)&&(length>ninf)) );
204+
}
205+
206+
207+
void
208+
generate_gf()
209+
/*
210+
* Generate field GF(2**m) from the irreducible polynomial p(X) with
211+
* coefficients in p[0]..p[m].
212+
*
213+
* Lookup tables:
214+
* index->polynomial form: alpha_to[] contains j=alpha^i;
215+
* polynomial form -> index form: index_of[j=alpha^i] = i
216+
*
217+
* alpha=2 is the primitive element of GF(2**m)
218+
*/
219+
{
220+
register int i, mask;
221+
222+
mask = 1;
223+
alpha_to[m] = 0;
224+
for (i = 0; i < m; i++) {
225+
alpha_to[i] = mask;
226+
index_of[alpha_to[i]] = i;
227+
if (p[i] != 0)
228+
alpha_to[m] ^= mask;
229+
mask <<= 1;
230+
}
231+
index_of[alpha_to[m]] = m;
232+
mask >>= 1;
233+
for (i = m + 1; i < n; i++) {
234+
if (alpha_to[i - 1] >= mask)
235+
alpha_to[i] = alpha_to[m] ^ ((alpha_to[i - 1] ^ mask) << 1);
236+
else
237+
alpha_to[i] = alpha_to[i - 1] << 1;
238+
index_of[alpha_to[i]] = i;
239+
}
240+
index_of[0] = -1;
241+
}
242+
243+
244+
void
245+
gen_poly()
246+
/*
247+
* Compute the generator polynomial of a binary BCH code. Fist generate the
248+
* cycle sets modulo 2**m - 1, cycle[][] = (i, 2*i, 4*i, ..., 2^l*i). Then
249+
* determine those cycle sets that contain integers in the set of (d-1)
250+
* consecutive integers {1..(d-1)}. The generator polynomial is calculated
251+
* as the product of linear factors of the form (x+alpha^i), for every i in
252+
* the above cycle sets.
253+
*/
254+
{
255+
register int ii, jj, ll, kaux;
256+
register int test, aux, nocycles, root, noterms, rdncy;
257+
int cycle[1024][21], size[1024], min[1024], zeros[1024];
258+
259+
/* Generate cycle sets modulo n, n = 2**m - 1 */
260+
cycle[0][0] = 0;
261+
size[0] = 1;
262+
cycle[1][0] = 1;
263+
size[1] = 1;
264+
jj = 1; /* cycle set index */
265+
if (m > 9) {
266+
printf("Computing cycle sets modulo %d\n", n);
267+
printf("(This may take some time)...\n");
268+
}
269+
do {
270+
/* Generate the jj-th cycle set */
271+
ii = 0;
272+
do {
273+
ii++;
274+
cycle[jj][ii] = (cycle[jj][ii - 1] * 2) % n;
275+
size[jj]++;
276+
aux = (cycle[jj][ii] * 2) % n;
277+
} while (aux != cycle[jj][0]);
278+
/* Next cycle set representative */
279+
ll = 0;
280+
do {
281+
ll++;
282+
test = 0;
283+
for (ii = 1; ((ii <= jj) && (!test)); ii++)
284+
/* Examine previous cycle sets */
285+
for (kaux = 0; ((kaux < size[ii]) && (!test)); kaux++)
286+
if (ll == cycle[ii][kaux])
287+
test = 1;
288+
} while ((test) && (ll < (n - 1)));
289+
if (!(test)) {
290+
jj++; /* next cycle set index */
291+
cycle[jj][0] = ll;
292+
size[jj] = 1;
293+
}
294+
} while (ll < (n - 1));
295+
nocycles = jj; /* number of cycle sets modulo n */
296+
297+
printf("Enter the error correcting capability, t: ");
298+
scanf("%d", &t);
299+
300+
d = 2 * t + 1;
301+
302+
/* Search for roots 1, 2, ..., d-1 in cycle sets */
303+
kaux = 0;
304+
rdncy = 0;
305+
for (ii = 1; ii <= nocycles; ii++) {
306+
min[kaux] = 0;
307+
test = 0;
308+
for (jj = 0; ((jj < size[ii]) && (!test)); jj++)
309+
for (root = 1; ((root < d) && (!test)); root++)
310+
if (root == cycle[ii][jj]) {
311+
test = 1;
312+
min[kaux] = ii;
313+
}
314+
if (min[kaux]) {
315+
rdncy += size[min[kaux]];
316+
kaux++;
317+
}
318+
}
319+
noterms = kaux;
320+
kaux = 1;
321+
for (ii = 0; ii < noterms; ii++)
322+
for (jj = 0; jj < size[min[ii]]; jj++) {
323+
zeros[kaux] = cycle[min[ii]][jj];
324+
kaux++;
325+
}
326+
327+
k = length - rdncy;
328+
329+
if (k<0)
330+
{
331+
printf("Parameters invalid!\n");
332+
exit(0);
333+
}
334+
335+
printf("This is a (%d, %d, %d) binary BCH code\n", length, k, d);
336+
337+
/* Compute the generator polynomial */
338+
g[0] = alpha_to[zeros[1]];
339+
g[1] = 1; /* g(x) = (X + zeros[1]) initially */
340+
for (ii = 2; ii <= rdncy; ii++) {
341+
g[ii] = 1;
342+
for (jj = ii - 1; jj > 0; jj--)
343+
if (g[jj] != 0)
344+
g[jj] = g[jj - 1] ^ alpha_to[(index_of[g[jj]] + zeros[ii]) % n];
345+
else
346+
g[jj] = g[jj - 1];
347+
g[0] = alpha_to[(index_of[g[0]] + zeros[ii]) % n];
348+
}
349+
printf("Generator polynomial:\ng(x) = ");
350+
for (ii = 0; ii <= rdncy; ii++) {
351+
printf("%d", g[ii]);
352+
if (ii && ((ii % 50) == 0))
353+
printf("\n");
354+
}
355+
printf("\n");
356+
}
357+
358+
359+
void
360+
encode_bch()
361+
/*
362+
* Compute redundacy bb[], the coefficients of b(x). The redundancy
363+
* polynomial b(x) is the remainder after dividing x^(length-k)*data(x)
364+
* by the generator polynomial g(x).
365+
*/
366+
{
367+
register int i, j;
368+
register int feedback;
369+
370+
for (i = 0; i < length - k; i++)
371+
bb[i] = 0;
372+
for (i = k - 1; i >= 0; i--) {
373+
feedback = data[i] ^ bb[length - k - 1];
374+
if (feedback != 0) {
375+
for (j = length - k - 1; j > 0; j--)
376+
if (g[j] != 0)
377+
bb[j] = bb[j - 1] ^ feedback;
378+
else
379+
bb[j] = bb[j - 1];
380+
bb[0] = g[0] && feedback;
381+
} else {
382+
for (j = length - k - 1; j > 0; j--)
383+
bb[j] = bb[j - 1];
384+
bb[0] = 0;
385+
}
386+
}
387+
}
388+
389+
390+
void
391+
decode_bch()
392+
/*
393+
* Simon Rockliff's implementation of Berlekamp's algorithm.
394+
*
395+
* Assume we have received bits in recd[i], i=0..(n-1).
396+
*
397+
* Compute the 2*t syndromes by substituting alpha^i into rec(X) and
398+
* evaluating, storing the syndromes in s[i], i=1..2t (leave s[0] zero) .
399+
* Then we use the Berlekamp algorithm to find the error location polynomial
400+
* elp[i].
401+
*
402+
* If the degree of the elp is >t, then we cannot correct all the errors, and
403+
* we have detected an uncorrectable error pattern. We output the information
404+
* bits uncorrected.
405+
*
406+
* If the degree of elp is <=t, we substitute alpha^i , i=1..n into the elp
407+
* to get the roots, hence the inverse roots, the error location numbers.
408+
* This step is usually called "Chien's search".
409+
*
410+
* If the number of errors located is not equal the degree of the elp, then
411+
* the decoder assumes that there are more than t errors and cannot correct
412+
* them, only detect them. We output the information bits uncorrected.
413+
*/
414+
{
415+
register int i, j, u, q, t2, count = 0, syn_error = 0;
416+
int elp[1026][1024], d[1026], l[1026], u_lu[1026], s[1025];
417+
int root[200], loc[200], err[1024], reg[201];
418+
419+
t2 = 2 * t;
420+
421+
/* first form the syndromes */
422+
printf("S(x) = ");
423+
for (i = 1; i <= t2; i++) {
424+
s[i] = 0;
425+
for (j = 0; j < length; j++)
426+
if (recd[j] != 0)
427+
s[i] ^= alpha_to[(i * j) % n];
428+
if (s[i] != 0)
429+
syn_error = 1; /* set error flag if non-zero syndrome */
430+
/*
431+
* Note: If the code is used only for ERROR DETECTION, then
432+
* exit program here indicating the presence of errors.
433+
*/
434+
/* convert syndrome from polynomial form to index form */
435+
s[i] = index_of[s[i]];
436+
printf("%3d ", s[i]);
437+
}
438+
printf("\n");
439+
440+
if (syn_error) { /* if there are errors, try to correct them */
441+
/*
442+
* Compute the error location polynomial via the Berlekamp
443+
* iterative algorithm. Following the terminology of Lin and
444+
* Costello's book : d[u] is the 'mu'th discrepancy, where
445+
* u='mu'+1 and 'mu' (the Greek letter!) is the step number
446+
* ranging from -1 to 2*t (see L&C), l[u] is the degree of
447+
* the elp at that step, and u_l[u] is the difference between
448+
* the step number and the degree of the elp.
449+
*/
450+
/* initialise table entries */
451+
d[0] = 0; /* index form */
452+
d[1] = s[1]; /* index form */
453+
elp[0][0] = 0; /* index form */
454+
elp[1][0] = 1; /* polynomial form */
455+
for (i = 1; i < t2; i++) {
456+
elp[0][i] = -1; /* index form */
457+
elp[1][i] = 0; /* polynomial form */
458+
}
459+
l[0] = 0;
460+
l[1] = 0;
461+
u_lu[0] = -1;
462+
u_lu[1] = 0;
463+
u = 0;
464+
465+
do {
466+
u++;
467+
if (d[u] == -1) {
468+
l[u + 1] = l[u];
469+
for (i = 0; i <= l[u]; i++) {
470+
elp[u + 1][i] = elp[u][i];
471+
elp[u][i] = index_of[elp[u][i]];
472+
}
473+
} else
474+
/*
475+
* search for words with greatest u_lu[q] for
476+
* which d[q]!=0
477+
*/
478+
{
479+
q = u - 1;
480+
while ((d[q] == -1) && (q > 0))
481+
q--;
482+
/* have found first non-zero d[q] */
483+
if (q > 0) {
484+
j = q;
485+
do {
486+
j--;
487+
if ((d[j] != -1) && (u_lu[q] < u_lu[j]))
488+
q = j;
489+
} while (j > 0);
490+
}
491+
492+
/*
493+
* have now found q such that d[u]!=0 and
494+
* u_lu[q] is maximum
495+
*/
496+
/* store degree of new elp polynomial */
497+
if (l[u] > l[q] + u - q)
498+
l[u + 1] = l[u];
499+
else
500+
l[u + 1] = l[q] + u - q;
501+
502+
/* form new elp(x) */
503+
for (i = 0; i < t2; i++)
504+
elp[u + 1][i] = 0;
505+
for (i = 0; i <= l[q]; i++)
506+
if (elp[q][i] != -1)
507+
elp[u + 1][i + u - q] =
508+
alpha_to[(d[u] + n - d[q] + elp[q][i]) % n];
509+
for (i = 0; i <= l[u]; i++) {
510+
elp[u + 1][i] ^= elp[u][i];
511+
elp[u][i] = index_of[elp[u][i]];
512+
}
513+
}
514+
u_lu[u + 1] = u - l[u + 1];
515+
516+
/* form (u+1)th discrepancy */
517+
if (u < t2) {
518+
/* no discrepancy computed on last iteration */
519+
if (s[u + 1] != -1)
520+
d[u + 1] = alpha_to[s[u + 1]];
521+
else
522+
d[u + 1] = 0;
523+
for (i = 1; i <= l[u + 1]; i++)
524+
if ((s[u + 1 - i] != -1) && (elp[u + 1][i] != 0))
525+
d[u + 1] ^= alpha_to[(s[u + 1 - i]
526+
+ index_of[elp[u + 1][i]]) % n];
527+
/* put d[u+1] into index form */
528+
d[u + 1] = index_of[d[u + 1]];
529+
}
530+
} while ((u < t2) && (l[u + 1] <= t));
531+
532+
u++;
533+
if (l[u] <= t) {/* Can correct errors */
534+
/* put elp into index form */
535+
for (i = 0; i <= l[u]; i++)
536+
elp[u][i] = index_of[elp[u][i]];
537+
538+
printf("sigma(x) = ");
539+
for (i = 0; i <= l[u]; i++)
540+
printf("%3d ", elp[u][i]);
541+
printf("\n");
542+
printf("Roots: ");
543+
544+
/* Chien search: find roots of the error location polynomial */
545+
for (i = 1; i <= l[u]; i++)
546+
reg[i] = elp[u][i];
547+
count = 0;
548+
for (i = 1; i <= n; i++) {
549+
q = 1;
550+
for (j = 1; j <= l[u]; j++)
551+
if (reg[j] != -1) {
552+
reg[j] = (reg[j] + j) % n;
553+
q ^= alpha_to[reg[j]];
554+
}
555+
if (!q) { /* store root and error
556+
* location number indices */
557+
root[count] = i;
558+
loc[count] = n - i;
559+
count++;
560+
printf("%3d ", n - i);
561+
}
562+
}
563+
printf("\n");
564+
if (count == l[u])
565+
/* no. roots = degree of elp hence <= t errors */
566+
for (i = 0; i < l[u]; i++)
567+
recd[loc[i]] ^= 1;
568+
else /* elp has degree >t hence cannot solve */
569+
printf("Incomplete decoding: errors detected\n");
570+
}
571+
}
572+
}
573+
574+
575+
576+
int main()
577+
{
578+
int i;
579+
580+
read_p(); /* Read m */
581+
generate_gf(); /* Construct the Galois Field GF(2**m) */
582+
gen_poly(); /* Compute the generator polynomial of BCH code */
583+
584+
#ifdef TEST
585+
for (int count = 0; count < sizeof(real_data) / sizeof(*real_data); count++) {
586+
memcpy(data, real_data[count], sizeof(*real_data));
587+
#else
588+
/* Randomly generate DATA */
589+
seed = 131073;
590+
srandom(seed);
591+
for (i = 0; i < k; i++)
592+
data[i] = ( random() & 65536 ) >> 16;
593+
#endif
594+
encode_bch(); /* encode data */
595+
596+
/*
597+
* recd[] are the coefficients of c(x) = x**(length-k)*data(x) + b(x)
598+
*/
599+
for (i = 0; i < length - k; i++)
600+
recd[i] = bb[i];
601+
for (i = 0; i < k; i++)
602+
recd[i + length - k] = data[i];
603+
printf("Code polynomial:\nc(x) = ");
604+
for (i = 0; i < length; i++) {
605+
printf("%1d", recd[i]);
606+
if (i && ((i % 50) == 0))
607+
printf("\n");
608+
}
609+
printf("\n");
610+
611+
printf("Enter the number of errors:\n");
612+
scanf("%d", &numerr); /* CHANNEL errors */
613+
printf("Enter error locations (integers between");
614+
printf(" 0 and %d): ", length-1);
615+
/*
616+
* recd[] are the coefficients of r(x) = c(x) + e(x)
617+
*/
618+
for (i = 0; i < numerr; i++)
619+
scanf("%d", &errpos[i]);
620+
if (numerr)
621+
for (i = 0; i < numerr; i++)
622+
recd[errpos[i]] ^= 1;
623+
printf("r(x) = ");
624+
for (i = 0; i < length; i++) {
625+
printf("%1d", recd[i]);
626+
if (i == length - k - 1) printf(" ");
627+
//if (i && ((i % 50) == 0))
628+
//printf("\n");
629+
}
630+
printf("\n");
631+
632+
decode_bch(); /* DECODE received codeword recv[] */
633+
634+
/*
635+
* print out original and decoded data
636+
*/
637+
printf("Results:\n");
638+
printf("original data = ");
639+
for (i = 0; i < k; i++) {
640+
printf("%1d", data[i]);
641+
if (i && ((i % 50) == 0))
642+
printf("\n");
643+
}
644+
printf("\nrecovered data = ");
645+
for (i = length - k; i < length; i++) {
646+
printf("%1d", recd[i]);
647+
if ((i-length+k) && (((i-length+k) % 50) == 0))
648+
printf("\n");
649+
}
650+
printf("\n");
651+
652+
int flag = 0;
653+
printf("n=%d, k=%d, length=%d\n", n, k, length);
654+
#ifdef TEST
655+
for (int jj = 0; jj < n - k; jj++) {
656+
if (expected[count][jj] != recd[jj]) {
657+
printf("bit %d: expected %d calc: %d\n", jj, expected[count][jj], recd[jj]);
658+
flag++;
659+
}
660+
}
661+
printf("%d ERRORS.\n", flag);
662+
#endif
663+
664+
/*
665+
* DECODING ERRORS? we compare only the data portion
666+
*/
667+
for (i = length - k; i < length; i++)
668+
if (data[i - length + k] != recd[i])
669+
decerror++;
670+
if (decerror)
671+
printf("There were %d decoding errors in message positions\n", decerror);
672+
else
673+
printf("Succesful decoding\n");
674+
#ifdef TEST
675+
}
676+
#endif
677+
}

‎src/bch3a.c

Lines changed: 620 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,620 @@
1+
/*
2+
* File: bch3.c
3+
* Title: Encoder/decoder for binary BCH codes in C (Version 3.1)
4+
* Author: Robert Morelos-Zaragoza
5+
* Date: August 1994
6+
* Revised: June 13, 1997
7+
*
8+
* =============== Encoder/Decoder for binary BCH codes in C =================
9+
*
10+
* Version 1: Original program. The user provides the generator polynomial
11+
* of the code (cumbersome!).
12+
* Version 2: Computes the generator polynomial of the code.
13+
* Version 3: No need to input the coefficients of a primitive polynomial of
14+
* degree m, used to construct the Galois Field GF(2**m). The
15+
* program now works for any binary BCH code of length such that:
16+
* 2**(m-1) - 1 < length <= 2**m - 1
17+
*
18+
* Note: You may have to change the size of the arrays to make it work.
19+
*
20+
* The encoding and decoding methods used in this program are based on the
21+
* book "Error Control Coding: Fundamentals and Applications", by Lin and
22+
* Costello, Prentice Hall, 1983.
23+
*
24+
* Thanks to Patrick Boyle (pboyle@era.com) for his observation that 'bch2.c'
25+
* did not work for lengths other than 2**m-1 which led to this new version.
26+
* Portions of this program are from 'rs.c', a Reed-Solomon encoder/decoder
27+
* in C, written by Simon Rockliff (simon@augean.ua.oz.au) on 21/9/89. The
28+
* previous version of the BCH encoder/decoder in C, 'bch2.c', was written by
29+
* Robert Morelos-Zaragoza (robert@spectra.eng.hawaii.edu) on 5/19/92.
30+
*
31+
* NOTE:
32+
* The author is not responsible for any malfunctioning of
33+
* this program, nor for any damage caused by it. Please include the
34+
* original program along with these comments in any redistribution.
35+
*
36+
* For more information, suggestions, or other ideas on implementing error
37+
* correcting codes, please contact me at:
38+
*
39+
* Robert Morelos-Zaragoza
40+
* 5120 Woodway, Suite 7036
41+
* Houston, Texas 77056
42+
*
43+
* email: r.morelos-zaragoza@ieee.org
44+
*
45+
* COPYRIGHT NOTICE: This computer program is free for non-commercial purposes.
46+
* You may implement this program for any non-commercial application. You may
47+
* also implement this program for commercial purposes, provided that you
48+
* obtain my written permission. Any modification of this program is covered
49+
* by this copyright.
50+
*
51+
* == Copyright (c) 1994-7, Robert Morelos-Zaragoza. All rights reserved. ==
52+
*
53+
* m = order of the Galois field GF(2**m)
54+
* n = 2**m - 1 = size of the multiplicative group of GF(2**m)
55+
* length = length of the BCH code
56+
* t = error correcting capability (max. no. of errors the code corrects)
57+
* d = 2*t + 1 = designed min. distance = no. of consecutive roots of g(x) + 1
58+
* k = n - deg(g(x)) = dimension (no. of information bits/codeword) of the code
59+
* p[] = coefficients of a primitive polynomial used to generate GF(2**m)
60+
* g[] = coefficients of the generator polynomial, g(x)
61+
* alpha_to [] = log table of GF(2**m)
62+
* index_of[] = antilog table of GF(2**m)
63+
* data[] = information bits = coefficients of data polynomial, i(x)
64+
* bb[] = coefficients of redundancy polynomial x^(length-k) i(x) modulo g(x)
65+
* numerr = number of errors
66+
* errpos[] = error positions
67+
* recd[] = coefficients of the received polynomial
68+
* decerror = number of decoding errors (in _message_ positions)
69+
*
70+
*/
71+
72+
#include <math.h>
73+
#include <stdio.h>
74+
#include <stdlib.h>
75+
#include <string.h>
76+
77+
int m, n, length, k, t, d;
78+
int p[21];
79+
int alpha_to[1048576], index_of[1048576], g[548576];
80+
int recd[1048576], data[1048576], bb[548576];
81+
int seed;
82+
int numerr, errpos[1024], decerror = 0;
83+
int orig_recd[1048576];
84+
85+
uint64_t packets[] = {
86+
87+
#define ALL_DATA
88+
/* GOOD TEST */ 0xb217a2b953ddc552, /* random example from bch3 program */
89+
#ifdef ALL_DATA
90+
0xf05a6a6a016333d0, /* g001-cut-lenthened_457.938M.wav */
91+
0xf081526b71a56308, /* 1st in eotd_received_data */
92+
/* 3 errors */ 0xf085506a01e56e84, /* 2nd in eotd_received_data */
93+
/* fixed */ 0xf085506a01e50684, /* 2nd, but with the bits fixed */
94+
#endif
95+
#ifdef ALL_DATA
96+
0xf085595a01e56e84, /* 3rd */
97+
0xf134501a01e566fe, /* 4th */
98+
0xf0eb10ea016e541c, /* 5th */
99+
0xf0ea5cea016e550e, /* 6th */
100+
0xe021101a0132bce4, /* Sun Mar 20 05:41:00 2022 */
101+
0xf042505bcfd564e4, /* Sun Mar 20 12:58:43 2022 */
102+
0xf08c10aa01737b1a, /* Sun Mar 20 13:35:48 2022 */
103+
0xf08c10b1c0e09064, /* Sun Mar 20 13:37:05 2022 */
104+
0xf08c106a01647ae8, /* Sun Mar 20 13:37:48 2022 */
105+
0x508c126a01647ae8,
106+
#endif
107+
};
108+
109+
void
110+
read_p()
111+
/*
112+
* Read m, the degree of a primitive polynomial p(x) used to compute the
113+
* Galois field GF(2**m). Get precomputed coefficients p[] of p(x). Read
114+
* the code length.
115+
*/
116+
{
117+
int i, ninf;
118+
119+
printf("bch3: An encoder/decoder for binary BCH codes\n");
120+
printf("Copyright (c) 1994-7. Robert Morelos-Zaragoza.\n");
121+
printf("This program is free, please read first the copyright notice.\n");
122+
printf("\nFirst, enter a value of m such that the code length is\n");
123+
printf("2**(m-1) - 1 < length <= 2**m - 1\n\n");
124+
do {
125+
//printf("Enter m (between 2 and 20): ");
126+
//scanf("%d", &m);
127+
m = 6;
128+
} while ( !(m>1) || !(m<21) );
129+
for (i=1; i<m; i++)
130+
p[i] = 0;
131+
p[0] = p[m] = 1;
132+
if (m == 2) p[1] = 1;
133+
else if (m == 3) p[1] = 1;
134+
else if (m == 4) p[1] = 1;
135+
else if (m == 5) p[2] = 1;
136+
else if (m == 6) p[1] = 1;
137+
else if (m == 7) p[1] = 1;
138+
else if (m == 8) p[4] = p[5] = p[6] = 1;
139+
else if (m == 9) p[4] = 1;
140+
else if (m == 10) p[3] = 1;
141+
else if (m == 11) p[2] = 1;
142+
else if (m == 12) p[3] = p[4] = p[7] = 1;
143+
else if (m == 13) p[1] = p[3] = p[4] = 1;
144+
else if (m == 14) p[1] = p[11] = p[12] = 1;
145+
else if (m == 15) p[1] = 1;
146+
else if (m == 16) p[2] = p[3] = p[5] = 1;
147+
else if (m == 17) p[3] = 1;
148+
else if (m == 18) p[7] = 1;
149+
else if (m == 19) p[1] = p[5] = p[6] = 1;
150+
else if (m == 20) p[3] = 1;
151+
printf("p(x) = ");
152+
n = 1;
153+
for (i = 0; i <= m; i++) {
154+
n *= 2;
155+
printf("%1d", p[i]);
156+
}
157+
printf("\n");
158+
n = n / 2 - 1;
159+
ninf = (n + 1) / 2 - 1;
160+
/* do {
161+
printf("Enter code length (%d < length <= %d): ", ninf, n);
162+
scanf("%d", &length);
163+
} while ( !((length <= n)&&(length>ninf)) ); */
164+
length = 63;
165+
}
166+
167+
168+
void
169+
generate_gf()
170+
/*
171+
* Generate field GF(2**m) from the irreducible polynomial p(X) with
172+
* coefficients in p[0]..p[m].
173+
*
174+
* Lookup tables:
175+
* index->polynomial form: alpha_to[] contains j=alpha^i;
176+
* polynomial form -> index form: index_of[j=alpha^i] = i
177+
*
178+
* alpha=2 is the primitive element of GF(2**m)
179+
*/
180+
{
181+
register int i, mask;
182+
183+
mask = 1;
184+
alpha_to[m] = 0;
185+
for (i = 0; i < m; i++) {
186+
alpha_to[i] = mask;
187+
index_of[alpha_to[i]] = i;
188+
if (p[i] != 0)
189+
alpha_to[m] ^= mask;
190+
mask <<= 1;
191+
}
192+
index_of[alpha_to[m]] = m;
193+
mask >>= 1;
194+
for (i = m + 1; i < n; i++) {
195+
if (alpha_to[i - 1] >= mask)
196+
alpha_to[i] = alpha_to[m] ^ ((alpha_to[i - 1] ^ mask) << 1);
197+
else
198+
alpha_to[i] = alpha_to[i - 1] << 1;
199+
index_of[alpha_to[i]] = i;
200+
}
201+
index_of[0] = -1;
202+
}
203+
204+
205+
void
206+
gen_poly()
207+
/*
208+
* Compute the generator polynomial of a binary BCH code. Fist generate the
209+
* cycle sets modulo 2**m - 1, cycle[][] = (i, 2*i, 4*i, ..., 2^l*i). Then
210+
* determine those cycle sets that contain integers in the set of (d-1)
211+
* consecutive integers {1..(d-1)}. The generator polynomial is calculated
212+
* as the product of linear factors of the form (x+alpha^i), for every i in
213+
* the above cycle sets.
214+
*/
215+
{
216+
register int ii, jj, ll, kaux;
217+
register int test, aux, nocycles, root, noterms, rdncy;
218+
int cycle[1024][21], size[1024], min[1024], zeros[1024];
219+
220+
/* Generate cycle sets modulo n, n = 2**m - 1 */
221+
cycle[0][0] = 0;
222+
size[0] = 1;
223+
cycle[1][0] = 1;
224+
size[1] = 1;
225+
jj = 1; /* cycle set index */
226+
if (m > 9) {
227+
printf("Computing cycle sets modulo %d\n", n);
228+
printf("(This may take some time)...\n");
229+
}
230+
do {
231+
/* Generate the jj-th cycle set */
232+
ii = 0;
233+
do {
234+
ii++;
235+
cycle[jj][ii] = (cycle[jj][ii - 1] * 2) % n;
236+
size[jj]++;
237+
aux = (cycle[jj][ii] * 2) % n;
238+
} while (aux != cycle[jj][0]);
239+
/* Next cycle set representative */
240+
ll = 0;
241+
do {
242+
ll++;
243+
test = 0;
244+
for (ii = 1; ((ii <= jj) && (!test)); ii++)
245+
/* Examine previous cycle sets */
246+
for (kaux = 0; ((kaux < size[ii]) && (!test)); kaux++)
247+
if (ll == cycle[ii][kaux])
248+
test = 1;
249+
} while ((test) && (ll < (n - 1)));
250+
if (!(test)) {
251+
jj++; /* next cycle set index */
252+
cycle[jj][0] = ll;
253+
size[jj] = 1;
254+
}
255+
} while (ll < (n - 1));
256+
nocycles = jj; /* number of cycle sets modulo n */
257+
258+
//printf("Enter the error correcting capability, t: ");
259+
//scanf("%d", &t);
260+
t =3;
261+
262+
d = 2 * t + 1;
263+
264+
/* Search for roots 1, 2, ..., d-1 in cycle sets */
265+
kaux = 0;
266+
rdncy = 0;
267+
for (ii = 1; ii <= nocycles; ii++) {
268+
min[kaux] = 0;
269+
test = 0;
270+
for (jj = 0; ((jj < size[ii]) && (!test)); jj++)
271+
for (root = 1; ((root < d) && (!test)); root++)
272+
if (root == cycle[ii][jj]) {
273+
test = 1;
274+
min[kaux] = ii;
275+
}
276+
if (min[kaux]) {
277+
rdncy += size[min[kaux]];
278+
kaux++;
279+
}
280+
}
281+
noterms = kaux;
282+
kaux = 1;
283+
for (ii = 0; ii < noterms; ii++)
284+
for (jj = 0; jj < size[min[ii]]; jj++) {
285+
zeros[kaux] = cycle[min[ii]][jj];
286+
kaux++;
287+
}
288+
289+
k = length - rdncy;
290+
291+
if (k<0)
292+
{
293+
printf("Parameters invalid!\n");
294+
exit(0);
295+
}
296+
297+
printf("This is a (%d, %d, %d) binary BCH code\n", length, k, d);
298+
299+
/* Compute the generator polynomial */
300+
g[0] = alpha_to[zeros[1]];
301+
g[1] = 1; /* g(x) = (X + zeros[1]) initially */
302+
for (ii = 2; ii <= rdncy; ii++) {
303+
g[ii] = 1;
304+
for (jj = ii - 1; jj > 0; jj--)
305+
if (g[jj] != 0)
306+
g[jj] = g[jj - 1] ^ alpha_to[(index_of[g[jj]] + zeros[ii]) % n];
307+
else
308+
g[jj] = g[jj - 1];
309+
g[0] = alpha_to[(index_of[g[0]] + zeros[ii]) % n];
310+
}
311+
printf("Generator polynomial:\ng(x) = ");
312+
for (ii = 0; ii <= rdncy; ii++) {
313+
printf("%d", g[ii]);
314+
}
315+
printf("\n");
316+
}
317+
318+
319+
void
320+
encode_bch()
321+
/*
322+
* Compute redundacy bb[], the coefficients of b(x). The redundancy
323+
* polynomial b(x) is the remainder after dividing x^(length-k)*data(x)
324+
* by the generator polynomial g(x).
325+
*/
326+
{
327+
register int i, j;
328+
register int feedback;
329+
330+
for (i = 0; i < length - k; i++)
331+
bb[i] = 0;
332+
for (i = k - 1; i >= 0; i--) {
333+
feedback = data[i] ^ bb[length - k - 1];
334+
if (feedback != 0) {
335+
for (j = length - k - 1; j > 0; j--)
336+
if (g[j] != 0)
337+
bb[j] = bb[j - 1] ^ feedback;
338+
else
339+
bb[j] = bb[j - 1];
340+
bb[0] = g[0] && feedback;
341+
} else {
342+
for (j = length - k - 1; j > 0; j--)
343+
bb[j] = bb[j - 1];
344+
bb[0] = 0;
345+
}
346+
}
347+
}
348+
349+
350+
/* zero = success */
351+
int
352+
decode_bch()
353+
/*
354+
* Simon Rockliff's implementation of Berlekamp's algorithm.
355+
*
356+
* Assume we have received bits in recd[i], i=0..(n-1).
357+
*
358+
* Compute the 2*t syndromes by substituting alpha^i into rec(X) and
359+
* evaluating, storing the syndromes in s[i], i=1..2t (leave s[0] zero) .
360+
* Then we use the Berlekamp algorithm to find the error location polynomial
361+
* elp[i].
362+
*
363+
* If the degree of the elp is >t, then we cannot correct all the errors, and
364+
* we have detected an uncorrectable error pattern. We output the information
365+
* bits uncorrected.
366+
*
367+
* If the degree of elp is <=t, we substitute alpha^i , i=1..n into the elp
368+
* to get the roots, hence the inverse roots, the error location numbers.
369+
* This step is usually called "Chien's search".
370+
*
371+
* If the number of errors located is not equal the degree of the elp, then
372+
* the decoder assumes that there are more than t errors and cannot correct
373+
* them, only detect them. We output the information bits uncorrected.
374+
*/
375+
{
376+
register int i, j, u, q, t2, count = 0, syn_error = 0;
377+
int elp[1026][1024], d[1026], l[1026], u_lu[1026], s[1025];
378+
int root[200], loc[200], err[1024], reg[201];
379+
380+
t2 = 2 * t;
381+
382+
/* first form the syndromes */
383+
printf("S(x) = ");
384+
for (i = 1; i <= t2; i++) {
385+
s[i] = 0;
386+
for (j = 0; j < length; j++)
387+
if (recd[j] != 0)
388+
s[i] ^= alpha_to[(i * j) % n];
389+
if (s[i] != 0)
390+
syn_error = 1; /* set error flag if non-zero syndrome */
391+
/*
392+
* Note: If the code is used only for ERROR DETECTION, then
393+
* exit program here indicating the presence of errors.
394+
*/
395+
/* convert syndrome from polynomial form to index form */
396+
s[i] = index_of[s[i]];
397+
printf("%3d ", s[i]);
398+
}
399+
printf("\n");
400+
401+
if (syn_error) { /* if there are errors, try to correct them */
402+
/*
403+
* Compute the error location polynomial via the Berlekamp
404+
* iterative algorithm. Following the terminology of Lin and
405+
* Costello's book : d[u] is the 'mu'th discrepancy, where
406+
* u='mu'+1 and 'mu' (the Greek letter!) is the step number
407+
* ranging from -1 to 2*t (see L&C), l[u] is the degree of
408+
* the elp at that step, and u_l[u] is the difference between
409+
* the step number and the degree of the elp.
410+
*/
411+
/* initialise table entries */
412+
d[0] = 0; /* index form */
413+
d[1] = s[1]; /* index form */
414+
elp[0][0] = 0; /* index form */
415+
elp[1][0] = 1; /* polynomial form */
416+
for (i = 1; i < t2; i++) {
417+
elp[0][i] = -1; /* index form */
418+
elp[1][i] = 0; /* polynomial form */
419+
}
420+
l[0] = 0;
421+
l[1] = 0;
422+
u_lu[0] = -1;
423+
u_lu[1] = 0;
424+
u = 0;
425+
426+
do {
427+
u++;
428+
if (d[u] == -1) {
429+
l[u + 1] = l[u];
430+
for (i = 0; i <= l[u]; i++) {
431+
elp[u + 1][i] = elp[u][i];
432+
elp[u][i] = index_of[elp[u][i]];
433+
}
434+
} else
435+
/*
436+
* search for words with greatest u_lu[q] for
437+
* which d[q]!=0
438+
*/
439+
{
440+
q = u - 1;
441+
while ((d[q] == -1) && (q > 0))
442+
q--;
443+
/* have found first non-zero d[q] */
444+
if (q > 0) {
445+
j = q;
446+
do {
447+
j--;
448+
if ((d[j] != -1) && (u_lu[q] < u_lu[j]))
449+
q = j;
450+
} while (j > 0);
451+
}
452+
453+
/*
454+
* have now found q such that d[u]!=0 and
455+
* u_lu[q] is maximum
456+
*/
457+
/* store degree of new elp polynomial */
458+
if (l[u] > l[q] + u - q)
459+
l[u + 1] = l[u];
460+
else
461+
l[u + 1] = l[q] + u - q;
462+
463+
/* form new elp(x) */
464+
for (i = 0; i < t2; i++)
465+
elp[u + 1][i] = 0;
466+
for (i = 0; i <= l[q]; i++)
467+
if (elp[q][i] != -1)
468+
elp[u + 1][i + u - q] =
469+
alpha_to[(d[u] + n - d[q] + elp[q][i]) % n];
470+
for (i = 0; i <= l[u]; i++) {
471+
elp[u + 1][i] ^= elp[u][i];
472+
elp[u][i] = index_of[elp[u][i]];
473+
}
474+
}
475+
u_lu[u + 1] = u - l[u + 1];
476+
477+
/* form (u+1)th discrepancy */
478+
if (u < t2) {
479+
/* no discrepancy computed on last iteration */
480+
if (s[u + 1] != -1)
481+
d[u + 1] = alpha_to[s[u + 1]];
482+
else
483+
d[u + 1] = 0;
484+
for (i = 1; i <= l[u + 1]; i++)
485+
if ((s[u + 1 - i] != -1) && (elp[u + 1][i] != 0))
486+
d[u + 1] ^= alpha_to[(s[u + 1 - i]
487+
+ index_of[elp[u + 1][i]]) % n];
488+
/* put d[u+1] into index form */
489+
d[u + 1] = index_of[d[u + 1]];
490+
}
491+
} while ((u < t2) && (l[u + 1] <= t));
492+
493+
u++;
494+
if (l[u] <= t) {/* Can correct errors */
495+
/* put elp into index form */
496+
for (i = 0; i <= l[u]; i++)
497+
elp[u][i] = index_of[elp[u][i]];
498+
499+
printf("sigma(x) = ");
500+
for (i = 0; i <= l[u]; i++)
501+
printf("%3d ", elp[u][i]);
502+
printf("\n");
503+
printf("Roots: ");
504+
505+
/* Chien search: find roots of the error location polynomial */
506+
for (i = 1; i <= l[u]; i++)
507+
reg[i] = elp[u][i];
508+
count = 0;
509+
for (i = 1; i <= n; i++) {
510+
q = 1;
511+
for (j = 1; j <= l[u]; j++)
512+
if (reg[j] != -1) {
513+
reg[j] = (reg[j] + j) % n;
514+
q ^= alpha_to[reg[j]];
515+
}
516+
if (!q) { /* store root and error
517+
* location number indices */
518+
root[count] = i;
519+
loc[count] = n - i;
520+
count++;
521+
printf("%3d ", n - i);
522+
}
523+
}
524+
printf("\n");
525+
if (count == l[u])
526+
/* no. roots = degree of elp hence <= t errors */
527+
for (i = 0; i < l[u]; i++)
528+
recd[loc[i]] ^= 1;
529+
else { /* elp has degree >t hence cannot solve */
530+
printf("Incomplete decoding: errors detected\n");
531+
return 1;
532+
}
533+
}
534+
}
535+
return 0;
536+
}
537+
538+
void setup(int idx) {
539+
uint64_t pkt = packets[idx];
540+
printf("\n\nPACKET %llx\n", pkt);
541+
printf("-------------------------\n");
542+
pkt >>= 1; // Lose dummy bit */
543+
544+
/* Move BCH code over */
545+
for (int i = length - k - 1; i >= 0; i--) {
546+
recd[i] = pkt & 0x01;
547+
pkt >>= 1;
548+
}
549+
550+
/* Move data over */
551+
for (int i = length - 1; i >= length - k; i--) {
552+
recd[i] = pkt & 0x01;
553+
data[i - length + k] = recd[i];
554+
pkt >>= 1;
555+
}
556+
}
557+
558+
int main()
559+
{
560+
int i;
561+
562+
read_p(); /* Read m */
563+
generate_gf(); /* Construct the Galois Field GF(2**m) */
564+
gen_poly(); /* Compute the generator polynomial of BCH code */
565+
566+
for (int count = 0; count < sizeof(packets) / sizeof(uint64_t); count++) {
567+
setup(count);
568+
memcpy(&orig_recd, recd, sizeof(recd));
569+
570+
int result = decode_bch(); /* DECODE received codeword recd[] */
571+
572+
printf("n=%d, k=%d, length=%d\n", n, k, length);
573+
printf("Results:\n");
574+
printf("\noriginal pkt: ");
575+
for (int i = 0; i < length; i++) {
576+
if (i == length - k) printf(" ");
577+
printf("%d", orig_recd[i]);
578+
}
579+
printf("\n");
580+
581+
if (result) continue;
582+
/*
583+
* DECODING ERRORS? we compare only the data portion
584+
*/
585+
for (i = length - k; i < length; i++)
586+
if (data[i - length + k] != recd[i])
587+
decerror++;
588+
if (decerror) {
589+
printf("There were %d decoding errors in message positions\n", decerror);
590+
continue;
591+
}
592+
593+
/*
594+
* print out original and decoded data
595+
*/
596+
printf("recovered pkt: ");
597+
for (i = 0; i < length; i++) {
598+
if (i == length - k) printf(" ");
599+
printf("%1d", recd[i]);
600+
}
601+
printf("\n");
602+
603+
604+
int flag = 0;
605+
printf("------------ ");
606+
for (int jj = 0; jj < length; jj++) {
607+
if (jj == length - k) printf(" ");
608+
if (orig_recd[jj] != recd[jj]) {
609+
printf("^");
610+
//printf("bit %d: expected %d calc: %d\n", jj, orig_recd[jj], recd[jj]);
611+
flag++;
612+
} else {
613+
printf(" ");
614+
}
615+
}
616+
printf("\n%d ERRORS.\n", flag);
617+
618+
619+
}
620+
}

0 commit comments

Comments
 (0)
Please sign in to comment.