IT++ Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 void BCH::decode(const bvec &coded_bits, bvec &decoded_bits)
209 {
210  int j, i, degree, kk, foundzeros, cisvalid;
211  int iterations = floor_i(static_cast<double>(coded_bits.length()) / n);
212  bvec rbin(n), mbin(k);
213  decoded_bits.set_size(iterations*k, false);
214 
215  GFX r(n + 1, n - 1), c(n + 1, n - 1), m(n + 1, k - 1), S(n + 1, 2*t), Lambda(n + 1),
216  OldLambda(n + 1), T(n + 1), Ohmega(n + 1), One(n + 1, (char*)"0");
217  GF delta(n + 1);
218  ivec errorpos;
219 
220  for (i = 0; i < iterations; i++) {
221  //Fix the received polynomial r(x)
222  rbin = coded_bits.mid(i * n, n);
223  for (j = 0; j < n; j++) {
224  degree = static_cast<int>(rbin(n - j - 1)) - 1;
225  // reverse mapping, see encode(.)
226  r[j] = GF(n + 1, degree);
227  }
228  //Fix the syndrome polynomial S(x).
229  S[0] = GF(n + 1, -1);
230  for (j = 1; j <= 2*t; j++) {
231  S[j] = r(GF(n + 1, j));
232  }
233  if (S.get_true_degree() >= 1) { //Errors in the received word
234  //Iterate to find Lambda(x).
235  kk = 0;
236  Lambda = GFX(n + 1, (char*)"0");
237  T = GFX(n + 1, (char*)"0");
238  while (kk < t) {
239  Ohmega = Lambda * (S + One);
240  delta = Ohmega[2*kk+1];
241  OldLambda = Lambda;
242  Lambda = OldLambda + delta * (GFX(n + 1, (char*)"-1 0") * T);
243  if ((delta == GF(n + 1, -1)) || (OldLambda.get_true_degree() > kk)) {
244  T = GFX(n + 1, (char*)"-1 -1 0") * T;
245  }
246  else {
247  T = (GFX(n + 1, (char*)"-1 0") * OldLambda) / delta;
248  }
249  kk = kk + 1;
250  }
251  //Find the zeros to Lambda(x).
252  errorpos.set_size(Lambda.get_true_degree(), true);
253  foundzeros = 0;
254  for (j = 0; j <= n - 1; j++) {
255  if (Lambda(GF(n + 1, j)) == GF(n + 1, -1)) {
256  errorpos(foundzeros) = (n - j) % n;
257  foundzeros += 1;
258  if (foundzeros >= Lambda.get_true_degree()) {
259  break;
260  }
261  }
262  }
263  //Correct the codeword.
264  for (j = 0; j < foundzeros; j++) {
265  rbin(n - errorpos(j) - 1) += 1; // again, reverse mapping
266  }
267  //Reconstruct the corrected codeword.
268  for (j = 0; j < n; j++) {
269  degree = static_cast<int>(rbin(n - j - 1)) - 1;
270  c[j] = GF(n + 1, degree);
271  }
272  //Code word validation.
273  S[0] = GF(n + 1, -1);
274  for (j = 1; j <= 2*t; j++) {
275  S[j] = c(GF(n + 1, j));
276  }
277  if (S.get_true_degree() <= 0) { //c(x) is a valid codeword.
278  cisvalid = true;
279  }
280  else {
281  cisvalid = false;
282  }
283  }
284  else {
285  c = r;
286  cisvalid = true;
287  }
288  //Construct the message bit vector.
289  if (cisvalid) { //c(x) is a valid codeword.
290  if (c.get_true_degree() > 1) {
291  if (systematic) {
292  for (j = 0; j < k; j++)
293  m[j] = c[n-k+j];
294  }
295  else {
296  m = divgfx(c, g);
297  }
298  mbin.clear();
299  for (j = 0; j <= m.get_true_degree(); j++) {
300  if (m[j] == GF(n + 1, 0)) {
301  mbin(k - j - 1) = 1;
302  }
303  }
304  }
305  else { //The zero word was transmitted
306  mbin = zeros_b(k);
307  m = GFX(n + 1, (char*)"-1");
308  }
309  }
310  else { //Decoder failure.
311  mbin = zeros_b(k);
312  m = GFX(n + 1, (char*)"-1");
313  }
314  decoded_bits.replace_mid(i*k, mbin);
315  }
316 }
317 
318 
319 bvec BCH::decode(const bvec &coded_bits)
320 {
321  bvec decoded_bits;
322  decode(coded_bits, decoded_bits);
323  return decoded_bits;
324 }
325 
326 
327 // --- Soft-decision decoding is not implemented ---
328 
329 void BCH::decode(const vec &, bvec &)
330 {
331  it_error("BCH::decode(): Soft-decision decoding is not implemented");
332 }
333 
334 bvec BCH::decode(const vec &)
335 {
336  it_error("BCH::decode(): Soft-decision decoding is not implemented");
337  return bvec();
338 }
339 
340 
341 } // namespace itpp
SourceForge Logo

Generated on Fri Mar 21 2014 17:14:13 for IT++ by Doxygen 1.8.1.2