IT++ Logo
ldpc.cpp
Go to the documentation of this file.
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
 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