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

Generated on Sun Dec 9 17:38:46 2007 for IT++ by Doxygen 1.5.4