IT++ Logo
bch.cpp
Go to the documentation of this file.
1 
29 #include <itpp/comm/bch.h>
30 #include <itpp/base/binary.h>
31 #include <itpp/base/specmat.h>
32 #include <itpp/base/array.h>
33 
34 namespace itpp
35 {
36 
37 //---------------------- BCH -----------------------------------
38 
39 BCH::BCH(int in_n, int in_k, int in_t, const ivec &genpolynom, bool sys) :
40  n(in_n), k(in_k), t(in_t), systematic(sys)
41 {
42  //fix the generator polynomial g(x).
43  ivec exponents = zeros_i(n - k + 1);
44  bvec temp = oct2bin(genpolynom, 0);
45  for (int i = 0; i < temp.length(); i++) {
46  exponents(i) = static_cast<int>(temp(temp.length() - i - 1)) - 1;
47  }
48  g.set(n + 1, exponents);
49 }
50 
51 BCH::BCH(int in_n, int in_t, bool sys) :
52  n(in_n), t(in_t), systematic(sys)
53 {
54  // step 1: determine cyclotomic cosets
55  // although we use elements in GF(n+1), we do not use GFX class, but ivec,
56  // since we have to multiply by 2 and need the exponents in clear notation
57  int m_tmp = int2bits(n);
58  int two_pow_m = 1 << m_tmp;
59 
60  it_assert(two_pow_m == n + 1, "BCH::BCH(): (in_n + 1) is not a power of 2");
61  it_assert((t > 0) && (2 * t < n), "BCH::BCH(): in_t must be positive and smaller than n/2");
62 
63  Array<ivec> cyclo_sets(2 * t + 1);
64  // unfortunately it is not obvious how many cyclotomic cosets exist (?)
65  // a bad guess is n/2, which can be a whole lot...
66  // but we only need 2*t + 1 at maximum for step 2.
67  // (since all elements are sorted ascending [cp. comment at 2.], the last
68  // coset we need is the one with coset leader 2t. + coset {0})
69 
70  // start with {0} as first set
71  int curr_coset_idx = 0;
72  cyclo_sets(curr_coset_idx) = zeros_i(1);
73 
74  int cycl_element = 1;
75 
76  do {
77  bool found = false;
78  // find next element, which is not in a previous coset
79  do {
80  int i = 0;
81  // we do not have to search the first coset, since this is always {0}
82  found = false;
83  while ((!found) && (i <= curr_coset_idx)) {
84  int j = 0;
85  while ((!found) && (j < cyclo_sets(i).length())) {
86  if (cycl_element == cyclo_sets(i)(j)) {
87  found = true;
88  }
89  j++;
90  }
91  i++;
92  }
93  cycl_element++;
94  }
95  while ((found) && (cycl_element <= 2 * t));
96 
97  if (!found) {
98  // found one
99  cyclo_sets(++curr_coset_idx).set_size(m_tmp);
100  // a first guess (we delete afterwards back to correct length):
101  // there should be no more than m elements in one coset
102 
103  int element_index = 0;
104  cyclo_sets(curr_coset_idx)(element_index) = cycl_element - 1;
105 
106  // multiply by two (mod 2^m - 1) as long as new elements are created
107  while ((((cyclo_sets(curr_coset_idx)(element_index) * 2) % n)
108  != cyclo_sets(curr_coset_idx)(0))
109  && (element_index < m_tmp - 1)) {
110  element_index++;
111  cyclo_sets(curr_coset_idx)(element_index)
112  = (cyclo_sets(curr_coset_idx)(element_index - 1) * 2) % n;
113  }
114  // delete unused digits
115  if (element_index + 1 < m_tmp - 1) {
116  cyclo_sets(curr_coset_idx).del(element_index + 1, m_tmp - 1);
117  }
118  }
119  }
120  while ((cycl_element <= 2 * t) && (curr_coset_idx <= 2 * t));
121 
122  // step 2: find all cosets that contain all the powers (1..2t) of alpha
123  // this is pretty easy, since the cosets are in ascending order
124  // (if regarding the first (=primitive) element for ordering) -
125  // all due to the method, they have been constructed
126  // Since we only calculated all cosets up to 2t, this is even trivial
127  // => we take all curr_coset_idx Cosets
128 
129  // maximum index of cosets to be considered
130  int max_coset_index = curr_coset_idx;
131 
132  // step 3: multiply the minimal polynomials corresponding to this sets
133  // of powers
134  g.set(two_pow_m, ivec("0")); // = alpha^0 = 1
135  ivec min_poly_exp(2);
136  min_poly_exp(1) = 0; // product of (x-alpha^cycl_element)
137 
138  for (int i = 1; i <= max_coset_index; i++) {
139  for (int j = 0; j < cyclo_sets(i).length(); j++) {
140  min_poly_exp(0) = cyclo_sets(i)(j);
141  g *= GFX(two_pow_m, min_poly_exp);
142  }
143  }
144 
145  // finally determine k
146  k = n - g.get_true_degree();
147 }
148 
149 
150 void BCH::encode(const bvec &uncoded_bits, bvec &coded_bits)
151 {
152  int i, j, degree;
153  int iterations = floor_i(static_cast<double>(uncoded_bits.length()) / k);
154  GFX m(n + 1, k);
155  GFX c(n + 1, n);
156  GFX r(n + 1, n - k);
157  GFX uncoded_shifted(n + 1, n);
158  coded_bits.set_size(iterations * n, false);
159  bvec mbit(k), cbit(n);
160 
161  if (systematic)
162  for (i = 0; i < n - k; i++)
163  uncoded_shifted[i] = GF(n + 1, -1);
164 
165  for (i = 0; i < iterations; i++) {
166  //Fix the message polynom m(x).
167  mbit = uncoded_bits.mid(i * k, k);
168  for (j = 0; j < k; j++) {
169  degree = static_cast<int>(mbit(k - j - 1)) - 1;
170  // all bits should be mapped first bit <-> highest degree,
171  // e.g. 1100 <-> m(x)=x^3 + x^2, but GFX indexes identically
172  // with exponent m[0] <-> coefficient of x^0
173  m[j] = GF(n + 1, degree);
174  if (systematic) {
175  uncoded_shifted[j + n - k] = m[j];
176  }
177  }
178  //Fix the outputbits cbit.
179  if (systematic) {
180  r = modgfx(uncoded_shifted, g);
181  c = uncoded_shifted - r;
182  // uncoded_shifted has coefficients from x^(n-k)..x^(n-1)
183  // and r has coefficients from x^0..x^(n-k-1).
184  // Thus, this sum perfectly fills c.
185  }
186  else {
187  c = g * m;
188  }
189  for (j = 0; j < n; j++) {
190  if (c[j] == GF(n + 1, 0)) {
191  cbit(n - j - 1) = 1; // again reverse mapping like mbit(.)
192  }
193  else {
194  cbit(n - j - 1) = 0;
195  }
196  }
197  coded_bits.replace_mid(i * n, cbit);
198  }
199 }
200 
201 bvec BCH::encode(const bvec &uncoded_bits)
202 {
203  bvec coded_bits;
204  encode(uncoded_bits, coded_bits);
205  return coded_bits;
206 }
207 
208 bool BCH::decode(const bvec &coded_bits, bvec &decoded_message, bvec &cw_isvalid)
209 {
210  bool decoderfailure, no_dec_failure;
211  int j, i, degree, kk, foundzeros;
212  ivec errorpos;
213  int iterations = floor_i(static_cast<double>(coded_bits.length()) / n);
214  bvec rbin(n), mbin(k);
215  decoded_message.set_size(iterations * k, false);
216  cw_isvalid.set_length(iterations);
217 
218  GFX r(n + 1, n - 1), c(n + 1, n - 1), m(n + 1, k - 1), S(n + 1, 2 * t), Lambda(n + 1),
219  OldLambda(n + 1), T(n + 1), Omega(n + 1), One(n + 1, (char*) "0");
220  GF delta(n + 1);
221 
222  no_dec_failure = true;
223  for (i = 0; i < iterations; i++) {
224  decoderfailure = false;
225  //Fix the received polynomial r(x)
226  rbin = coded_bits.mid(i * n, n);
227  for (j = 0; j < n; j++) {
228  degree = static_cast<int>(rbin(n - j - 1)) - 1;
229  // reverse mapping, see encode(.)
230  r[j] = GF(n + 1, degree);
231  }
232  //Fix the syndrome polynomial S(x).
233  S.clear();
234  for (j = 1; j <= 2 * t; j++) {
235  S[j] = r(GF(n + 1, j));
236  }
237  if (S.get_true_degree() >= 1) { //Errors in the received word
238  //Iterate to find Lambda(x).
239  kk = 0;
240  Lambda = GFX(n + 1, (char*) "0");
241  T = GFX(n + 1, (char*) "0");
242  while (kk < t) {
243  Omega = Lambda * (S + One);
244  delta = Omega[2 * kk + 1];
245  OldLambda = Lambda;
246  Lambda = OldLambda + delta * (GFX(n + 1, (char*) "-1 0") * T);
247  if ((delta == GF(n + 1, -1)) || (OldLambda.get_true_degree() > kk)) {
248  T = GFX(n + 1, (char*) "-1 -1 0") * T;
249  }
250  else {
251  T = (GFX(n + 1, (char*) "-1 0") * OldLambda) / delta;
252  }
253  kk = kk + 1;
254  }
255  //Find the zeros to Lambda(x).
256  errorpos.set_size(Lambda.get_true_degree());
257  foundzeros = 0;
258  for (j = 0; j <= n - 1; j++) {
259  if (Lambda(GF(n + 1, j)) == GF(n + 1, -1)) {
260  errorpos(foundzeros) = (n - j) % n;
261  foundzeros += 1;
262  if (foundzeros >= Lambda.get_true_degree()) {
263  break;
264  }
265  }
266  }
267  if (foundzeros != Lambda.get_true_degree()) {
268  decoderfailure = true;
269  }
270  else {
271  //Correct the codeword.
272  for (j = 0; j < foundzeros; j++) {
273  rbin(n - errorpos(j) - 1) += 1; // again, reverse mapping
274  }
275  //Reconstruct the corrected codeword.
276  for (j = 0; j < n; j++) {
277  degree = static_cast<int>(rbin(n - j - 1)) - 1;
278  c[j] = GF(n + 1, degree);
279  }
280  //Code word validation.
281  S.clear();
282  for (j = 1; j <= 2 * t; j++) {
283  S[j] = c(GF(n + 1, j));
284  }
285  if (S.get_true_degree() <= 0) { //c(x) is a valid codeword.
286  decoderfailure = false;
287  }
288  else {
289  decoderfailure = true;
290  }
291  }
292  }
293  else {
294  c = r;
295  decoderfailure = false;
296  }
297  //Construct the message bit vector.
298  if (decoderfailure == false) { //c(x) is a valid codeword.
299  if (c.get_true_degree() > 1) {
300  if (systematic) {
301  for (j = 0; j < k; j++)
302  m[j] = c[n - k + j];
303  }
304  else {
305  m = divgfx(c, g);
306  }
307  mbin.clear();
308  for (j = 0; j <= m.get_true_degree(); j++) {
309  if (m[j] == GF(n + 1, 0)) {
310  mbin(k - j - 1) = 1;
311  }
312  }
313  }
314  else { //The zero word was transmitted
315  mbin = zeros_b(k);
316  m = GFX(n + 1, (char*) "-1");
317  }
318  }
319  else { //Decoder failure.
320  // for a systematic code it is better to extract the undecoded message
321  // from the received code word, i.e. obtaining a bit error
322  // prob. p_b << 1/2, than setting all-zero (p_b = 1/2)
323  if (systematic) {
324  mbin = coded_bits.mid(i * n, k);
325  }
326  else {
327  mbin = zeros_b(k);
328  }
329  no_dec_failure = false;
330  }
331  decoded_message.replace_mid(i * k, mbin);
332  cw_isvalid(i) = (!decoderfailure);
333  }
334  return no_dec_failure;
335 }
336 
337 void BCH::decode(const bvec &coded_bits, bvec &decoded_bits)
338 {
339  bvec cw_isvalid;
340  if (!decode(coded_bits, decoded_bits, cw_isvalid)) {
341  for (int i = 0; i < cw_isvalid.length(); i++) {
342  if (!cw_isvalid(i)) {
343  decoded_bits.replace_mid(i * k, zeros_b(k));
344  }
345  }
346  }
347 }
348 
349 bvec BCH::decode(const bvec &coded_bits)
350 {
351  bvec decoded_bits;
352  decode(coded_bits, decoded_bits);
353  return decoded_bits;
354 }
355 
356 
357 // --- Soft-decision decoding is not implemented ---
358 
359 void BCH::decode(const vec &, bvec &)
360 {
361  it_error("BCH::decode(): Soft-decision decoding is not implemented");
362 }
363 
364 bvec BCH::decode(const vec &)
365 {
366  it_error("BCH::decode(): Soft-decision decoding is not implemented");
367  return bvec();
368 }
369 
370 
371 } // namespace itpp
SourceForge Logo

Generated on Sat May 25 2013 16:32:20 for IT++ by Doxygen 1.8.2