00001 00029 #include <itpp/comm/ldpc.h> 00030 #include <iomanip> 00031 #include <sstream> 00032 00033 namespace itpp 00034 { 00035 00042 static const int LDPC_binary_file_version = 2; 00043 00044 // --------------------------------------------------------------------------- 00045 // LDPC_Parity 00046 // --------------------------------------------------------------------------- 00047 00048 // public methods 00049 00050 LDPC_Parity::LDPC_Parity(int nc, int nv): init_flag(false) 00051 { 00052 initialize(nc, nv); 00053 } 00054 00055 LDPC_Parity::LDPC_Parity(const std::string& filename, 00056 const std::string& format): init_flag(false) 00057 { 00058 if (format == "alist") { 00059 load_alist(filename); 00060 } 00061 else { 00062 it_error("LDPC_Parity::LDPC_Parity(): Only 'alist' format is supported"); 00063 } 00064 } 00065 00066 LDPC_Parity::LDPC_Parity(const GF2mat_sparse_alist &alist): 00067 init_flag(false) 00068 { 00069 import_alist(alist); 00070 } 00071 00072 void LDPC_Parity::initialize(int nc, int nv) 00073 { 00074 ncheck = nc; 00075 nvar = nv; 00076 H = GF2mat_sparse(ncheck, nvar); 00077 Ht = GF2mat_sparse(nvar, ncheck); 00078 sumX1 = zeros_i(nvar); 00079 sumX2 = zeros_i(ncheck); 00080 init_flag = true; 00081 } 00082 00083 void LDPC_Parity::set(int i, int j, bin x) 00084 { 00085 it_assert(init_flag, "LDPC_Parity::set(): Object not initialized"); 00086 it_assert_debug((i >= 0) && (i < ncheck), 00087 "LDPC_Parity::set(): Wrong index i"); 00088 it_assert_debug((j >= 0) && (j < nvar), 00089 "LDPC_Parity::set(): Wrong index j"); 00090 it_assert_debug(H(i, j) == Ht(j, i), "LDPC_Parity:set(): Internal error"); 00091 00092 int diff = static_cast<int>(x) - static_cast<int>(H(i, j)); 00093 sumX1(j) += diff; 00094 sumX2(i) += diff; 00095 00096 if (x == 1) { 00097 H.set(i, j, 1); 00098 Ht.set(j, i, 1); 00099 } 00100 else { 00101 H.clear_elem(i, j); 00102 Ht.clear_elem(j, i); 00103 } 00104 00105 it_assert_debug(H(i, j) == x, "LDPC_Parity::set(): Internal error"); 00106 it_assert_debug(Ht(j, i) == x, "LDPC_Parity::set(): Internal error"); 00107 } 00108 00109 void LDPC_Parity::display_stats() const 00110 { 00111 it_assert(init_flag, 00112 "LDPC_Parity::display_stats(): Object not initialized"); 00113 int cmax = max(sumX1); 00114 int vmax = max(sumX2); 00115 vec vdeg = zeros(cmax + 1); // number of variable nodes with n neighbours 00116 vec cdeg = zeros(vmax + 1); // number of check nodes with n neighbours 00117 for (int col = 0; col < nvar; col++) { 00118 vdeg(length(get_col(col).get_nz_indices()))++; 00119 it_assert(sumX1(col) == length(get_col(col).get_nz_indices()), 00120 "LDPC_Parity::display_stats(): Internal error"); 00121 } 00122 for (int row = 0; row < ncheck; row++) { 00123 cdeg(length(get_row(row).get_nz_indices()))++; 00124 it_assert(sumX2(row) == length(get_row(row).get_nz_indices()), 00125 "LDPC_Parity::display_stats(): Internal error"); 00126 } 00127 00128 // from edge perspective 00129 // number of edges connected to vnodes of degree n 00130 vec vdegedge = elem_mult(vdeg, linspace(0, vdeg.length() - 1, 00131 vdeg.length())); 00132 // number of edges connected to cnodes of degree n 00133 vec cdegedge = elem_mult(cdeg, linspace(0, cdeg.length() - 1, 00134 cdeg.length())); 00135 00136 int edges = sum(elem_mult(to_ivec(linspace(0, vdeg.length() - 1, 00137 vdeg.length())), 00138 to_ivec(vdeg))); 00139 00140 it_info("--- LDPC parity check matrix ---"); 00141 it_info("Dimension [ncheck x nvar]: " << ncheck << " x " << nvar); 00142 it_info("Variable node degree distribution from node perspective:\n" 00143 << vdeg / nvar); 00144 it_info("Check node degree distribution from node perspective:\n" 00145 << cdeg / ncheck); 00146 it_info("Variable node degree distribution from edge perspective:\n" 00147 << vdegedge / edges); 00148 it_info("Check node degree distribution from edge perspective:\n" 00149 << cdegedge / edges); 00150 it_info("Rate: " << get_rate()); 00151 it_info("--------------------------------"); 00152 } 00153 00154 00155 void LDPC_Parity::load_alist(const std::string& alist_file) 00156 { 00157 import_alist(GF2mat_sparse_alist(alist_file)); 00158 } 00159 00160 void LDPC_Parity::save_alist(const std::string& alist_file) const 00161 { 00162 GF2mat_sparse_alist alist = export_alist(); 00163 alist.write(alist_file); 00164 } 00165 00166 00167 void LDPC_Parity::import_alist(const GF2mat_sparse_alist& alist) 00168 { 00169 GF2mat_sparse X = alist.to_sparse(); 00170 00171 initialize(X.rows(), X.cols()); 00172 // brute force copy from X to this 00173 for (int i = 0; i < ncheck; i++) { 00174 for (int j = 0; j < nvar; j++) { 00175 if (X(i, j)) { 00176 set(i, j, 1); 00177 } 00178 } 00179 } 00180 } 00181 00182 GF2mat_sparse_alist LDPC_Parity::export_alist() const 00183 { 00184 it_assert(init_flag, 00185 "LDPC_Parity::export_alist(): Object not initialized"); 00186 GF2mat_sparse_alist alist; 00187 alist.from_sparse(H); 00188 return alist; 00189 } 00190 00191 00192 int LDPC_Parity::check_connectivity(int from_i, int from_j, int to_i, 00193 int to_j, int godir, int L) const 00194 { 00195 it_assert(init_flag, 00196 "LDPC_Parity::check_connectivity(): Object not initialized"); 00197 int i, j, result; 00198 00199 if (L < 0) { // unable to reach coordinate with given L 00200 return (-3); 00201 } 00202 00203 // check if reached destination 00204 if ((from_i == to_i) && (from_j == to_j) && (godir != 0)) { 00205 return L; 00206 } 00207 00208 if (get(from_i, from_j) == 0) { // meaningless search 00209 return (-2); 00210 } 00211 00212 if (L == 2) { // Treat this case separately for efficiency 00213 if (godir == 2) { // go horizontally 00214 if (get(from_i, to_j) == 1) { return 0; } 00215 } 00216 if (godir == 1) { // go vertically 00217 if (get(to_i, from_j) == 1) { return 0; } 00218 } 00219 return (-3); 00220 } 00221 00222 if ((godir == 1) || (godir == 0)) { // go vertically 00223 ivec cj = get_col(from_j).get_nz_indices(); 00224 for (i = 0; i < length(cj); i++) { 00225 if (cj(i) != from_i) { 00226 result = check_connectivity(cj(i), from_j, to_i, to_j, 2, L - 1); 00227 if (result >= 0) { 00228 return (result); 00229 } 00230 } 00231 } 00232 } 00233 00234 if (godir == 2) { // go horizontally 00235 ivec ri = get_row(from_i).get_nz_indices(); 00236 for (j = 0; j < length(ri); j++) { 00237 if (ri(j) != from_j) { 00238 result = check_connectivity(from_i, ri(j), to_i, to_j, 1, L - 1); 00239 if (result >= 0) { 00240 return (result); 00241 } 00242 } 00243 } 00244 } 00245 00246 return (-1); 00247 }; 00248 00249 int LDPC_Parity::check_for_cycles(int L) const 00250 { 00251 it_assert(init_flag, 00252 "LDPC_Parity::check_for_cycles(): Object not initialized"); 00253 // looking for odd length cycles does not make sense 00254 if ((L&1) == 1) { return (-1); } 00255 if (L == 0) { return (-4); } 00256 00257 int cycles = 0; 00258 for (int i = 0; i < nvar; i++) { 00259 ivec ri = get_col(i).get_nz_indices(); 00260 for (int j = 0; j < length(ri); j++) { 00261 if (check_connectivity(ri(j), i, ri(j), i, 0, L) >= 0) { 00262 cycles++; 00263 } 00264 } 00265 } 00266 return cycles; 00267 }; 00268 00269 // ivec LDPC_Parity::get_rowdegree() const 00270 // { 00271 // ivec rdeg = zeros_i(Nmax); 00272 // for (int i=0; i<ncheck; i++) { 00273 // rdeg(sumX2(i))++; 00274 // } 00275 // return rdeg; 00276 // }; 00277 00278 // ivec LDPC_Parity::get_coldegree() const 00279 // { 00280 // ivec cdeg = zeros_i(Nmax); 00281 // for (int j=0; j<nvar; j++) { 00282 // cdeg(sumX1(j))++; 00283 // } 00284 // return cdeg; 00285 // }; 00286 00287 00288 // ---------------------------------------------------------------------- 00289 // LDPC_Parity_Unstructured 00290 // ---------------------------------------------------------------------- 00291 00292 int LDPC_Parity_Unstructured::cycle_removal_MGW(int Maxcyc) 00293 { 00294 it_assert(init_flag, 00295 "LDPC_Parity::cycle_removal_MGW(): Object not initialized"); 00296 typedef Sparse_Mat<short> Ssmat; 00297 typedef Sparse_Vec<short> Ssvec; 00298 00299 Maxcyc -= 2; 00300 00301 // Construct the adjacency matrix of the graph 00302 Ssmat G(ncheck + nvar, ncheck + nvar, 5); 00303 for (int j = 0; j < nvar; j++) { 00304 GF2vec_sparse col = get_col(j); 00305 for (int i = 0; i < col.nnz(); i++) { 00306 if (get(col.get_nz_index(i), j) == 1) { 00307 G.set(col.get_nz_index(i), j + ncheck, 1); 00308 G.set(ncheck + j, col.get_nz_index(i), 1); 00309 } 00310 } 00311 } 00312 00313 Array<Ssmat> Gpow(Maxcyc); 00314 Gpow(0).set_size(ncheck + nvar, ncheck + nvar, 1); 00315 Gpow(0).clear(); 00316 for (int i = 0; i < ncheck + nvar; i++) { 00317 Gpow(0).set(i, i, 1); 00318 } 00319 Gpow(1) = G; 00320 00321 /* Main cycle elimination loop starts here. Note that G and all 00322 powers of G are symmetric matrices. This fact is exploited in 00323 the code. 00324 */ 00325 int r; 00326 int cycles_found = 0; 00327 int scl = Maxcyc; 00328 for (r = 4; r <= Maxcyc; r += 2) { 00329 // compute the next power of the adjacency matrix 00330 Gpow(r / 2) = Gpow(r / 2 - 1) * G; 00331 bool traverse_again; 00332 do { 00333 traverse_again = false; 00334 cycles_found = 0; 00335 it_info_debug("Starting new pass of cycle elimination, target girth " 00336 << (r + 2) << "..."); 00337 int pdone = 0; 00338 for (int j = 0; j < ncheck + nvar; j++) { // loop over elements of G 00339 for (int i = 0; i < ncheck + nvar; i++) { 00340 int ptemp = floor_i(100.0 * (i + j * (ncheck + nvar)) / 00341 ((nvar + ncheck) * (nvar + ncheck))); 00342 if (ptemp > pdone + 10) { 00343 it_info_debug(ptemp << "% done."); 00344 pdone = ptemp; 00345 } 00346 00347 if (((Gpow(r / 2))(i, j) >= 2) && ((Gpow(r / 2 - 2))(i, j) == 0)) { 00348 // Found a cycle. 00349 cycles_found++; 00350 00351 // choose k 00352 ivec tmpi = (elem_mult(Gpow(r / 2 - 1).get_col(i), 00353 G.get_col(j))).get_nz_indices(); 00354 // int k = tmpi(rand()%length(tmpi)); 00355 int k = tmpi(randi(0, length(tmpi) - 1)); 00356 it_assert_debug(G(j, k) == 1 && G(k, j) == 1, 00357 "LDPC_Parity_Unstructured::cycle_removal_MGW(): " 00358 "Internal error"); 00359 00360 // determine candidate edges for an edge swap 00361 Ssvec rowjk = Gpow(r / 2) * (Gpow(r / 2 - 1).get_col(j) 00362 + Gpow(r / 2 - 1).get_col(k)); 00363 int p, l; 00364 ivec Ce_ind = sort_index(randu(nvar + ncheck)); // random order 00365 00366 for (int s = 0; s < nvar + ncheck; s++) { 00367 l = Ce_ind(s); 00368 if (rowjk(l) != 0) { continue; } 00369 ivec colcandi = G.get_col(l).get_nz_indices(); 00370 if (length(colcandi) > 0) { 00371 // select a node p which is a member of Ce 00372 for (int u = 0; u < length(colcandi); u++) { 00373 p = colcandi(u); 00374 if (p != l) { 00375 if (rowjk(p) == 0) { 00376 goto found_candidate_vector; 00377 } 00378 } 00379 } 00380 } 00381 } 00382 continue; // go to the next entry (i,j) 00383 00384 found_candidate_vector: 00385 // swap edges 00386 00387 if (p >= ncheck) { int z = l; l = p; p = z; } 00388 if (j >= ncheck) { int z = k; k = j; j = z; } 00389 00390 // Swap endpoints of edges (p,l) and (j,k) 00391 // cout << "(" << j << "," << k << ")<->(" 00392 // << p << "," << l << ") " ; 00393 // cout << "."; 00394 // cout.flush(); 00395 00396 // Update the matrix 00397 it_assert_debug((get(j, k - ncheck) == 1) && (get(p, l - ncheck) == 1), 00398 "LDPC_Parity_Unstructured::cycle_removal_MGW(): " 00399 "Internal error"); 00400 set(j, k - ncheck, 0); 00401 set(p, l - ncheck, 0); 00402 it_assert_debug((get(j, l - ncheck) == 0) && (get(p, k - ncheck) == 0), 00403 "LDPC_Parity_Unstructured::cycle_removal_MGW(): " 00404 "Internal error"); 00405 set(j, l - ncheck, 1); 00406 set(p, k - ncheck, 1); 00407 00408 // Update adjacency matrix 00409 it_assert_debug(G(p, l) == 1 && G(l, p) == 1 && G(j, k) == 1 00410 && G(k, j) == 1, "G"); 00411 it_assert_debug(G(j, l) == 0 && G(l, j) == 0 && G(p, k) == 0 00412 && G(k, p) == 0, "G"); 00413 00414 // Delta is the update term to G 00415 Ssmat Delta(ncheck + nvar, ncheck + nvar, 2); 00416 Delta.set(j, k, -1); 00417 Delta.set(k, j, -1); 00418 Delta.set(p, l, -1); 00419 Delta.set(l, p, -1); 00420 Delta.set(j, l, 1); 00421 Delta.set(l, j, 1); 00422 Delta.set(p, k, 1); 00423 Delta.set(k, p, 1); 00424 00425 // update G and its powers 00426 G = G + Delta; 00427 it_assert_debug(G(p, l) == 0 && G(l, p) == 0 && G(j, k) == 0 00428 && G(k, j) == 0, "G"); 00429 it_assert_debug(G(j, l) == 1 && G(l, j) == 1 && G(p, k) == 1 00430 && G(k, p) == 1, "G"); 00431 00432 Gpow(1) = G; 00433 Gpow(2) = G * G; 00434 for (int z = 3; z <= r / 2; z++) { 00435 Gpow(z) = Gpow(z - 1) * G; 00436 } 00437 00438 traverse_again = true; 00439 } // if G()... 00440 } // loop over i 00441 } // loop over j 00442 if ((!traverse_again) && (cycles_found > 0)) { // no point continue 00443 scl = r - 2; 00444 goto finished; 00445 } 00446 } 00447 while (cycles_found != 0); 00448 scl = r; // there were no cycles of length r; move on to next r 00449 it_info_debug("Achieved girth " << (scl + 2) 00450 << ". Proceeding to next level."); 00451 } // loop over r 00452 00453 finished: 00454 int girth = scl + 2; // scl=length of smallest cycle 00455 it_info_debug("Cycle removal (MGW algoritm) finished. Graph girth: " 00456 << girth << ". Cycles remaining on next girth level: " 00457 << cycles_found); 00458 00459 return girth; 00460 } 00461 00462 void LDPC_Parity_Unstructured::generate_random_H(const ivec& C, 00463 const ivec& R, 00464 const ivec& cycopt) 00465 { 00466 // Method based on random permutation. Attempts to avoid placing new 00467 // edges so that cycles are created. More aggressive cycle avoidance 00468 // for degree-2 nodes. EGL January 2007. 00469 00470 initialize(sum(R), sum(C)); 00471 // C, R: Target number of columns/rows with certain number of ones 00472 00473 // compute number of edges 00474 int Ne = 0; 00475 for (int i = 0;i < C.length();i++) { 00476 for (int j = 0; j < C(i); j++) { 00477 for (int m = 0; m < i; m++) Ne++; 00478 } 00479 } 00480 00481 // compute connectivity matrix 00482 ivec vcon(Ne); 00483 ivec ccon(Ne); 00484 ivec vd(nvar); 00485 ivec cd(ncheck); 00486 int k = 0; 00487 int l = 0; 00488 for (int i = 0;i < C.length();i++) { 00489 for (int j = 0; j < C(i); j++) { 00490 for (int m = 0; m < i; m++) { 00491 vcon(k) = l; 00492 vd(l) = i; 00493 k++; 00494 } 00495 l++; 00496 } 00497 } 00498 k = 0; 00499 l = 0; 00500 for (int i = 0;i < R.length();i++) { 00501 for (int j = 0; j < R(i); j++) { 00502 for (int m = 0; m < i; m++) { 00503 ccon(k) = l; 00504 cd(l) = i; 00505 k++; 00506 } 00507 l++; 00508 } 00509 } 00510 it_assert(k == Ne, "C/R mismatch"); 00511 00512 // compute random permutations 00513 ivec ind = sort_index(randu(Ne)); 00514 ivec cp = sort_index(randu(nvar)); 00515 ivec rp = sort_index(randu(ncheck)); 00516 00517 // set the girth goal for various variable node degrees 00518 ivec Laim = zeros_i(Nmax); 00519 for (int i = 0; i < length(cycopt); i++) { 00520 Laim(i + 2) = cycopt(i); 00521 } 00522 for (int i = length(cycopt); i < Nmax - 2; i++) { 00523 Laim(i + 2) = cycopt(length(cycopt) - 1); 00524 } 00525 it_info_debug("Running with Laim=" << Laim.left(25)); 00526 00527 int failures = 0; 00528 const int Max_attempts = 100; 00529 const int apcl = 10; // attempts before reducing girth target 00530 for (int k = 0; k < Ne; k++) { 00531 const int el = Ne - k - 2; 00532 if (k % 250 == 0) { 00533 it_info_debug("Processing edge: " << k << " out of " << Ne 00534 << ". Variable node degree: " << vd(vcon(k)) 00535 << ". Girth target: " << Laim(vd(vcon(k))) 00536 << ". Accumulated failures: " << failures); 00537 } 00538 const int c = cp(vcon(k)); 00539 int L = Laim(vd(vcon(k))); 00540 int attempt = 0; 00541 while (true) { 00542 if (attempt > 0 && attempt % apcl == 0 && L >= 6) { L -= 2; }; 00543 int r = rp(ccon(ind(k))); 00544 if (get(r, c)) { // double edge 00545 // set(r,c,0); 00546 if (el > 0) { 00547 // int t=k+1+rand()%el; 00548 int t = k + 1 + randi(0, el - 1); 00549 int x = ind(t); 00550 ind(t) = ind(k); 00551 ind(k) = x; 00552 attempt++; 00553 if (attempt == Max_attempts) { 00554 failures++; 00555 break; 00556 } 00557 } 00558 else { // almost at the last edge 00559 break; 00560 } 00561 } 00562 else { 00563 set(r, c, 1); 00564 if (L > 0) { // attempt to avoid cycles 00565 if (check_connectivity(r, c, r, c, 0, L) >= 0) { // detected cycle 00566 set(r, c, 0); 00567 if (el > 0) { 00568 // make a swap in the index permutation 00569 // int t=k+1+rand()%el; 00570 int t = k + 1 + randi(0, el - 1); 00571 int x = ind(t); 00572 ind(t) = ind(k); 00573 ind(k) = x; 00574 attempt++; 00575 if (attempt == Max_attempts) { // give up 00576 failures++; 00577 set(r, c, 1); 00578 break; 00579 } 00580 } 00581 else { // no edges left 00582 set(r, c, 1); 00583 break; 00584 } 00585 } 00586 else { 00587 break; 00588 } 00589 } 00590 else { 00591 break; 00592 } 00593 } 00594 } 00595 } 00596 } 00597 00598 void LDPC_Parity_Unstructured::compute_CR(const vec& var_deg, const vec& chk_deg, const int Nvar, 00599 ivec &C, ivec &R) 00600 { 00601 // compute the degree distributions from a node perspective 00602 vec Vi = linspace(1, length(var_deg), length(var_deg)); 00603 vec Ci = linspace(1, length(chk_deg), length(chk_deg)); 00604 // Compute number of cols with n 1's 00605 // C, R: Target number of columns/rows with certain number of ones 00606 C = to_ivec(round(Nvar * elem_div(var_deg, Vi) 00607 / sum(elem_div(var_deg, Vi)))); 00608 C = concat(0, C); 00609 int edges = sum(elem_mult(to_ivec(linspace(0, C.length() - 1, 00610 C.length())), C)); 00611 R = to_ivec(round(edges * elem_div(chk_deg, Ci))); 00612 R = concat(0, R); 00613 vec Ri = linspace(0, length(R) - 1, length(R)); 00614 vec Coli = linspace(0, length(C) - 1, length(C)); 00615 00616 // trim to deal with inconsistencies due to rounding errors 00617 if (sum(C) != Nvar) { 00618 ivec ind = find(C == max(C)); 00619 C(ind(0)) = C(ind(0)) - (sum(C) - Nvar); 00620 } 00621 00622 //the number of edges calculated from R must match the number of 00623 //edges calculated from C 00624 while (sum(elem_mult(to_vec(R), Ri)) != 00625 sum(elem_mult(to_vec(C), Coli))) { 00626 //we're only changing R, this is probably(?) better for irac codes 00627 if (sum(elem_mult(to_vec(R), Ri)) > sum(elem_mult(to_vec(C), Coli))) { 00628 //remove an edge from R 00629 ivec ind = find(R == max(R)); 00630 int old = R(ind(0)); 00631 R.set(ind(0), old - 1); 00632 old = R(ind(0) - 1); 00633 R.set(ind(0) - 1, old + 1); 00634 } 00635 else { 00636 ivec ind = find(R == max(R)); 00637 if (ind(0) == R.length() - 1) { 00638 R = concat(R, 0); 00639 Ri = linspace(0, length(R) - 1, length(R)); 00640 } 00641 int old = R(ind(0)); 00642 R.set(ind(0), old - 1); 00643 old = R(ind(0) + 1); 00644 R.set(ind(0) + 1, old + 1); 00645 } 00646 } 00647 00648 C = concat(C, zeros_i(Nmax - length(C))); 00649 R = concat(R, zeros_i(Nmax - length(R))); 00650 00651 it_info_debug("C=" << C << std::endl); 00652 it_info_debug("R=" << R << std::endl); 00653 00654 } 00655 00656 // ---------------------------------------------------------------------- 00657 // LDPC_Parity_Regular 00658 // ---------------------------------------------------------------------- 00659 00660 LDPC_Parity_Regular::LDPC_Parity_Regular(int Nvar, int k, int l, 00661 const std::string& method, 00662 const ivec& options) 00663 { 00664 generate(Nvar, k, l, method, options); 00665 } 00666 00667 void LDPC_Parity_Regular::generate(int Nvar, int k, int l, 00668 const std::string& method, 00669 const ivec& options) 00670 { 00671 vec var_deg = zeros(k); 00672 vec chk_deg = zeros(l); 00673 var_deg(k - 1) = 1; 00674 chk_deg(l - 1) = 1; 00675 00676 ivec C, R; 00677 compute_CR(var_deg, chk_deg, Nvar, C, R); 00678 it_info_debug("sum(C)=" << sum(C) << " Nvar=" << Nvar); 00679 it_info_debug("sum(R)=" << sum(R) << " approximate target=" << round_i(Nvar * k / static_cast<double>(l))); 00680 00681 if (method == "rand") { 00682 generate_random_H(C, R, options); 00683 } 00684 else { 00685 it_error("not implemented"); 00686 }; 00687 } 00688 00689 00690 // ---------------------------------------------------------------------- 00691 // LDPC_Parity_Irregular 00692 // ---------------------------------------------------------------------- 00693 00694 LDPC_Parity_Irregular::LDPC_Parity_Irregular(int Nvar, 00695 const vec& var_deg, 00696 const vec& chk_deg, 00697 const std::string& method, 00698 const ivec& options) 00699 { 00700 generate(Nvar, var_deg, chk_deg, method, options); 00701 } 00702 00703 void LDPC_Parity_Irregular::generate(int Nvar, const vec& var_deg, 00704 const vec& chk_deg, 00705 const std::string& method, 00706 const ivec& options) 00707 { 00708 ivec C, R; 00709 compute_CR(var_deg, chk_deg, Nvar, C, R); 00710 00711 if (method == "rand") { 00712 generate_random_H(C, R, options); 00713 } 00714 else { 00715 it_error("not implemented"); 00716 }; 00717 } 00718 00719 // ---------------------------------------------------------------------- 00720 // BLDPC_Parity 00721 // ---------------------------------------------------------------------- 00722 00723 BLDPC_Parity::BLDPC_Parity(const imat& base_matrix, int exp_factor) 00724 { 00725 expand_base(base_matrix, exp_factor); 00726 } 00727 00728 BLDPC_Parity::BLDPC_Parity(const std::string& filename, int exp_factor) 00729 { 00730 load_base_matrix(filename); 00731 expand_base(H_b, exp_factor); 00732 } 00733 00734 void BLDPC_Parity::expand_base(const imat& base_matrix, int exp_factor) 00735 { 00736 Z = exp_factor; 00737 H_b = base_matrix; 00738 H_b_valid = true; 00739 00740 initialize(H_b.rows() * Z, H_b.cols() * Z); 00741 for (int r = 0; r < H_b.rows(); r++) { 00742 for (int c = 0; c < H_b.cols(); c++) { 00743 int rz = r * Z; 00744 int cz = c * Z; 00745 switch (H_b(r, c)) { 00746 case -1: 00747 break; 00748 case 0: 00749 for (int i = 0; i < Z; ++i) 00750 set(rz + i, cz + i, 1); 00751 break; 00752 default: 00753 for (int i = 0; i < Z; ++i) 00754 set(rz + i, cz + (i + H_b(r, c)) % Z, 1); 00755 break; 00756 } 00757 } 00758 } 00759 } 00760 00761 00762 int BLDPC_Parity::get_exp_factor() const 00763 { 00764 return Z; 00765 } 00766 00767 imat BLDPC_Parity::get_base_matrix() const 00768 { 00769 return H_b; 00770 } 00771 00772 void BLDPC_Parity::set_exp_factor(int exp_factor) 00773 { 00774 Z = exp_factor; 00775 if (H_b_valid) { 00776 expand_base(H_b, exp_factor); 00777 } 00778 else { 00779 calculate_base_matrix(); 00780 } 00781 } 00782 00783 00784 void BLDPC_Parity::load_base_matrix(const std::string& filename) 00785 { 00786 std::ifstream bm_file(filename.c_str()); 00787 it_assert(bm_file.is_open(), "BLDPC_Parity::load_base_matrix(): Could not " 00788 "open file \"" << filename << "\" for reading"); 00789 // delete old base matrix content 00790 H_b.set_size(0, 0); 00791 // parse text file content, row by row 00792 std::string line; 00793 int line_counter = 0; 00794 getline(bm_file, line); 00795 while (!bm_file.eof()) { 00796 line_counter++; 00797 std::stringstream ss(line); 00798 ivec row(0); 00799 while (ss.good()) { 00800 int val; 00801 ss >> val; 00802 row = concat(row, val); 00803 } 00804 if ((H_b.rows() == 0) || (row.size() == H_b.cols())) 00805 H_b.append_row(row); 00806 else 00807 it_warning("BLDPC_Parity::load_base_matrix(): Wrong size of " 00808 "a parsed row number " << line_counter); 00809 getline(bm_file, line); 00810 } 00811 bm_file.close(); 00812 00813 // transpose parsed base matrix if necessary 00814 if (H_b.rows() > H_b.cols()) 00815 H_b = H_b.transpose(); 00816 00817 H_b_valid = true; 00818 init_flag = false; 00819 } 00820 00821 void BLDPC_Parity::save_base_matrix(const std::string& filename) const 00822 { 00823 it_assert(H_b_valid, "BLDPC_Parity::save_base_matrix(): Base matrix is " 00824 "not valid"); 00825 std::ofstream bm_file(filename.c_str()); 00826 it_assert(bm_file.is_open(), "BLDPC_Parity::save_base_matrix(): Could not " 00827 "open file \"" << filename << "\" for writing"); 00828 00829 for (int r = 0; r < H_b.rows(); r++) { 00830 for (int c = 0; c < H_b.cols(); c++) { 00831 bm_file << std::setw(3) << H_b(r, c); 00832 } 00833 bm_file << "\n"; 00834 } 00835 00836 bm_file.close(); 00837 } 00838 00839 00840 void BLDPC_Parity::calculate_base_matrix() 00841 { 00842 std::string error_str = "BLDPC_Parity::calculate_base_matrix(): " 00843 "Invalid BLDPC matrix. Cannot calculate base matrix from it."; 00844 int rows = H.rows() / Z; 00845 int cols = H.cols() / Z; 00846 it_assert((rows * Z == H.rows()) && (cols * Z == H.cols()), error_str); 00847 H_b.set_size(rows, cols); 00848 00849 for (int r = 0; r < rows; ++r) { 00850 int rz = r * Z; 00851 for (int c = 0; c < cols; ++c) { 00852 int cz = c * Z; 00853 GF2mat_sparse H_Z = H.get_submatrix(rz, rz + Z - 1, cz, cz + Z - 1); 00854 00855 if (H_Z.nnz() == 0) { 00856 H_b(r, c) = -1; 00857 } 00858 else if (H_Z.nnz() == Z) { 00859 // check for cyclic-shifted ZxZ matrix 00860 int shift = 0; 00861 while ((shift < Z) && (H_Z(0, shift) != 1)) 00862 ++shift; 00863 it_assert(shift < Z, error_str); 00864 for (int i = 1; i < Z; ++i) 00865 it_assert(H_Z(0, shift) == H_Z(i, (i + shift) % Z), error_str); 00866 H_b(r, c) = shift; 00867 } 00868 else { 00869 it_error(error_str); 00870 } 00871 } // for (int c = 0; c < cols; ++c) 00872 } // for (int r = 0; r < rows; ++r) 00873 00874 it_info("Base matrix calculated"); 00875 H_b_valid = true; 00876 } 00877 00878 00879 00880 // ---------------------------------------------------------------------- 00881 // LDPC_Generator_Systematic 00882 // ---------------------------------------------------------------------- 00883 00884 LDPC_Generator_Systematic::LDPC_Generator_Systematic(LDPC_Parity* const H, 00885 bool natural_ordering, 00886 const ivec& ind): 00887 LDPC_Generator("systematic"), G() 00888 { 00889 ivec tmp; 00890 tmp = construct(H, natural_ordering, ind); 00891 } 00892 00893 00894 ivec LDPC_Generator_Systematic::construct(LDPC_Parity* const H, 00895 bool natural_ordering, 00896 const ivec& avoid_cols) 00897 { 00898 int nvar = H->get_nvar(); 00899 int ncheck = H->get_ncheck(); 00900 00901 // create dense representation of parity check matrix 00902 GF2mat Hd(H->get_H()); 00903 00904 // -- Determine initial column ordering -- 00905 ivec col_order(nvar); 00906 if (natural_ordering) { 00907 for (int i = 0; i < nvar; i++) { 00908 col_order(i) = i; 00909 } 00910 } 00911 else { 00912 // take the columns in random order, but the ones to avoid at last 00913 vec col_importance = randu(nvar); 00914 for (int i = 0; i < length(avoid_cols); i++) { 00915 col_importance(avoid_cols(i)) = (-col_importance(avoid_cols(i))); 00916 } 00917 col_order = sort_index(-col_importance); 00918 } 00919 00920 ivec actual_ordering(nvar); 00921 00922 // Now partition P as P=[P1 P2]. Then find G so [P1 P2][I G]'=0. 00923 00924 // -- Create P1 and P2 -- 00925 GF2mat P1; //(ncheck,nvar-ncheck); // non-invertible part 00926 GF2mat P2; //(ncheck,ncheck); // invertible part 00927 00928 it_info_debug("Computing a systematic generator matrix..."); 00929 00930 int j1 = 0, j2 = 0; 00931 int rank; 00932 ivec perm; 00933 GF2mat T, U; 00934 for (int k = 0; k < nvar; k++) { 00935 it_error_if(j1 >= nvar - ncheck, "LDPC_Generator_Systematic::construct(): " 00936 "Unable to obtain enough independent columns."); 00937 00938 bvec c = Hd.get_col(col_order(k)); 00939 if (j2 == 0) { // first column in P2 is number col_order(k) 00940 P2 = GF2mat(c); 00941 rank = P2.T_fact(T, U, perm); 00942 actual_ordering(k) = nvar - ncheck; 00943 j2++; 00944 } 00945 else { 00946 if (j2 < ncheck) { 00947 if (P2.T_fact_update_addcol(T, U, perm, c)) { 00948 P2 = P2.concatenate_horizontal(c); 00949 actual_ordering(k) = nvar - ncheck + j2; 00950 j2++; 00951 continue; 00952 } 00953 } 00954 if (j1 == 0) { 00955 P1 = GF2mat(c); 00956 actual_ordering(k) = j1; 00957 } 00958 else { 00959 P1 = P1.concatenate_horizontal(c); 00960 actual_ordering(k) = j1; 00961 } 00962 j1++; 00963 } 00964 } 00965 00966 it_info_debug("Rank of parity check matrix: " << j2); 00967 00968 // -- Compute the systematic part of the generator matrix -- 00969 G = (P2.inverse() * P1).transpose(); 00970 00971 // -- Permute the columns of the parity check matrix -- 00972 GF2mat P = P1.concatenate_horizontal(P2); 00973 *H = LDPC_Parity(ncheck, nvar); 00974 // brute force copy from P to H 00975 for (int i = 0; i < ncheck; i++) { 00976 for (int j = 0; j < nvar; j++) { 00977 if (P.get(i, j)) { 00978 H->set(i, j, 1); 00979 } 00980 } 00981 } 00982 00983 // -- Check that the result was correct -- 00984 it_assert_debug((GF2mat(H->get_H()) 00985 * (gf2dense_eye(nvar - ncheck).concatenate_horizontal(G)).transpose()).is_zero(), 00986 "LDPC_Generator_Systematic::construct(): Incorrect generator matrix G"); 00987 00988 G = G.transpose(); // store the generator matrix in transposed form 00989 it_info_debug("Systematic generator matrix computed."); 00990 00991 init_flag = true; 00992 00993 return actual_ordering; 00994 } 00995 00996 00997 void LDPC_Generator_Systematic::save(const std::string& filename) const 00998 { 00999 it_file f(filename); 01000 int ver; 01001 f >> Name("Fileversion") >> ver; 01002 it_assert(ver == LDPC_binary_file_version, 01003 "LDPC_Generator_Systematic::save(): Unsupported file format"); 01004 f << Name("G_type") << type; 01005 f << Name("G") << G; 01006 f.close(); 01007 } 01008 01009 01010 void LDPC_Generator_Systematic::load(const std::string& filename) 01011 { 01012 it_ifile f(filename); 01013 int ver; 01014 f >> Name("Fileversion") >> ver; 01015 it_assert(ver == LDPC_binary_file_version, 01016 "LDPC_Generator_Systematic::load(): Unsupported file format"); 01017 std::string gen_type; 01018 f >> Name("G_type") >> gen_type; 01019 it_assert(gen_type == type, 01020 "LDPC_Generator_Systematic::load(): Wrong generator type"); 01021 f >> Name("G") >> G; 01022 f.close(); 01023 01024 init_flag = true; 01025 } 01026 01027 01028 void LDPC_Generator_Systematic::encode(const bvec &input, bvec &output) 01029 { 01030 it_assert(init_flag, "LDPC_Generator_Systematic::encode(): Systematic " 01031 "generator not set up"); 01032 it_assert(input.size() == G.cols(), "LDPC_Generator_Systematic::encode(): " 01033 "Improper input vector size (" << input.size() << " != " 01034 << G.cols() << ")"); 01035 01036 output = concat(input, G * input); 01037 } 01038 01039 01040 // ---------------------------------------------------------------------- 01041 // BLDPC_Generator 01042 // ---------------------------------------------------------------------- 01043 01044 BLDPC_Generator::BLDPC_Generator(const BLDPC_Parity* const H, 01045 const std::string type): 01046 LDPC_Generator(type), H_enc(), N(0), M(0), K(0), Z(0) 01047 { 01048 construct(H); 01049 } 01050 01051 01052 void BLDPC_Generator::encode(const bvec &input, bvec &output) 01053 { 01054 it_assert(init_flag, "BLDPC_Generator::encode(): Cannot encode with not " 01055 "initialized generator"); 01056 it_assert(input.size() == K, "BLDPC_Generator::encode(): Input vector " 01057 "length is not equal to K"); 01058 01059 // copy systematic bits first 01060 output = input; 01061 output.set_size(N, true); 01062 01063 // backward substitution to obtain the first Z parity check bits 01064 for (int k = 0; k < Z; k++) { 01065 for (int j = 0; j < K; j++) { 01066 output(K + k) += H_enc(M - 1 - k, j) * input(j); 01067 } 01068 for (int j = 0; j < k; j++) { 01069 output(K + k) += H_enc(M - 1 - k, K + j) * output(K + j); 01070 } 01071 } 01072 01073 // forward substitution to obtain the next M-Z parity check bits 01074 for (int k = 0; k < M - Z; k++) { 01075 for (int j = 0; j < K; j++) { 01076 output(K + Z + k) += H_enc(k, j) * input(j); 01077 } 01078 for (int j = K; j < K + Z; j++) { 01079 output(K + Z + k) += H_enc(k, j) * output(j); 01080 } 01081 for (int j = K + Z; j < K + Z + k; j++) { 01082 output(K + Z + k) += H_enc(k, j) * output(j); 01083 } 01084 } 01085 } 01086 01087 01088 void BLDPC_Generator::construct(const BLDPC_Parity* const H) 01089 { 01090 if (H != 0 && H->is_valid()) { 01091 H_enc = H->get_H(); // sparse to dense conversion 01092 Z = H->get_exp_factor(); 01093 N = H_enc.cols(); 01094 M = H_enc.rows(); 01095 K = N - M; 01096 01097 // ---------------------------------------------------------------------- 01098 // STEP 1 01099 // ---------------------------------------------------------------------- 01100 // loop over last M-Z columns of matrix H 01101 for (int i = 0; i < M - Z; i += Z) { 01102 // loop over last Z rows of matrix H 01103 for (int j = 0; j < Z; j++) { 01104 // Gaussian elimination by adding particular rows 01105 H_enc.add_rows(M - 1 - j, M - Z - 1 - j - i); 01106 } 01107 } 01108 01109 // ---------------------------------------------------------------------- 01110 // STEP 2 01111 // ---------------------------------------------------------------------- 01112 // set first processed row index to M-Z 01113 int r1 = M - Z; 01114 // loop over columns with indexes K .. K+Z-1 01115 for (int c1 = K + Z - 1; c1 >= K; c1--) { 01116 int r2 = r1; 01117 // find the first '1' in column c1 01118 while (H_enc(r2, c1) == 0 && r2 < M - 1) 01119 r2++; 01120 // if necessary, swap rows r1 and r2 01121 if (r2 != r1) 01122 H_enc.swap_rows(r1, r2); 01123 01124 // look for the other '1' in column c1 and get rid of them 01125 for (r2 = r1 + 1; r2 < M; r2++) { 01126 if (H_enc(r2, c1) == 1) { 01127 // Gaussian elimination by adding particular rows 01128 H_enc.add_rows(r2, r1); 01129 } 01130 } 01131 01132 // increase first processed row index 01133 r1++; 01134 } 01135 01136 init_flag = true; 01137 } 01138 } 01139 01140 01141 void BLDPC_Generator::save(const std::string& filename) const 01142 { 01143 it_assert(init_flag, 01144 "BLDPC_Generator::save(): Can not save not initialized generator"); 01145 // Every Z-th row of H_enc until M-Z 01146 GF2mat H_T(M / Z - 1, N); 01147 for (int i = 0; i < M / Z - 1; i++) { 01148 H_T.set_row(i, H_enc.get_row(i*Z)); 01149 } 01150 // Last Z preprocessed rows of H_enc 01151 GF2mat H_Z = H_enc.get_submatrix(M - Z, 0, M - 1, N - 1); 01152 01153 it_file f(filename); 01154 int ver; 01155 f >> Name("Fileversion") >> ver; 01156 it_assert(ver == LDPC_binary_file_version, "BLDPC_Generator::save(): " 01157 "Unsupported file format"); 01158 f << Name("G_type") << type; 01159 f << Name("H_T") << H_T; 01160 f << Name("H_Z") << H_Z; 01161 f << Name("Z") << Z; 01162 f.close(); 01163 } 01164 01165 void BLDPC_Generator::load(const std::string& filename) 01166 { 01167 GF2mat H_T, H_Z; 01168 01169 it_ifile f(filename); 01170 int ver; 01171 f >> Name("Fileversion") >> ver; 01172 it_assert(ver == LDPC_binary_file_version, "BLDPC_Generator::load(): " 01173 "Unsupported file format"); 01174 std::string gen_type; 01175 f >> Name("G_type") >> gen_type; 01176 it_assert(gen_type == type, 01177 "BLDPC_Generator::load(): Wrong generator type"); 01178 f >> Name("H_T") >> H_T; 01179 f >> Name("H_Z") >> H_Z; 01180 f >> Name("Z") >> Z; 01181 f.close(); 01182 01183 N = H_T.cols(); 01184 M = (H_T.rows() + 1) * Z; 01185 K = N - M; 01186 01187 H_enc = GF2mat(M - Z, N); 01188 for (int i = 0; i < H_T.rows(); i++) { 01189 for (int j = 0; j < Z; j++) { 01190 for (int k = 0; k < N; k++) { 01191 if (H_T(i, (k / Z)*Z + (k + Z - j) % Z)) { 01192 H_enc.set(i*Z + j, k, 1); 01193 } 01194 } 01195 } 01196 } 01197 H_enc = H_enc.concatenate_vertical(H_Z); 01198 01199 init_flag = true; 01200 } 01201 01202 01203 // ---------------------------------------------------------------------- 01204 // LDPC_Code 01205 // ---------------------------------------------------------------------- 01206 01207 LDPC_Code::LDPC_Code(): H_defined(false), G_defined(false), dec_method("BP"), 01208 max_iters(50), psc(true), pisc(false), 01209 llrcalc(LLR_calc_unit()) { } 01210 01211 LDPC_Code::LDPC_Code(const LDPC_Parity* const H, 01212 LDPC_Generator* const G_in): 01213 H_defined(false), G_defined(false), dec_method("BP"), max_iters(50), 01214 psc(true), pisc(false), llrcalc(LLR_calc_unit()) 01215 { 01216 set_code(H, G_in); 01217 } 01218 01219 LDPC_Code::LDPC_Code(const std::string& filename, 01220 LDPC_Generator* const G_in): 01221 H_defined(false), G_defined(false), dec_method("BP"), max_iters(50), 01222 psc(true), pisc(false), llrcalc(LLR_calc_unit()) 01223 { 01224 load_code(filename, G_in); 01225 } 01226 01227 01228 void LDPC_Code::set_code(const LDPC_Parity* const H, 01229 LDPC_Generator* const G_in) 01230 { 01231 decoder_parameterization(H); 01232 setup_decoder(); 01233 G = G_in; 01234 if (G != 0) { 01235 G_defined = true; 01236 integrity_check(); 01237 } 01238 } 01239 01240 void LDPC_Code::load_code(const std::string& filename, 01241 LDPC_Generator* const G_in) 01242 { 01243 it_info_debug("LDPC_Code::load_code(): Loading LDPC codec from " 01244 << filename); 01245 01246 it_ifile f(filename); 01247 int ver; 01248 f >> Name("Fileversion") >> ver; 01249 it_assert(ver == LDPC_binary_file_version, "LDPC_Code::load_code(): " 01250 "Unsupported file format"); 01251 f >> Name("H_defined") >> H_defined; 01252 f >> Name("G_defined") >> G_defined; 01253 f >> Name("nvar") >> nvar; 01254 f >> Name("ncheck") >> ncheck; 01255 f >> Name("C") >> C; 01256 f >> Name("V") >> V; 01257 f >> Name("sumX1") >> sumX1; 01258 f >> Name("sumX2") >> sumX2; 01259 f >> Name("iind") >> iind; 01260 f >> Name("jind") >> jind; 01261 f.close(); 01262 01263 // load generator data 01264 if (G_defined) { 01265 it_assert(G_in != 0, "LDPC_Code::load_code(): Generator object is " 01266 "missing. Can not load the generator data from a file."); 01267 G = G_in; 01268 G->load(filename); 01269 } 01270 else { 01271 G = 0; 01272 it_info_debug("LDPC_Code::load_code(): Generator data not loaded. " 01273 "Generator object will not be used."); 01274 } 01275 01276 it_info_debug("LDPC_Code::load_code(): Successfully loaded LDPC codec " 01277 "from " << filename); 01278 01279 setup_decoder(); 01280 } 01281 01282 void LDPC_Code::save_code(const std::string& filename) const 01283 { 01284 it_assert(H_defined, "LDPC_Code::save_to_file(): There is no parity " 01285 "check matrix"); 01286 it_info_debug("LDPC_Code::save_to_file(): Saving LDPC codec to " 01287 << filename); 01288 01289 it_file f; 01290 f.open(filename, true); 01291 f << Name("Fileversion") << LDPC_binary_file_version; 01292 f << Name("H_defined") << H_defined; 01293 f << Name("G_defined") << G_defined; 01294 f << Name("nvar") << nvar; 01295 f << Name("ncheck") << ncheck; 01296 f << Name("C") << C; 01297 f << Name("V") << V; 01298 f << Name("sumX1") << sumX1; 01299 f << Name("sumX2") << sumX2; 01300 f << Name("iind") << iind; 01301 f << Name("jind") << jind; 01302 f.close(); 01303 01304 // save generator data; 01305 if (G_defined) 01306 G->save(filename); 01307 else 01308 it_info_debug("LDPC_Code::save_code(): Missing generator object - " 01309 "generator data not saved"); 01310 01311 it_info_debug("LDPC_Code::save_code(): Successfully saved LDPC codec to " 01312 << filename); 01313 } 01314 01315 01316 void LDPC_Code::set_decoding_method(const std::string& method_in) 01317 { 01318 it_assert((method_in == "bp") || (method_in == "BP"), 01319 "LDPC_Code::set_decoding_method(): Not implemented decoding method"); 01320 dec_method = method_in; 01321 } 01322 01323 void LDPC_Code::set_exit_conditions(int max_iters_in, 01324 bool syndr_check_each_iter, 01325 bool syndr_check_at_start) 01326 { 01327 it_assert(max_iters >= 0, "LDPC_Code::set_nrof_iterations(): Maximum " 01328 "number of iterations can not be negative"); 01329 max_iters = max_iters_in; 01330 psc = syndr_check_each_iter; 01331 pisc = syndr_check_at_start; 01332 } 01333 01334 void LDPC_Code::set_llrcalc(const LLR_calc_unit& llrcalc_in) 01335 { 01336 llrcalc = llrcalc_in; 01337 } 01338 01339 01340 void LDPC_Code::encode(const bvec &input, bvec &output) 01341 { 01342 it_assert(G_defined, "LDPC_Code::encode(): LDPC Generator is required " 01343 "for encoding"); 01344 G->encode(input, output); 01345 it_assert_debug(syndrome_check(output), "LDPC_Code::encode(): Syndrome " 01346 "check failed"); 01347 } 01348 01349 bvec LDPC_Code::encode(const bvec &input) 01350 { 01351 bvec result; 01352 encode(input, result); 01353 return result; 01354 } 01355 01356 void LDPC_Code::decode(const vec &llr_in, bvec &syst_bits) 01357 { 01358 QLLRvec qllrin = llrcalc.to_qllr(llr_in); 01359 QLLRvec qllrout; 01360 bp_decode(qllrin, qllrout); 01361 syst_bits = (qllrout.left(nvar - ncheck) < 0); 01362 } 01363 01364 bvec LDPC_Code::decode(const vec &llr_in) 01365 { 01366 bvec syst_bits; 01367 decode(llr_in, syst_bits); 01368 return syst_bits; 01369 } 01370 01371 void LDPC_Code::decode_soft_out(const vec &llr_in, vec &llr_out) 01372 { 01373 QLLRvec qllrin = llrcalc.to_qllr(llr_in); 01374 QLLRvec qllrout; 01375 bp_decode(qllrin, qllrout); 01376 llr_out = llrcalc.to_double(qllrout); 01377 } 01378 01379 vec LDPC_Code::decode_soft_out(const vec &llr_in) 01380 { 01381 vec llr_out; 01382 decode_soft_out(llr_in, llr_out); 01383 return llr_out; 01384 } 01385 01386 int LDPC_Code::bp_decode(const QLLRvec &LLRin, QLLRvec &LLRout) 01387 { 01388 // Note the IT++ convention that a sure zero corresponds to 01389 // LLR=+infinity 01390 it_assert(H_defined, "LDPC_Code::bp_decode(): Parity check matrix not " 01391 "defined"); 01392 it_assert((LLRin.size() == nvar) && (sumX1.size() == nvar) 01393 && (sumX2.size() == ncheck), "LDPC_Code::bp_decode(): Wrong " 01394 "input dimensions"); 01395 01396 if (pisc && syndrome_check(LLRin)) { 01397 LLRout = LLRin; 01398 return 0; 01399 } 01400 01401 LLRout.set_size(LLRin.size()); 01402 01403 // allocate temporary variables used for the check node update 01404 ivec jj(max_cnd); 01405 QLLRvec m(max_cnd); 01406 QLLRvec ml(max_cnd); 01407 QLLRvec mr(max_cnd); 01408 01409 // initial step 01410 for (int i = 0; i < nvar; i++) { 01411 int index = i; 01412 for (int j = 0; j < sumX1(i); j++) { 01413 mvc[index] = LLRin(i); 01414 index += nvar; 01415 } 01416 } 01417 01418 bool is_valid_codeword = false; 01419 int iter = 0; 01420 do { 01421 iter++; 01422 if (nvar >= 100000) { it_info_no_endl_debug("."); } 01423 // --------- Step 1: check to variable nodes ---------- 01424 for (int j = 0; j < ncheck; j++) { 01425 // The check node update calculations are hardcoded for degrees 01426 // up to 6. For larger degrees, a general algorithm is used. 01427 switch (sumX2(j)) { 01428 case 0: 01429 it_error("LDPC_Code::bp_decode(): sumX2(j)=0"); 01430 case 1: 01431 it_error("LDPC_Code::bp_decode(): sumX2(j)=1"); 01432 case 2: { 01433 mcv[j+ncheck] = mvc[jind[j]]; 01434 mcv[j] = mvc[jind[j+ncheck]]; 01435 break; 01436 } 01437 case 3: { 01438 int j0 = j; 01439 QLLR m0 = mvc[jind[j0]]; 01440 int j1 = j0 + ncheck; 01441 QLLR m1 = mvc[jind[j1]]; 01442 int j2 = j1 + ncheck; 01443 QLLR m2 = mvc[jind[j2]]; 01444 mcv[j0] = llrcalc.Boxplus(m1, m2); 01445 mcv[j1] = llrcalc.Boxplus(m0, m2); 01446 mcv[j2] = llrcalc.Boxplus(m0, m1); 01447 break; 01448 } 01449 case 4: { 01450 int j0 = j; 01451 QLLR m0 = mvc[jind[j0]]; 01452 int j1 = j0 + ncheck; 01453 QLLR m1 = mvc[jind[j1]]; 01454 int j2 = j1 + ncheck; 01455 QLLR m2 = mvc[jind[j2]]; 01456 int j3 = j2 + ncheck; 01457 QLLR m3 = mvc[jind[j3]]; 01458 QLLR m01 = llrcalc.Boxplus(m0, m1); 01459 QLLR m23 = llrcalc.Boxplus(m2, m3); 01460 mcv[j0] = llrcalc.Boxplus(m1, m23); 01461 mcv[j1] = llrcalc.Boxplus(m0, m23); 01462 mcv[j2] = llrcalc.Boxplus(m01, m3); 01463 mcv[j3] = llrcalc.Boxplus(m01, m2); 01464 break; 01465 } 01466 case 5: { 01467 int j0 = j; 01468 QLLR m0 = mvc[jind[j0]]; 01469 int j1 = j0 + ncheck; 01470 QLLR m1 = mvc[jind[j1]]; 01471 int j2 = j1 + ncheck; 01472 QLLR m2 = mvc[jind[j2]]; 01473 int j3 = j2 + ncheck; 01474 QLLR m3 = mvc[jind[j3]]; 01475 int j4 = j3 + ncheck; 01476 QLLR m4 = mvc[jind[j4]]; 01477 QLLR m01 = llrcalc.Boxplus(m0, m1); 01478 QLLR m02 = llrcalc.Boxplus(m01, m2); 01479 QLLR m34 = llrcalc.Boxplus(m3, m4); 01480 QLLR m24 = llrcalc.Boxplus(m2, m34); 01481 mcv[j0] = llrcalc.Boxplus(m1, m24); 01482 mcv[j1] = llrcalc.Boxplus(m0, m24); 01483 mcv[j2] = llrcalc.Boxplus(m01, m34); 01484 mcv[j3] = llrcalc.Boxplus(m02, m4); 01485 mcv[j4] = llrcalc.Boxplus(m02, m3); 01486 break; 01487 } 01488 case 6: { 01489 int j0 = j; 01490 QLLR m0 = mvc[jind[j0]]; 01491 int j1 = j0 + ncheck; 01492 QLLR m1 = mvc[jind[j1]]; 01493 int j2 = j1 + ncheck; 01494 QLLR m2 = mvc[jind[j2]]; 01495 int j3 = j2 + ncheck; 01496 QLLR m3 = mvc[jind[j3]]; 01497 int j4 = j3 + ncheck; 01498 QLLR m4 = mvc[jind[j4]]; 01499 int j5 = j4 + ncheck; 01500 QLLR m5 = mvc[jind[j5]]; 01501 QLLR m01 = llrcalc.Boxplus(m0, m1); 01502 QLLR m23 = llrcalc.Boxplus(m2, m3); 01503 QLLR m45 = llrcalc.Boxplus(m4, m5); 01504 QLLR m03 = llrcalc.Boxplus(m01, m23); 01505 QLLR m25 = llrcalc.Boxplus(m23, m45); 01506 QLLR m0145 = llrcalc.Boxplus(m01, m45); 01507 mcv[j0] = llrcalc.Boxplus(m1, m25); 01508 mcv[j1] = llrcalc.Boxplus(m0, m25); 01509 mcv[j2] = llrcalc.Boxplus(m0145, m3); 01510 mcv[j3] = llrcalc.Boxplus(m0145, m2); 01511 mcv[j4] = llrcalc.Boxplus(m03, m5); 01512 mcv[j5] = llrcalc.Boxplus(m03, m4); 01513 break; 01514 } 01515 default: { 01516 int nodes = sumX2(j); 01517 if( nodes > max_cnd ) { 01518 std::ostringstream m_sout; 01519 m_sout << "check node degrees >" << max_cnd << " not supported in this version"; 01520 it_error( m_sout.str() ); 01521 } 01522 01523 nodes--; 01524 jj[0] = j; 01525 m[0] = mvc[jind[jj[0]]]; 01526 for(int i = 1; i <= nodes; i++ ) { 01527 jj[i] = jj[i-1] + ncheck; 01528 m[i] = mvc[jind[jj[i]]]; 01529 } 01530 01531 // compute partial sums from the left and from the right 01532 ml[0] = m[0]; 01533 mr[0] = m[nodes]; 01534 for(int i = 1; i < nodes; i++ ) { 01535 ml[i] = llrcalc.Boxplus( ml[i-1], m[i] ); 01536 mr[i] = llrcalc.Boxplus( mr[i-1], m[nodes-i] ); 01537 } 01538 01539 // merge partial sums 01540 mcv[jj[0]] = mr[nodes-1]; 01541 mcv[jj[nodes]] = ml[nodes-1]; 01542 for(int i = 1; i < nodes; i++ ) 01543 mcv[jj[i]] = llrcalc.Boxplus( ml[i-1], mr[nodes-1-i] ); 01544 } 01545 } // switch statement 01546 } 01547 01548 // step 2: variable to check nodes 01549 for (int i = 0; i < nvar; i++) { 01550 switch (sumX1(i)) { 01551 case 0: 01552 it_error("LDPC_Code::bp_decode(): sumX1(i)=0"); 01553 case 1: { 01554 // Degenerate case-should not occur for good coded. A lonely 01555 // variable node only sends its incoming message 01556 mvc[i] = LLRin(i); 01557 LLRout(i) = LLRin(i); 01558 break; 01559 } 01560 case 2: { 01561 QLLR m0 = mcv[iind[i]]; 01562 int i1 = i + nvar; 01563 QLLR m1 = mcv[iind[i1]]; 01564 mvc[i] = LLRin(i) + m1; 01565 mvc[i1] = LLRin(i) + m0; 01566 LLRout(i) = mvc[i1] + m1; 01567 break; 01568 } 01569 case 3: { 01570 int i0 = i; 01571 QLLR m0 = mcv[iind[i0]]; 01572 int i1 = i0 + nvar; 01573 QLLR m1 = mcv[iind[i1]]; 01574 int i2 = i1 + nvar; 01575 QLLR m2 = mcv[iind[i2]]; 01576 LLRout(i) = LLRin(i) + m0 + m1 + m2; 01577 mvc[i0] = LLRout(i) - m0; 01578 mvc[i1] = LLRout(i) - m1; 01579 mvc[i2] = LLRout(i) - m2; 01580 break; 01581 } 01582 case 4: { 01583 int i0 = i; 01584 QLLR m0 = mcv[iind[i0]]; 01585 int i1 = i0 + nvar; 01586 QLLR m1 = mcv[iind[i1]]; 01587 int i2 = i1 + nvar; 01588 QLLR m2 = mcv[iind[i2]]; 01589 int i3 = i2 + nvar; 01590 QLLR m3 = mcv[iind[i3]]; 01591 LLRout(i) = LLRin(i) + m0 + m1 + m2 + m3; 01592 mvc[i0] = LLRout(i) - m0; 01593 mvc[i1] = LLRout(i) - m1; 01594 mvc[i2] = LLRout(i) - m2; 01595 mvc[i3] = LLRout(i) - m3; 01596 break; 01597 } 01598 default: { // differential update 01599 QLLR mvc_temp = LLRin(i); 01600 int index_iind = i; // tracks i+jp*nvar 01601 for (int jp = 0; jp < sumX1(i); jp++) { 01602 mvc_temp += mcv[iind[index_iind]]; 01603 index_iind += nvar; 01604 } 01605 LLRout(i) = mvc_temp; 01606 index_iind = i; // tracks i+j*nvar 01607 for (int j = 0; j < sumX1[i]; j++) { 01608 mvc[index_iind] = mvc_temp - mcv[iind[index_iind]]; 01609 index_iind += nvar; 01610 } 01611 } 01612 } 01613 } 01614 01615 if (psc && syndrome_check(LLRout)) { 01616 is_valid_codeword = true; 01617 break; 01618 } 01619 } 01620 while (iter < max_iters); 01621 01622 if (nvar >= 100000) { it_info_debug(""); } 01623 return (is_valid_codeword ? iter : -iter); 01624 } 01625 01626 01627 bool LDPC_Code::syndrome_check(const bvec &x) const 01628 { 01629 QLLRvec llr = 1 - 2 * to_ivec(x); 01630 return syndrome_check(llr); 01631 } 01632 01633 bool LDPC_Code::syndrome_check(const QLLRvec &LLR) const 01634 { 01635 // Please note the IT++ convention that a sure zero corresponds to 01636 // LLR=+infinity 01637 int i, j, synd, vi; 01638 01639 for (j = 0; j < ncheck; j++) { 01640 synd = 0; 01641 int vind = j; // tracks j+i*ncheck 01642 for (i = 0; i < sumX2(j); i++) { 01643 vi = V(vind); 01644 if (LLR(vi) < 0) { 01645 synd++; 01646 } 01647 vind += ncheck; 01648 } 01649 if ((synd&1) == 1) { 01650 return false; // codeword is invalid 01651 } 01652 } 01653 return true; // codeword is valid 01654 }; 01655 01656 01657 // ---------------------------------------------------------------------- 01658 // LDPC_Code private methods 01659 // ---------------------------------------------------------------------- 01660 01661 void LDPC_Code::decoder_parameterization(const LDPC_Parity* const Hmat) 01662 { 01663 // copy basic parameters 01664 sumX1 = Hmat->sumX1; 01665 sumX2 = Hmat->sumX2; 01666 nvar = Hmat->nvar; //get_nvar(); 01667 ncheck = Hmat->ncheck; //get_ncheck(); 01668 int cmax = max(sumX1); 01669 int vmax = max(sumX2); 01670 01671 // decoder parameterization 01672 V = zeros_i(ncheck * vmax); 01673 C = zeros_i(cmax * nvar); 01674 jind = zeros_i(ncheck * vmax); 01675 iind = zeros_i(nvar * cmax); 01676 01677 it_info_debug("LDPC_Code::decoder_parameterization(): Computations " 01678 "- phase 1"); 01679 for (int i = 0; i < nvar; i++) { 01680 ivec coli = Hmat->get_col(i).get_nz_indices(); 01681 for (int j0 = 0; j0 < length(coli); j0++) { 01682 C(j0 + cmax*i) = coli(j0); 01683 } 01684 } 01685 01686 it_info_debug("LDPC_Code::decoder_parameterization(): Computations " 01687 "- phase 2"); 01688 it_info_debug("Computing decoder parameterization. Phase 2"); 01689 for (int j = 0; j < ncheck; j++) { 01690 ivec rowj = Hmat->get_row(j).get_nz_indices(); 01691 for (int i0 = 0; i0 < length(rowj); i0++) { 01692 V(j + ncheck*i0) = rowj(i0); 01693 } 01694 } 01695 01696 it_info_debug("LDPC_Code::decoder_parameterization(): Computations " 01697 "- phase 3"); 01698 it_info_debug("Computing decoder parameterization. Phase 3."); 01699 for (int j = 0; j < ncheck; j++) { 01700 for (int ip = 0; ip < sumX2(j); ip++) { 01701 int vip = V(j + ip * ncheck); 01702 int k = 0; 01703 while (1 == 1) { 01704 if (C(k + vip*cmax) == j) { 01705 break; 01706 } 01707 k++; 01708 } 01709 jind(j + ip*ncheck) = vip + k * nvar; 01710 } 01711 } 01712 01713 it_info_debug("LDPC_Code::decoder_parameterization(): Computations " 01714 "- phase 4"); 01715 for (int i = 0; i < nvar; i++) { 01716 for (int jp = 0; jp < sumX1(i); jp++) { 01717 int cjp = C(jp + i * cmax); 01718 int k = 0; 01719 while (1 == 1) { 01720 if (V(cjp + k*ncheck) == i) {break; } 01721 k++; 01722 } 01723 iind(i + jp*nvar) = cjp + k * ncheck; 01724 } 01725 } 01726 01727 H_defined = true; 01728 } 01729 01730 01731 void LDPC_Code::setup_decoder() 01732 { 01733 if (H_defined) { 01734 mcv.set_size(max(sumX2) * ncheck); 01735 mvc.set_size(max(sumX1) * nvar); 01736 } 01737 } 01738 01739 01740 void LDPC_Code::integrity_check() 01741 { 01742 if (G_defined) { 01743 it_info_debug("LDPC_Code::integrity_check(): Checking integrity of " 01744 "the LDPC_Parity and LDPC_Generator data"); 01745 bvec bv(nvar - ncheck), cw; 01746 bv.clear(); 01747 bv(0) = 1; 01748 for (int i = 0; i < nvar - ncheck; i++) { 01749 G->encode(bv, cw); 01750 it_assert(syndrome_check(cw), 01751 "LDPC_Code::integrity_check(): Syndrome check failed"); 01752 bv.shift_right(bv(nvar - ncheck - 1)); 01753 } 01754 } 01755 else { 01756 it_info_debug("LDPC_Code::integrity_check(): No generator defined " 01757 "- no check performed"); 01758 } 01759 } 01760 01761 // ---------------------------------------------------------------------- 01762 // Related functions 01763 // ---------------------------------------------------------------------- 01764 01765 std::ostream &operator<<(std::ostream &os, const LDPC_Code &C) 01766 { 01767 ivec rdeg = zeros_i(max(C.sumX2) + 1); 01768 for (int i = 0; i < C.ncheck; i++) { 01769 rdeg(C.sumX2(i))++; 01770 } 01771 01772 ivec cdeg = zeros_i(max(C.sumX1) + 1); 01773 for (int j = 0; j < C.nvar; j++) { 01774 cdeg(C.sumX1(j))++; 01775 } 01776 01777 os << "--- LDPC codec ----------------------------------\n" 01778 << "Nvar : " << C.get_nvar() << "\n" 01779 << "Ncheck : " << C.get_ncheck() << "\n" 01780 << "Rate : " << C.get_rate() << "\n" 01781 << "Column degrees (node perspective): " << cdeg << "\n" 01782 << "Row degrees (node perspective): " << rdeg << "\n" 01783 << "-------------------------------------------------\n" 01784 << "Decoder parameters:\n" 01785 << " - method : " << C.dec_method << "\n" 01786 << " - max. iterations : " << C.max_iters << "\n" 01787 << " - syndrome check at each iteration : " << C.psc << "\n" 01788 << " - syndrome check at start : " << C.pisc << "\n" 01789 << "-------------------------------------------------\n" 01790 << C.llrcalc << "\n"; 01791 return os; 01792 } 01793 01794 } // namespace itpp
Generated on Sat Jul 9 2011 15:21:31 for IT++ by Doxygen 1.7.4