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
Generated on Wed Jul 27 2011 16:27:04 for IT++ by Doxygen 1.7.4