IT++ Logo
bch.cpp
Go to the documentation of this file.
00001 
00029 #include <itpp/comm/bch.h>
00030 #include <itpp/base/binary.h>
00031 #include <itpp/base/specmat.h>
00032 #include <itpp/base/array.h>
00033 
00034 namespace itpp
00035 {
00036 
00037 //---------------------- BCH -----------------------------------
00038 
00039 BCH::BCH(int in_n, int in_k, int in_t, const ivec &genpolynom, bool sys):
00040     n(in_n), k(in_k), t(in_t), systematic(sys)
00041 {
00042   //fix the generator polynomial g(x).
00043   ivec exponents = zeros_i(n - k + 1);
00044   bvec temp = oct2bin(genpolynom, 0);
00045   for (int i = 0; i < temp.length(); i++) {
00046     exponents(i) = static_cast<int>(temp(temp.length() - i - 1)) - 1;
00047   }
00048   g.set(n + 1, exponents);
00049 }
00050 
00051 BCH::BCH(int in_n, int in_t, bool sys):
00052     n(in_n), t(in_t), systematic(sys)
00053 {
00054   // step 1: determine cyclotomic cosets
00055   // although we use elements in GF(n+1), we do not use GFX class, but ivec,
00056   // since we have to multiply by 2 and need the exponents in clear notation
00057   int m_tmp = int2bits(n);
00058   int two_pow_m = 1 << m_tmp;
00059 
00060   it_assert(two_pow_m == n + 1, "BCH::BCH(): (in_n + 1) is not a power of 2");
00061   it_assert((t > 0) && (2*t < n), "BCH::BCH(): in_t must be positive and smaller than n/2");
00062 
00063   Array<ivec> cyclo_sets(2*t + 1);
00064   // unfortunately it is not obvious how many cyclotomic cosets exist (?)
00065   // a bad guess is n/2, which can be a whole lot...
00066   // but we only need 2*t + 1 at maximum for step 2.
00067   // (since all elements are sorted ascending [cp. comment at 2.], the last
00068   // coset we need is the one with coset leader 2t. + coset {0})
00069 
00070   // start with {0} as first set
00071   int curr_coset_idx = 0;
00072   cyclo_sets(curr_coset_idx) = zeros_i(1);
00073 
00074   int cycl_element = 1;
00075 
00076   do {
00077     bool found = false;
00078     // find next element, which is not in a previous coset
00079     do {
00080       int i = 0;
00081       // we do not have to search the first coset, since this is always {0}
00082       found = false;
00083       while ((!found) && (i <= curr_coset_idx)) {
00084         int j = 0;
00085         while ((!found) && (j < cyclo_sets(i).length())) {
00086           if (cycl_element == cyclo_sets(i)(j)) {
00087             found = true;
00088           }
00089           j++;
00090         }
00091         i++;
00092       }
00093       cycl_element++;
00094     }
00095     while ((found) && (cycl_element <= 2*t));
00096 
00097     if (!found) {
00098       // found one
00099       cyclo_sets(++curr_coset_idx).set_size(m_tmp);
00100       // a first guess (we delete afterwards back to correct length):
00101       // there should be no more than m elements in one coset
00102 
00103       int element_index = 0;
00104       cyclo_sets(curr_coset_idx)(element_index) = cycl_element - 1;
00105 
00106       // multiply by two (mod 2^m - 1) as long as new elements are created
00107       while ((((cyclo_sets(curr_coset_idx)(element_index) * 2) % n)
00108               != cyclo_sets(curr_coset_idx)(0))
00109              && (element_index < m_tmp - 1)) {
00110         element_index++;
00111         cyclo_sets(curr_coset_idx)(element_index)
00112           = (cyclo_sets(curr_coset_idx)(element_index - 1) * 2) % n;
00113       }
00114       // delete unused digits
00115       if (element_index + 1 < m_tmp - 1) {
00116         cyclo_sets(curr_coset_idx).del(element_index + 1, m_tmp - 1);
00117       }
00118     }
00119   }
00120   while ((cycl_element <= 2*t) && (curr_coset_idx <= 2*t));
00121 
00122   // step 2: find all cosets that contain all the powers (1..2t) of alpha
00123   // this is pretty easy, since the cosets are in ascending order
00124   // (if regarding the first (=primitive) element for ordering) -
00125   // all due to the method, they have been constructed
00126   // Since we only calculated all cosets up to 2t, this is even trivial
00127   // => we take all curr_coset_idx Cosets
00128 
00129   // maximum index of cosets to be considered
00130   int max_coset_index = curr_coset_idx;
00131 
00132   // step 3: multiply the minimal polynomials corresponding to this sets
00133   // of powers
00134   g.set(two_pow_m, ivec("0"));   // = alpha^0 = 1
00135   ivec min_poly_exp(2);
00136   min_poly_exp(1) = 0;        // product of (x-alpha^cycl_element)
00137 
00138   for (int i = 1; i <= max_coset_index; i++) {
00139     for (int j = 0; j < cyclo_sets(i).length(); j++) {
00140       min_poly_exp(0) = cyclo_sets(i)(j);
00141       g *= GFX(two_pow_m, min_poly_exp);
00142     }
00143   }
00144 
00145   // finally determine k
00146   k = n - g.get_true_degree();
00147 }
00148 
00149 
00150 void BCH::encode(const bvec &uncoded_bits, bvec &coded_bits)
00151 {
00152   int i, j, degree;
00153   int iterations = floor_i(static_cast<double>(uncoded_bits.length()) / k);
00154   GFX m(n + 1, k);
00155   GFX c(n + 1, n);
00156   GFX r(n + 1, n - k);
00157   GFX uncoded_shifted(n + 1, n);
00158   coded_bits.set_size(iterations*n, false);
00159   bvec mbit(k), cbit(n);
00160 
00161   if (systematic)
00162     for (i = 0; i < n - k; i++)
00163       uncoded_shifted[i] = GF(n + 1, -1);
00164 
00165   for (i = 0; i < iterations; i++) {
00166     //Fix the message polynom m(x).
00167     mbit = uncoded_bits.mid(i * k, k);
00168     for (j = 0; j < k; j++) {
00169       degree = static_cast<int>(mbit(k - j - 1)) - 1;
00170       // all bits should be mapped first bit <-> highest degree,
00171       // e.g. 1100 <-> m(x)=x^3 + x^2, but GFX indexes identically
00172       // with exponent m[0] <-> coefficient of x^0
00173       m[j] = GF(n + 1, degree);
00174       if (systematic) {
00175         uncoded_shifted[j+n-k] = m[j];
00176       }
00177     }
00178     //Fix the outputbits cbit.
00179     if (systematic) {
00180       r = modgfx(uncoded_shifted, g);
00181       c = uncoded_shifted - r;
00182       // uncoded_shifted has coefficients from x^(n-k)..x^(n-1)
00183       // and r has coefficients from x^0..x^(n-k-1).
00184       // Thus, this sum perfectly fills c.
00185     }
00186     else {
00187       c = g * m;
00188     }
00189     for (j = 0; j < n; j++) {
00190       if (c[j] == GF(n + 1, 0)) {
00191         cbit(n - j - 1) = 1;  // again reverse mapping like mbit(.)
00192       }
00193       else {
00194         cbit(n - j - 1) = 0;
00195       }
00196     }
00197     coded_bits.replace_mid(i*n, cbit);
00198   }
00199 }
00200 
00201 bvec BCH::encode(const bvec &uncoded_bits)
00202 {
00203   bvec coded_bits;
00204   encode(uncoded_bits, coded_bits);
00205   return coded_bits;
00206 }
00207 
00208 void BCH::decode(const bvec &coded_bits, bvec &decoded_bits)
00209 {
00210   int j, i, degree, kk, foundzeros, cisvalid;
00211   int iterations = floor_i(static_cast<double>(coded_bits.length()) / n);
00212   bvec rbin(n), mbin(k);
00213   decoded_bits.set_size(iterations*k, false);
00214 
00215   GFX r(n + 1, n - 1), c(n + 1, n - 1), m(n + 1, k - 1), S(n + 1, 2*t), Lambda(n + 1),
00216   OldLambda(n + 1), T(n + 1), Ohmega(n + 1), One(n + 1, (char*)"0");
00217   GF delta(n + 1);
00218   ivec errorpos;
00219 
00220   for (i = 0; i < iterations; i++) {
00221     //Fix the received polynomial r(x)
00222     rbin = coded_bits.mid(i * n, n);
00223     for (j = 0; j < n; j++) {
00224       degree = static_cast<int>(rbin(n - j - 1)) - 1;
00225       // reverse mapping, see encode(.)
00226       r[j] = GF(n + 1, degree);
00227     }
00228     //Fix the syndrome polynomial S(x).
00229     S[0] = GF(n + 1, -1);
00230     for (j = 1; j <= 2*t; j++) {
00231       S[j] =  r(GF(n + 1, j));
00232     }
00233     if (S.get_true_degree() >= 1) { //Errors in the received word
00234       //Iterate to find Lambda(x).
00235       kk = 0;
00236       Lambda = GFX(n + 1, (char*)"0");
00237       T = GFX(n + 1, (char*)"0");
00238       while (kk < t) {
00239         Ohmega = Lambda * (S + One);
00240         delta = Ohmega[2*kk+1];
00241         OldLambda = Lambda;
00242         Lambda = OldLambda + delta * (GFX(n + 1, (char*)"-1 0") * T);
00243         if ((delta == GF(n + 1, -1)) || (OldLambda.get_true_degree() > kk)) {
00244           T = GFX(n + 1, (char*)"-1 -1 0") * T;
00245         }
00246         else {
00247           T = (GFX(n + 1, (char*)"-1 0") * OldLambda) / delta;
00248         }
00249         kk = kk + 1;
00250       }
00251       //Find the zeros to Lambda(x).
00252       errorpos.set_size(Lambda.get_true_degree(), true);
00253       foundzeros = 0;
00254       for (j = 0; j <= n - 1; j++) {
00255         if (Lambda(GF(n + 1, j)) == GF(n + 1, -1)) {
00256           errorpos(foundzeros) = (n - j) % n;
00257           foundzeros += 1;
00258           if (foundzeros >= Lambda.get_true_degree()) {
00259             break;
00260           }
00261         }
00262       }
00263       //Correct the codeword.
00264       for (j = 0; j < foundzeros; j++) {
00265         rbin(n - errorpos(j) - 1) += 1;  // again, reverse mapping
00266       }
00267       //Reconstruct the corrected codeword.
00268       for (j = 0; j < n; j++) {
00269         degree = static_cast<int>(rbin(n - j - 1)) - 1;
00270         c[j] = GF(n + 1, degree);
00271       }
00272       //Code word validation.
00273       S[0] = GF(n + 1, -1);
00274       for (j = 1; j <= 2*t; j++) {
00275         S[j] =  c(GF(n + 1, j));
00276       }
00277       if (S.get_true_degree() <= 0) { //c(x) is a valid codeword.
00278         cisvalid = true;
00279       }
00280       else {
00281         cisvalid = false;
00282       }
00283     }
00284     else {
00285       c = r;
00286       cisvalid = true;
00287     }
00288     //Construct the message bit vector.
00289     if (cisvalid) { //c(x) is a valid codeword.
00290       if (c.get_true_degree() > 1) {
00291         if (systematic) {
00292           for (j = 0; j < k; j++)
00293             m[j] = c[n-k+j];
00294         }
00295         else {
00296           m = divgfx(c, g);
00297         }
00298         mbin.clear();
00299         for (j = 0; j <= m.get_true_degree(); j++) {
00300           if (m[j] == GF(n + 1, 0)) {
00301             mbin(k - j - 1) = 1;
00302           }
00303         }
00304       }
00305       else { //The zero word was transmitted
00306         mbin = zeros_b(k);
00307         m = GFX(n + 1, (char*)"-1");
00308       }
00309     }
00310     else { //Decoder failure.
00311       mbin = zeros_b(k);
00312       m = GFX(n + 1, (char*)"-1");
00313     }
00314     decoded_bits.replace_mid(i*k, mbin);
00315   }
00316 }
00317 
00318 
00319 bvec BCH::decode(const bvec &coded_bits)
00320 {
00321   bvec decoded_bits;
00322   decode(coded_bits, decoded_bits);
00323   return decoded_bits;
00324 }
00325 
00326 
00327 // --- Soft-decision decoding is not implemented ---
00328 
00329 void BCH::decode(const vec &, bvec &)
00330 {
00331   it_error("BCH::decode(): Soft-decision decoding is not implemented");
00332 }
00333 
00334 bvec BCH::decode(const vec &)
00335 {
00336   it_error("BCH::decode(): Soft-decision decoding is not implemented");
00337   return bvec();
00338 }
00339 
00340 
00341 } // namespace itpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
SourceForge Logo

Generated on Sat Jul 9 2011 15:21:31 for IT++ by Doxygen 1.7.4