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

Generated on Tue Dec 6 2011 16:51:54 for IT++ by Doxygen 1.7.4