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

Generated on Sat May 3 16:10:42 2008 for IT++ by Doxygen 1.5.5