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     bmat H_dense(H_b.rows() * Z, H_b.cols() * Z);
00727 
00728     for (int r = 0; r < H_b.rows(); r++)
00729       for (int c = 0; c < H_b.cols(); c++)
00730         switch (H_b(r, c)) {
00731         case -1:
00732           H_dense.set_submatrix(r*Z, c*Z, zeros_b(Z, Z));
00733           break;
00734         case 0:
00735           H_dense.set_submatrix(r*Z, c*Z, eye_b(Z));
00736           break;
00737         default:
00738           H_dense.set_submatrix(r*Z, c*Z, circular_eye_b(Z, H_b(r, c)));
00739           break;
00740         }
00741 
00742     initialize(H_dense.rows(), H_dense.cols());
00743 
00744     for (int r = 0; r < ncheck; r++)
00745       for (int c = 0; c < nvar; c++)
00746         if (H_dense(r, c))
00747           set(r, c, 1);
00748   }
00749 
00750 
00751   int BLDPC_Parity::get_exp_factor() const
00752   {
00753     return Z;
00754   }
00755 
00756   imat BLDPC_Parity::get_base_matrix() const
00757   {
00758     return H_b;
00759   }
00760 
00761   void BLDPC_Parity::set_exp_factor(int exp_factor)
00762   {
00763     Z = exp_factor;
00764     if (H_b_valid) {
00765       expand_base(H_b, exp_factor);
00766     }
00767     else {
00768       calculate_base_matrix();
00769     }
00770   }
00771 
00772 
00773   void BLDPC_Parity::load_base_matrix(const std::string& filename)
00774   {
00775     std::ifstream bm_file(filename.c_str());
00776     it_assert(bm_file.is_open(), "BLDPC_Parity::load_base_matrix(): Could not "
00777               "open file \"" << filename << "\" for reading");
00778     // delete old base matrix content
00779     H_b.set_size(0, 0);
00780     // parse text file content, row by row
00781     std::string line;
00782     int line_counter = 0;
00783     getline(bm_file, line);
00784     while (!bm_file.eof()) {
00785       line_counter++;
00786       std::stringstream ss(line);
00787       ivec row(0);
00788       while (ss.good()) {
00789         int val;
00790         ss >> val;
00791         row = concat(row, val);
00792       }
00793       if ((H_b.rows() == 0) || (row.size() == H_b.cols()))
00794         H_b.append_row(row);
00795       else
00796         it_warning("BLDPC_Parity::load_base_matrix(): Wrong size of "
00797                    "a parsed row number " << line_counter);
00798       getline(bm_file, line);
00799     }
00800     bm_file.close();
00801 
00802     // transpose parsed base matrix if necessary
00803     if (H_b.rows() > H_b.cols())
00804       H_b = H_b.transpose();
00805 
00806     H_b_valid = true;
00807     init_flag = false;
00808   }
00809 
00810   void BLDPC_Parity::save_base_matrix(const std::string& filename) const
00811   {
00812     it_assert(H_b_valid, "BLDPC_Parity::save_base_matrix(): Base matrix is "
00813               "not valid");
00814     std::ofstream bm_file(filename.c_str());
00815     it_assert(bm_file.is_open(), "BLDPC_Parity::save_base_matrix(): Could not "
00816               "open file \"" << filename << "\" for writing");
00817 
00818     for (int r = 0; r < H_b.rows(); r++) {
00819       for (int c = 0; c < H_b.cols(); c++) {
00820         bm_file << std::setw(3) << H_b(r, c);
00821       }
00822       bm_file << "\n";
00823     }
00824 
00825     bm_file.close();
00826   }
00827 
00828 
00829   bmat BLDPC_Parity::circular_eye_b(int size, int shift)
00830   {
00831     bmat I = eye_b(size);
00832     int c = I.cols();
00833     int s = shift % c;
00834     if (s > 0)
00835       return concat_horizontal(I.get_cols(c-s, c-1), I.get_cols(0, c-1-s));
00836     else
00837       return I;
00838   }
00839 
00840   void BLDPC_Parity::calculate_base_matrix()
00841   {
00842     bmat H_dense = H.full();
00843     int rows = H_dense.rows() / Z;
00844     int cols = H_dense.cols() / Z;
00845     it_assert((rows * Z == H_dense.rows()) && (cols * Z == H_dense.cols()),
00846               "BLDPC_Parity::calculate_base_matrix(): Parity check matrix "
00847               "does not match an expansion factor");
00848     H_b.set_size(rows, cols);
00849 
00850     for (int r = 0; r < rows; r++) {
00851       for (int c = 0; c < cols; c++) {
00852         bmat tmp = H_dense.get(r*Z, (r+1)*Z-1, c*Z, (c+1)*Z-1);
00853         if (tmp == zeros_b(Z, Z)) {
00854           H_b(r, c) = -1;
00855         }
00856         else {
00857           int shift = 0;
00858           while (tmp != circular_eye_b(Z, shift) && shift < Z) {
00859             shift++;
00860           }
00861           it_assert(shift < Z, "BLDPC_Parity::calculate_base_matrix(): "
00862                     "Parity check matrix was not constructed with expansion "
00863                     "factor Z = " << Z);
00864           H_b(r, c) = shift;
00865         }
00866       }
00867     }
00868 
00869     H_b_valid = true;
00870   }
00871 
00872 
00873 
00874   // ----------------------------------------------------------------------
00875   // LDPC_Generator_Systematic
00876   // ----------------------------------------------------------------------
00877 
00878   LDPC_Generator_Systematic::LDPC_Generator_Systematic(LDPC_Parity* const H,
00879                                                        bool natural_ordering,
00880                                                        const ivec& ind):
00881     LDPC_Generator("systematic"), G()
00882   {
00883     ivec tmp;
00884     tmp = construct(H, natural_ordering, ind);
00885   }
00886 
00887 
00888   ivec LDPC_Generator_Systematic::construct(LDPC_Parity* const H,
00889                                             bool natural_ordering,
00890                                             const ivec& avoid_cols)
00891   {
00892     int nvar = H->get_nvar();
00893     int ncheck = H->get_ncheck();
00894 
00895     // create dense representation of parity check matrix
00896     GF2mat Hd(H->get_H());
00897 
00898     // -- Determine initial column ordering --
00899     ivec col_order(nvar);
00900     if (natural_ordering) {
00901       for (int i=0; i<nvar; i++) {
00902         col_order(i)=i;
00903       }
00904     }
00905     else {
00906       // take the columns in random order, but the ones to avoid at last
00907       vec col_importance = randu(nvar);
00908       for (int i=0; i<length(avoid_cols); i++) {
00909         col_importance(avoid_cols(i)) = (-col_importance(avoid_cols(i)));
00910       }
00911       col_order = sort_index(-col_importance);
00912     }
00913 
00914     ivec actual_ordering(nvar);
00915 
00916     // Now partition P as P=[P1 P2]. Then find G so [P1 P2][I G]'=0.
00917 
00918     // -- Create P1 and P2 --
00919     GF2mat P1; //(ncheck,nvar-ncheck);      // non-invertible part
00920     GF2mat P2; //(ncheck,ncheck);           // invertible part
00921 
00922     it_info_debug("Computing a systematic generator matrix...");
00923 
00924     int j1=0, j2=0;
00925     int rank;
00926     ivec perm;
00927     GF2mat T, U;
00928     for (int k=0; k<nvar; k++) {
00929       it_error_if(j1 >= nvar-ncheck, "LDPC_Generator_Systematic::construct(): "
00930                   "Unable to obtain enough independent columns.");
00931 
00932       bvec c = Hd.get_col(col_order(k));
00933       if (j2==0) {       // first column in P2 is number col_order(k)
00934         P2 = GF2mat(c);
00935         rank = P2.T_fact(T,U,perm);
00936         actual_ordering(k)=nvar-ncheck;
00937         j2++;
00938       }
00939       else {
00940         if (j2<ncheck) {
00941           if (P2.T_fact_update_addcol(T,U,perm,c)) {
00942             P2 = P2.concatenate_horizontal(c);
00943             actual_ordering(k)=nvar-ncheck+j2;
00944             j2++;
00945             continue;
00946           }
00947         }
00948         if (j1==0) {
00949           P1 = GF2mat(c);
00950           actual_ordering(k)=j1;
00951         }
00952         else {
00953           P1 = P1.concatenate_horizontal(c);
00954           actual_ordering(k)=j1;
00955         }
00956         j1++;
00957       }
00958     }
00959 
00960     it_info_debug("Rank of parity check matrix: " << j2);
00961 
00962     // -- Compute the systematic part of the generator matrix --
00963     G = (P2.inverse()*P1).transpose();
00964 
00965     // -- Permute the columns of the parity check matrix --
00966     GF2mat P = P1.concatenate_horizontal(P2);
00967     *H = LDPC_Parity(ncheck, nvar);
00968     // brute force copy from P to H
00969     for (int i=0; i<ncheck; i++) {
00970       for (int j=0; j<nvar; j++) {
00971         if (P.get(i,j)) {
00972           H->set(i,j,1);
00973         }
00974       }
00975     }
00976 
00977     // -- Check that the result was correct --
00978     it_assert_debug((GF2mat(H->get_H())
00979                      * (gf2dense_eye(nvar-ncheck).concatenate_horizontal(G)).transpose()).is_zero(),
00980                     "LDPC_Generator_Systematic::construct(): Incorrect generator matrix G");
00981 
00982     G = G.transpose();  // store the generator matrix in transposed form
00983     it_info_debug("Systematic generator matrix computed.");
00984 
00985     init_flag = true;
00986 
00987     return actual_ordering;
00988   }
00989 
00990 
00991   void LDPC_Generator_Systematic::save(const std::string& filename) const
00992   {
00993     it_file f(filename);
00994     int ver;
00995     f >> Name("Fileversion") >> ver;
00996     it_assert(ver == LDPC_binary_file_version,
00997               "LDPC_Generator_Systematic::save(): Unsupported file format");
00998     f << Name("G_type") << type;
00999     f << Name("G") << G;
01000     f.close();
01001   }
01002 
01003 
01004   void LDPC_Generator_Systematic::load(const std::string& filename)
01005   {
01006     it_ifile f(filename);
01007     int ver;
01008     f >> Name("Fileversion") >> ver;
01009     it_assert(ver == LDPC_binary_file_version,
01010               "LDPC_Generator_Systematic::load(): Unsupported file format");
01011     std::string gen_type;
01012     f >> Name("G_type") >> gen_type;
01013     it_assert(gen_type == type,
01014               "LDPC_Generator_Systematic::load(): Wrong generator type");
01015     f >> Name("G") >> G;
01016     f.close();
01017 
01018     init_flag = true;
01019   }
01020 
01021 
01022   void LDPC_Generator_Systematic::encode(const bvec &input, bvec &output)
01023   {
01024     it_assert(init_flag, "LDPC_Generator_Systematic::encode(): Systematic "
01025               "generator not set up");
01026     it_assert(input.size() == G.cols(), "LDPC_Generator_Systematic::encode(): "
01027               "Improper input vector size (" << input.size() << " != "
01028               << G.cols() << ")");
01029 
01030     output = concat(input, G * input);
01031   }
01032 
01033 
01034   // ----------------------------------------------------------------------
01035   // BLDPC_Generator
01036   // ----------------------------------------------------------------------
01037 
01038   BLDPC_Generator::BLDPC_Generator(const BLDPC_Parity* const H,
01039                                    const std::string type):
01040     LDPC_Generator(type), H_enc(), N(0), M(0), K(0), Z(0)
01041   {
01042     construct(H);
01043   }
01044 
01045 
01046   void BLDPC_Generator::encode(const bvec &input, bvec &output)
01047   {
01048     it_assert(init_flag, "BLDPC_Generator::encode(): Cannot encode with not "
01049               "initialized generator");
01050     it_assert(input.size() == K, "BLDPC_Generator::encode(): Input vector "
01051               "length is not equal to K");
01052 
01053     // copy systematic bits first
01054     output = input;
01055     output.set_size(N, true);
01056 
01057     // backward substitution to obtain the first Z parity check bits
01058     for (int k = 0; k < Z; k++) {
01059       for (int j = 0; j < K; j++) {
01060         output(K+k) += H_enc(M-1-k, j) * input(j);
01061       }
01062       for (int j = 0; j < k; j++) {
01063         output(K+k) += H_enc(M-1-k, K+j) * output(K+j);
01064       }
01065     }
01066 
01067     // forward substitution to obtain the next M-Z parity check bits
01068     for (int k = 0; k < M-Z; k++) {
01069       for (int j = 0; j < K; j++) {
01070         output(K+Z+k) += H_enc(k, j) * input(j);
01071       }
01072       for (int j = K; j < K+Z; j++) {
01073         output(K+Z+k) += H_enc(k, j) * output(j);
01074       }
01075       for (int j = K+Z; j < K+Z+k; j++) {
01076         output(K+Z+k) += H_enc(k, j) * output(j);
01077       }
01078     }
01079   }
01080 
01081 
01082   void BLDPC_Generator::construct(const BLDPC_Parity* const H)
01083   {
01084     if (H != 0 && H->is_valid()) {
01085       H_enc = H->get_H(); // sparse to dense conversion
01086       Z = H->get_exp_factor();
01087       N = H_enc.cols();
01088       M = H_enc.rows();
01089       K = N - M;
01090 
01091       // ----------------------------------------------------------------------
01092       // STEP 1
01093       // ----------------------------------------------------------------------
01094       // loop over last M-Z columns of matrix H
01095       for (int i = 0; i < M-Z; i += Z) {
01096         // loop over last Z rows of matrix H
01097         for (int j = 0; j < Z; j++) {
01098           // Gaussian elimination by adding particular rows
01099           H_enc.add_rows(M-1-j, M-Z-1-j-i);
01100         }
01101       }
01102 
01103       // ----------------------------------------------------------------------
01104       // STEP 2
01105       // ----------------------------------------------------------------------
01106       // set first processed row index to M-Z
01107       int r1 = M-Z;
01108       // loop over columns with indexes K .. K+Z-1
01109       for (int c1 = K+Z-1; c1 >= K; c1--) {
01110         int r2 = r1;
01111         // find the first '1' in column c1
01112         while (H_enc(r2, c1) == 0 && r2 < M-1)
01113           r2++;
01114         // if necessary, swap rows r1 and r2
01115         if (r2 != r1)
01116           H_enc.swap_rows(r1, r2);
01117 
01118         // look for the other '1' in column c1 and get rid of them
01119         for (r2 = r1+1; r2 < M; r2++) {
01120           if (H_enc(r2, c1) == 1) {
01121             // Gaussian elimination by adding particular rows
01122             H_enc.add_rows(r2, r1);
01123           }
01124         }
01125 
01126         // increase first processed row index
01127         r1++;
01128       }
01129 
01130       init_flag = true;
01131     }
01132   }
01133 
01134 
01135   void BLDPC_Generator::save(const std::string& filename) const
01136   {
01137     it_assert(init_flag,
01138               "BLDPC_Generator::save(): Can not save not initialized generator");
01139     // Every Z-th row of H_enc until M-Z
01140     GF2mat H_T(M/Z-1, N);
01141     for (int i = 0; i < M/Z-1; i++) {
01142       H_T.set_row(i, H_enc.get_row(i*Z));
01143     }
01144     // Last Z preprocessed rows of H_enc
01145     GF2mat H_Z = H_enc.get_submatrix(M-Z, 0, M-1, N-1);
01146 
01147     it_file f(filename);
01148     int ver;
01149     f >> Name("Fileversion") >> ver;
01150     it_assert(ver == LDPC_binary_file_version, "BLDPC_Generator::save(): "
01151               "Unsupported file format");
01152     f << Name("G_type") << type;
01153     f << Name("H_T") << H_T;
01154     f << Name("H_Z") << H_Z;
01155     f << Name("Z") << Z;
01156     f.close();
01157   }
01158 
01159   void BLDPC_Generator::load(const std::string& filename)
01160   {
01161     GF2mat H_T, H_Z;
01162 
01163     it_ifile f(filename);
01164     int ver;
01165     f >> Name("Fileversion") >> ver;
01166     it_assert(ver == LDPC_binary_file_version, "BLDPC_Generator::load(): "
01167               "Unsupported file format");
01168     std::string gen_type;
01169     f >> Name("G_type") >> gen_type;
01170     it_assert(gen_type == type,
01171               "BLDPC_Generator::load(): Wrong generator type");
01172     f >> Name("H_T") >> H_T;
01173     f >> Name("H_Z") >> H_Z;
01174     f >> Name("Z") >> Z;
01175     f.close();
01176 
01177     N = H_T.cols();
01178     M = (H_T.rows() + 1) * Z;
01179     K = N-M;
01180 
01181     H_enc = GF2mat(M-Z, N);
01182     for (int i = 0; i < H_T.rows(); i++) {
01183       for (int j = 0; j < Z; j++) {
01184         for (int k = 0; k < N; k++) {
01185           if (H_T(i, (k/Z)*Z + (k+Z-j)%Z)) {
01186             H_enc.set(i*Z + j, k, 1);
01187           }
01188         }
01189       }
01190     }
01191     H_enc = H_enc.concatenate_vertical(H_Z);
01192 
01193     init_flag = true;
01194   }
01195 
01196 
01197   // ----------------------------------------------------------------------
01198   // LDPC_Code
01199   // ----------------------------------------------------------------------
01200 
01201   LDPC_Code::LDPC_Code(): H_defined(false), G_defined(false), dec_method("BP"),
01202                           max_iters(50), psc(true), pisc(false),
01203                           llrcalc(LLR_calc_unit()) { }
01204 
01205   LDPC_Code::LDPC_Code(const LDPC_Parity* const H,
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     set_code(H, G_in);
01211   }
01212 
01213   LDPC_Code::LDPC_Code(const std::string& filename,
01214                        LDPC_Generator* const G_in):
01215     H_defined(false), G_defined(false), dec_method("BP"), max_iters(50),
01216     psc(true), pisc(false), llrcalc(LLR_calc_unit())
01217   {
01218     load_code(filename, G_in);
01219   }
01220 
01221 
01222   void LDPC_Code::set_code(const LDPC_Parity* const H,
01223                            LDPC_Generator* const G_in)
01224   {
01225     decoder_parameterization(H);
01226     setup_decoder();
01227     G = G_in;
01228     if (G != 0) {
01229       G_defined = true;
01230       integrity_check();
01231     }
01232   }
01233 
01234   void LDPC_Code::load_code(const std::string& filename,
01235                             LDPC_Generator* const G_in)
01236   {
01237     it_info_debug("LDPC_Code::load_code(): Loading LDPC codec from "
01238                   << filename);
01239 
01240     it_ifile f(filename);
01241     int ver;
01242     f >> Name("Fileversion") >> ver;
01243     it_assert(ver == LDPC_binary_file_version,"LDPC_Code::load_code(): "
01244               "Unsupported file format");
01245     f >> Name("H_defined") >> H_defined;
01246     f >> Name("G_defined") >> G_defined;
01247     f >> Name("nvar") >> nvar;
01248     f >> Name("ncheck") >> ncheck;
01249     f >> Name("C") >> C;
01250     f >> Name("V") >> V;
01251     f >> Name("sumX1") >> sumX1;
01252     f >> Name("sumX2") >> sumX2;
01253     f >> Name("iind") >> iind;
01254     f >> Name("jind") >> jind;
01255     f.close();
01256 
01257     // load generator data
01258     if (G_defined) {
01259       it_assert(G_in != 0, "LDPC_Code::load_code(): Generator object is "
01260                 "missing. Can not load the generator data from a file.");
01261       G = G_in;
01262       G->load(filename);
01263     }
01264     else {
01265       G = 0;
01266       it_info_debug("LDPC_Code::load_code(): Generator data not loaded. "
01267                     "Generator object will not be used.");
01268     }
01269 
01270     it_info_debug("LDPC_Code::load_code(): Successfully loaded LDPC codec "
01271                   "from " << filename);
01272 
01273     setup_decoder();
01274   }
01275 
01276   void LDPC_Code::save_code(const std::string& filename) const
01277   {
01278     it_assert(H_defined, "LDPC_Code::save_to_file(): There is no parity "
01279               "check matrix");
01280     it_info_debug("LDPC_Code::save_to_file(): Saving LDPC codec to "
01281                   << filename);
01282 
01283     it_file f;
01284     f.open(filename,true);
01285     f << Name("Fileversion") << LDPC_binary_file_version;
01286     f << Name("H_defined") << H_defined;
01287     f << Name("G_defined") << G_defined;
01288     f << Name("nvar") << nvar;
01289     f << Name("ncheck") << ncheck;
01290     f << Name("C") << C;
01291     f << Name("V") << V;
01292     f << Name("sumX1") << sumX1;
01293     f << Name("sumX2") << sumX2;
01294     f << Name("iind") << iind;
01295     f << Name("jind") << jind;
01296     f.close();
01297 
01298     // save generator data;
01299     if (G_defined)
01300       G->save(filename);
01301     else
01302       it_info_debug("LDPC_Code::save_code(): Missing generator object - "
01303                     "generator data not saved");
01304 
01305     it_info_debug("LDPC_Code::save_code(): Successfully saved LDPC codec to "
01306                   << filename);
01307   }
01308 
01309 
01310   void LDPC_Code::set_decoding_method(const std::string& method_in)
01311   {
01312     it_assert((method_in == "bp") || (method_in == "BP"),
01313               "LDPC_Code::set_decoding_method(): Not implemented decoding method");
01314     dec_method = method_in;
01315   }
01316 
01317   void LDPC_Code::set_exit_conditions(int max_iters_in,
01318                                       bool syndr_check_each_iter,
01319                                       bool syndr_check_at_start)
01320   {
01321     it_assert(max_iters >= 0, "LDPC_Code::set_nrof_iterations(): Maximum "
01322               "number of iterations can not be negative");
01323     max_iters = max_iters_in;
01324     psc = syndr_check_each_iter;
01325     pisc = syndr_check_at_start;
01326   }
01327 
01328   void LDPC_Code::set_llrcalc(const LLR_calc_unit& llrcalc_in)
01329   {
01330     llrcalc = llrcalc_in;
01331   }
01332 
01333 
01334   void LDPC_Code::encode(const bvec &input, bvec &output)
01335   {
01336     it_assert(G_defined, "LDPC_Code::encode(): LDPC Generator is required "
01337               "for encoding");
01338     G->encode(input, output);
01339     it_assert_debug(syndrome_check(output), "LDPC_Code::encode(): Syndrome "
01340                     "check failed");
01341   }
01342 
01343   bvec LDPC_Code::encode(const bvec &input)
01344   {
01345     bvec result;
01346     encode(input, result);
01347     return result;
01348   }
01349 
01350   void LDPC_Code::decode(const vec &llr_in, bvec &syst_bits)
01351   {
01352     QLLRvec qllrin = llrcalc.to_qllr(llr_in);
01353     QLLRvec qllrout;
01354     bp_decode(qllrin, qllrout);
01355     syst_bits = (qllrout.left(ncheck) < 0);
01356   }
01357 
01358   bvec LDPC_Code::decode(const vec &llr_in)
01359   {
01360     bvec syst_bits;
01361     decode(llr_in, syst_bits);
01362     return syst_bits;
01363   }
01364 
01365   void LDPC_Code::decode_soft_out(const vec &llr_in, vec &llr_out)
01366   {
01367     QLLRvec qllrin = llrcalc.to_qllr(llr_in);
01368     QLLRvec qllrout;
01369     bp_decode(qllrin, qllrout);
01370     llr_out = llrcalc.to_double(qllrout);
01371   }
01372 
01373   vec LDPC_Code::decode_soft_out(const vec &llr_in)
01374   {
01375     vec llr_out;
01376     decode_soft_out(llr_in, llr_out);
01377     return llr_out;
01378   }
01379 
01380   int LDPC_Code::bp_decode(const QLLRvec &LLRin, QLLRvec &LLRout)
01381   {
01382     // Note the IT++ convention that a sure zero corresponds to
01383     // LLR=+infinity
01384     it_assert(H_defined, "LDPC_Code::bp_decode(): Parity check matrix not "
01385               "defined");
01386     it_assert((LLRin.size() == nvar) && (sumX1.size() == nvar)
01387               && (sumX2.size() == ncheck), "LDPC_Code::bp_decode(): Wrong "
01388               "input dimensions");
01389 
01390     if (pisc && syndrome_check(LLRin)) {
01391       LLRout = LLRin;
01392       return 0;
01393     }
01394 
01395     LLRout.set_size(LLRin.size());
01396 
01397     // initial step
01398     for (int i=0; i<nvar; i++) {
01399       int index = i;
01400       for (int j=0; j<sumX1(i); j++) {
01401         mvc[index] = LLRin(i);
01402         index += nvar;
01403       }
01404     }
01405 
01406     bool is_valid_codeword=false;
01407     int iter =0;
01408     do {
01409       iter++;
01410       if (nvar>=100000) { it_info_no_endl_debug("."); }
01411       // --------- Step 1: check to variable nodes ----------
01412       for (int j=0; j<ncheck; j++) {
01413         switch (sumX2(j)) {
01414         case 0: it_error("LDPC_Code::bp_decode(): sumX2(j)=0");
01415         case 1: it_error("LDPC_Code::bp_decode(): sumX2(j)=1");
01416         case 2: {
01417           mcv[j+ncheck]=mvc[jind[j]];
01418           mcv[j]=mvc[jind[j+ncheck]];
01419           break;
01420         }
01421         case 3: {
01422           int j0=j;
01423           QLLR m0=mvc[jind[j0]];
01424           int j1=j0+ncheck;
01425           QLLR m1=mvc[jind[j1]];
01426           int j2=j1+ncheck;
01427           QLLR m2=mvc[jind[j2]];
01428           mcv[j0]=llrcalc.Boxplus(m1,m2);
01429           mcv[j1]=llrcalc.Boxplus(m0,m2);
01430           mcv[j2]=llrcalc.Boxplus(m0,m1);
01431           break;
01432         }
01433         case 4: {
01434           int j0=j;
01435           QLLR m0=mvc[jind[j0]];
01436           int j1=j0+ncheck;
01437           QLLR m1=mvc[jind[j1]];
01438           int j2=j1+ncheck;
01439           QLLR m2=mvc[jind[j2]];
01440           int j3=j2+ncheck;
01441           QLLR m3=mvc[jind[j3]];
01442           QLLR m01=llrcalc.Boxplus(m0,m1);
01443           QLLR m23=llrcalc.Boxplus(m2,m3);
01444           mcv[j0]=llrcalc.Boxplus(m1,m23);
01445           mcv[j1]=llrcalc.Boxplus(m0,m23);
01446           mcv[j2]=llrcalc.Boxplus(m01,m3);
01447           mcv[j3]=llrcalc.Boxplus(m01,m2);
01448           break;
01449         }
01450         case 5: {
01451           int j0=j;
01452           QLLR m0=mvc[jind[j0]];
01453           int j1=j0+ncheck;
01454           QLLR m1=mvc[jind[j1]];
01455           int j2=j1+ncheck;
01456           QLLR m2=mvc[jind[j2]];
01457           int j3=j2+ncheck;
01458           QLLR m3=mvc[jind[j3]];
01459           int j4=j3+ncheck;
01460           QLLR m4=mvc[jind[j4]];
01461           QLLR m01=llrcalc.Boxplus(m0,m1);
01462           QLLR m02=llrcalc.Boxplus(m01,m2);
01463           QLLR m34=llrcalc.Boxplus(m3,m4);
01464           QLLR m24=llrcalc.Boxplus(m2,m34);
01465           mcv[j0]=llrcalc.Boxplus(m1,m24);
01466           mcv[j1]=llrcalc.Boxplus(m0,m24);
01467           mcv[j2]=llrcalc.Boxplus(m01,m34);
01468           mcv[j3]=llrcalc.Boxplus(m02,m4);
01469           mcv[j4]=llrcalc.Boxplus(m02,m3);
01470           break;
01471         }
01472         case 6: {
01473           int j0=j;
01474           QLLR m0=mvc[jind[j0]];
01475           int j1=j0+ncheck;
01476           QLLR m1=mvc[jind[j1]];
01477           int j2=j1+ncheck;
01478           QLLR m2=mvc[jind[j2]];
01479           int j3=j2+ncheck;
01480           QLLR m3=mvc[jind[j3]];
01481           int j4=j3+ncheck;
01482           QLLR m4=mvc[jind[j4]];
01483           int j5=j4+ncheck;
01484           QLLR m5=mvc[jind[j5]];
01485           QLLR m01=llrcalc.Boxplus(m0,m1);
01486           QLLR m23=llrcalc.Boxplus(m2,m3);
01487           QLLR m45=llrcalc.Boxplus(m4,m5);
01488           QLLR m03=llrcalc.Boxplus(m01,m23);
01489           QLLR m25=llrcalc.Boxplus(m23,m45);
01490           QLLR m0145=llrcalc.Boxplus(m01,m45);
01491           mcv[j0]=llrcalc.Boxplus(m1,m25);
01492           mcv[j1]=llrcalc.Boxplus(m0,m25);
01493           mcv[j2]=llrcalc.Boxplus(m0145,m3);
01494           mcv[j3]=llrcalc.Boxplus(m0145,m2);
01495           mcv[j4]=llrcalc.Boxplus(m03,m5);
01496           mcv[j5]=llrcalc.Boxplus(m03,m4);
01497           break;
01498         }
01499         case 7: {
01500           int j0=j;
01501           QLLR m0=mvc[jind[j0]];
01502           int j1=j0+ncheck;
01503           QLLR m1=mvc[jind[j1]];
01504           int j2=j1+ncheck;
01505           QLLR m2=mvc[jind[j2]];
01506           int j3=j2+ncheck;
01507           QLLR m3=mvc[jind[j3]];
01508           int j4=j3+ncheck;
01509           QLLR m4=mvc[jind[j4]];
01510           int j5=j4+ncheck;
01511           QLLR m5=mvc[jind[j5]];
01512           int j6=j5+ncheck;
01513           QLLR m6=mvc[jind[j6]];
01514           QLLR m01=llrcalc.Boxplus(m0,m1);
01515           QLLR m23=llrcalc.Boxplus(m2,m3);
01516           QLLR m45=llrcalc.Boxplus(m4,m5);
01517           QLLR m46=llrcalc.Boxplus(m45,m6);
01518           QLLR m03=llrcalc.Boxplus(m01,m23);
01519           QLLR m25=llrcalc.Boxplus(m23,m45);
01520           QLLR m26=llrcalc.Boxplus(m25,m6);
01521           QLLR m04=llrcalc.Boxplus(m03,m4);
01522           mcv[j0]=llrcalc.Boxplus(m26,m1);
01523           mcv[j1]=llrcalc.Boxplus(m26,m0);
01524           mcv[j2]=llrcalc.Boxplus(m01,llrcalc.Boxplus(m3,m46));
01525           mcv[j3]=llrcalc.Boxplus(m2,llrcalc.Boxplus(m01,m46));
01526           mcv[j4]=llrcalc.Boxplus(m6,llrcalc.Boxplus(m03,m5));
01527           mcv[j5]=llrcalc.Boxplus(m6,m04);
01528           mcv[j6]=llrcalc.Boxplus(m5,m04);
01529           break;
01530         }
01531         case 8: {
01532           int j0=j;
01533           QLLR m0=mvc[jind[j0]];
01534           int j1=j0+ncheck;
01535           QLLR m1=mvc[jind[j1]];
01536           int j2=j1+ncheck;
01537           QLLR m2=mvc[jind[j2]];
01538           int j3=j2+ncheck;
01539           QLLR m3=mvc[jind[j3]];
01540           int j4=j3+ncheck;
01541           QLLR m4=mvc[jind[j4]];
01542           int j5=j4+ncheck;
01543           QLLR m5=mvc[jind[j5]];
01544           int j6=j5+ncheck;
01545           QLLR m6=mvc[jind[j6]];
01546           int j7=j6+ncheck;
01547           QLLR m7=mvc[jind[j7]];
01548           QLLR m01=llrcalc.Boxplus(m0,m1);
01549           QLLR m23=llrcalc.Boxplus(m2,m3);
01550           QLLR m45=llrcalc.Boxplus(m4,m5);
01551           QLLR m67=llrcalc.Boxplus(m6,m7);
01552           QLLR m47=llrcalc.Boxplus(m45,m67);
01553           QLLR m03=llrcalc.Boxplus(m01,m23);
01554           QLLR m25=llrcalc.Boxplus(m23,m45);
01555           mcv[j0]=llrcalc.Boxplus(m67,llrcalc.Boxplus(m1,m25));
01556           mcv[j1]=llrcalc.Boxplus(m67,llrcalc.Boxplus(m0,m25));
01557           mcv[j2]=llrcalc.Boxplus(m3,llrcalc.Boxplus(m01,m47));
01558           mcv[j3]=llrcalc.Boxplus(m2,llrcalc.Boxplus(m01,m47));
01559           mcv[j4]=llrcalc.Boxplus(m67,llrcalc.Boxplus(m03,m5));
01560           mcv[j5]=llrcalc.Boxplus(m67,llrcalc.Boxplus(m03,m4));
01561           mcv[j6]=llrcalc.Boxplus(m45,llrcalc.Boxplus(m03,m7));
01562           mcv[j7]=llrcalc.Boxplus(m03,llrcalc.Boxplus(m45,m6));
01563           break;
01564         }
01565         case 9: {
01566           int j0=j;
01567           QLLR m0=mvc[jind[j0]];
01568           int j1=j0+ncheck;
01569           QLLR m1=mvc[jind[j1]];
01570           int j2=j1+ncheck;
01571           QLLR m2=mvc[jind[j2]];
01572           int j3=j2+ncheck;
01573           QLLR m3=mvc[jind[j3]];
01574           int j4=j3+ncheck;
01575           QLLR m4=mvc[jind[j4]];
01576           int j5=j4+ncheck;
01577           QLLR m5=mvc[jind[j5]];
01578           int j6=j5+ncheck;
01579           QLLR m6=mvc[jind[j6]];
01580           int j7=j6+ncheck;
01581           QLLR m7=mvc[jind[j7]];
01582           int j8=j7+ncheck;
01583           QLLR m8=mvc[jind[j8]];
01584           QLLR m01=llrcalc.Boxplus(m0,m1);
01585           QLLR m23=llrcalc.Boxplus(m2,m3);
01586           QLLR m45=llrcalc.Boxplus(m4,m5);
01587           QLLR m67=llrcalc.Boxplus(m6,m7);
01588           QLLR m68=llrcalc.Boxplus(m67,m8);
01589           QLLR m03=llrcalc.Boxplus(m01,m23);
01590           QLLR m25=llrcalc.Boxplus(m23,m45);
01591           QLLR m05=llrcalc.Boxplus(m03,m45);
01592           mcv[j0]=llrcalc.Boxplus(m68,llrcalc.Boxplus(m1,m25));
01593           mcv[j1]=llrcalc.Boxplus(m68,llrcalc.Boxplus(m0,m25));
01594           mcv[j2]=llrcalc.Boxplus(llrcalc.Boxplus(m01,m68),
01595                                   llrcalc.Boxplus(m3,m45));
01596           mcv[j3]=llrcalc.Boxplus(llrcalc.Boxplus(m01,m68),
01597                                   llrcalc.Boxplus(m2,m45));
01598           mcv[j4]=llrcalc.Boxplus(m68,llrcalc.Boxplus(m03,m5));
01599           mcv[j5]=llrcalc.Boxplus(m68,llrcalc.Boxplus(m03,m4));
01600           mcv[j6]=llrcalc.Boxplus(llrcalc.Boxplus(m7,m8),m05);
01601           mcv[j7]=llrcalc.Boxplus(llrcalc.Boxplus(m05,m6),m8);
01602           mcv[j8]=llrcalc.Boxplus(m05,m67);
01603           break;
01604         }
01605         case 10: {
01606           int j0=j;
01607           QLLR m0=mvc[jind[j0]];
01608           int j1=j0+ncheck;
01609           QLLR m1=mvc[jind[j1]];
01610           int j2=j1+ncheck;
01611           QLLR m2=mvc[jind[j2]];
01612           int j3=j2+ncheck;
01613           QLLR m3=mvc[jind[j3]];
01614           int j4=j3+ncheck;
01615           QLLR m4=mvc[jind[j4]];
01616           int j5=j4+ncheck;
01617           QLLR m5=mvc[jind[j5]];
01618           int j6=j5+ncheck;
01619           QLLR m6=mvc[jind[j6]];
01620           int j7=j6+ncheck;
01621           QLLR m7=mvc[jind[j7]];
01622           int j8=j7+ncheck;
01623           QLLR m8=mvc[jind[j8]];
01624           int j9=j8+ncheck;
01625           QLLR m9=mvc[jind[j9]];
01626           QLLR m01=llrcalc.Boxplus(m0,m1);
01627           QLLR m23=llrcalc.Boxplus(m2,m3);
01628           QLLR m03=llrcalc.Boxplus(m01,m23);
01629           QLLR m45=llrcalc.Boxplus(m4,m5);
01630           QLLR m67=llrcalc.Boxplus(m6,m7);
01631           QLLR m89=llrcalc.Boxplus(m8,m9);
01632           QLLR m69=llrcalc.Boxplus(m67,m89);
01633           QLLR m25=llrcalc.Boxplus(m23,m45);
01634           QLLR m05=llrcalc.Boxplus(m03,m45);
01635           QLLR m07=llrcalc.Boxplus(m05,m67);
01636           mcv[j0]=llrcalc.Boxplus(m69,llrcalc.Boxplus(m1,m25));
01637           mcv[j1]=llrcalc.Boxplus(m69,llrcalc.Boxplus(m0,m25));
01638           mcv[j2]=llrcalc.Boxplus(llrcalc.Boxplus(m01,m69),
01639                                   llrcalc.Boxplus(m3,m45));
01640           mcv[j3]=llrcalc.Boxplus(llrcalc.Boxplus(m01,m69),
01641                                   llrcalc.Boxplus(m2,m45));
01642           mcv[j4]=llrcalc.Boxplus(m69,llrcalc.Boxplus(m03,m5));
01643           mcv[j5]=llrcalc.Boxplus(m69,llrcalc.Boxplus(m03,m4));
01644           mcv[j6]=llrcalc.Boxplus(llrcalc.Boxplus(m7,m89),m05);
01645           mcv[j7]=llrcalc.Boxplus(llrcalc.Boxplus(m05,m6),m89);
01646           mcv[j8]=llrcalc.Boxplus(m07,m9);
01647           mcv[j9]=llrcalc.Boxplus(m07,m8);
01648           break;
01649         }
01650         case 11: {
01651           int j0=j;
01652           QLLR m0=mvc[jind[j0]];
01653           int j1=j0+ncheck;
01654           QLLR m1=mvc[jind[j1]];
01655           int j2=j1+ncheck;
01656           QLLR m2=mvc[jind[j2]];
01657           int j3=j2+ncheck;
01658           QLLR m3=mvc[jind[j3]];
01659           int j4=j3+ncheck;
01660           QLLR m4=mvc[jind[j4]];
01661           int j5=j4+ncheck;
01662           QLLR m5=mvc[jind[j5]];
01663           int j6=j5+ncheck;
01664           QLLR m6=mvc[jind[j6]];
01665           int j7=j6+ncheck;
01666           QLLR m7=mvc[jind[j7]];
01667           int j8=j7+ncheck;
01668           QLLR m8=mvc[jind[j8]];
01669           int j9=j8+ncheck;
01670           QLLR m9=mvc[jind[j9]];
01671           int j10=j9+ncheck;
01672           QLLR m10=mvc[jind[j10]];
01673           QLLR m01=llrcalc.Boxplus(m0,m1);
01674           QLLR m23=llrcalc.Boxplus(m2,m3);
01675           QLLR m03=llrcalc.Boxplus(m01,m23);
01676           QLLR m45=llrcalc.Boxplus(m4,m5);
01677           QLLR m67=llrcalc.Boxplus(m6,m7);
01678           QLLR m89=llrcalc.Boxplus(m8,m9);
01679           QLLR m69=llrcalc.Boxplus(m67,m89);
01680           QLLR m6_10=llrcalc.Boxplus(m69,m10);
01681           QLLR m25=llrcalc.Boxplus(m23,m45);
01682           QLLR m05=llrcalc.Boxplus(m03,m45);
01683           QLLR m07=llrcalc.Boxplus(m05,m67);
01684           QLLR m8_10=llrcalc.Boxplus(m89,m10);
01685           mcv[j0]=llrcalc.Boxplus(m6_10,llrcalc.Boxplus(m1,m25));
01686           mcv[j1]=llrcalc.Boxplus(m6_10,llrcalc.Boxplus(m0,m25));
01687           mcv[j2]=llrcalc.Boxplus(llrcalc.Boxplus(m01,m6_10),
01688                                   llrcalc.Boxplus(m3,m45));
01689           mcv[j3]=llrcalc.Boxplus(llrcalc.Boxplus(m01,m6_10),
01690                                   llrcalc.Boxplus(m2,m45));
01691           mcv[j4]=llrcalc.Boxplus(m6_10,llrcalc.Boxplus(m03,m5));
01692           mcv[j5]=llrcalc.Boxplus(m6_10,llrcalc.Boxplus(m03,m4));
01693           mcv[j6]=llrcalc.Boxplus(llrcalc.Boxplus(m7,m8_10),m05);
01694           mcv[j7]=llrcalc.Boxplus(llrcalc.Boxplus(m05,m6),m8_10);
01695           mcv[j8]=llrcalc.Boxplus(m10,llrcalc.Boxplus(m07,m9));
01696           mcv[j9]=llrcalc.Boxplus(m10,llrcalc.Boxplus(m07,m8));
01697           mcv[j10]=llrcalc.Boxplus(m07,m89);
01698           break;
01699         }
01700         case 12: {
01701           int j0=j;
01702           QLLR m0=mvc[jind[j0]];
01703           int j1=j0+ncheck;
01704           QLLR m1=mvc[jind[j1]];
01705           int j2=j1+ncheck;
01706           QLLR m2=mvc[jind[j2]];
01707           int j3=j2+ncheck;
01708           QLLR m3=mvc[jind[j3]];
01709           int j4=j3+ncheck;
01710           QLLR m4=mvc[jind[j4]];
01711           int j5=j4+ncheck;
01712           QLLR m5=mvc[jind[j5]];
01713           int j6=j5+ncheck;
01714           QLLR m6=mvc[jind[j6]];
01715           int j7=j6+ncheck;
01716           QLLR m7=mvc[jind[j7]];
01717           int j8=j7+ncheck;
01718           QLLR m8=mvc[jind[j8]];
01719           int j9=j8+ncheck;
01720           QLLR m9=mvc[jind[j9]];
01721           int j10=j9+ncheck;
01722           QLLR m10=mvc[jind[j10]];
01723           int j11=j10+ncheck;
01724           QLLR m11=mvc[jind[j11]];
01725           QLLR m01=llrcalc.Boxplus(m0,m1);
01726           QLLR m23=llrcalc.Boxplus(m2,m3);
01727           QLLR m03=llrcalc.Boxplus(m01,m23);
01728           QLLR m45=llrcalc.Boxplus(m4,m5);
01729           QLLR m67=llrcalc.Boxplus(m6,m7);
01730           QLLR m89=llrcalc.Boxplus(m8,m9);
01731           QLLR m69=llrcalc.Boxplus(m67,m89);
01732           QLLR m10_11=llrcalc.Boxplus(m10,m11);
01733           QLLR m6_11=llrcalc.Boxplus(m69,m10_11);
01734           QLLR m25=llrcalc.Boxplus(m23,m45);
01735           QLLR m05=llrcalc.Boxplus(m03,m45);
01736           QLLR m07=llrcalc.Boxplus(m05,m67);
01737           QLLR m8_10=llrcalc.Boxplus(m89,m10);
01738           mcv[j0]=llrcalc.Boxplus(m6_11,llrcalc.Boxplus(m1,m25));
01739           mcv[j1]=llrcalc.Boxplus(m6_11,llrcalc.Boxplus(m0,m25));
01740           mcv[j2]=llrcalc.Boxplus(llrcalc.Boxplus(m01,m6_11),
01741                                   llrcalc.Boxplus(m3,m45));
01742           mcv[j3]=llrcalc.Boxplus(llrcalc.Boxplus(m01,m6_11),
01743                                   llrcalc.Boxplus(m2,m45));
01744           mcv[j4]=llrcalc.Boxplus(m6_11,llrcalc.Boxplus(m03,m5));
01745           mcv[j5]=llrcalc.Boxplus(m6_11,llrcalc.Boxplus(m03,m4));
01746           mcv[j6]=llrcalc.Boxplus(m11,llrcalc.Boxplus(llrcalc.Boxplus(m7,m8_10),m05));
01747           mcv[j7]=llrcalc.Boxplus(m11,llrcalc.Boxplus(llrcalc.Boxplus(m05,m6),m8_10));
01748           mcv[j8]=llrcalc.Boxplus(m10_11,llrcalc.Boxplus(m07,m9));
01749           mcv[j9]=llrcalc.Boxplus(m10_11,llrcalc.Boxplus(m07,m8));
01750           mcv[j10]=llrcalc.Boxplus(llrcalc.Boxplus(m07,m89),m11);
01751           mcv[j11]=llrcalc.Boxplus(llrcalc.Boxplus(m07,m89),m10);
01752           break;
01753         }
01754         default:
01755           it_error("check node degrees >12 not supported in this version");
01756         }  // switch statement
01757       }
01758 
01759       // step 2: variable to check nodes
01760       for (int i=0; i<nvar; i++) {
01761         switch (sumX1(i)) {
01762         case 0: it_error("LDPC_Code::bp_decode(): sumX1(i)=0");
01763         case 1: {
01764           // Degenerate case-should not occur for good coded. A lonely
01765           // variable node only sends its incoming message
01766           mvc[i] = LLRin(i);
01767           LLRout(i)=LLRin(i);
01768           break;
01769         }
01770         case 2: {
01771           QLLR m0=mcv[iind[i]];
01772           int i1=i+nvar;
01773           QLLR m1=mcv[iind[i1]];
01774           mvc[i] = LLRin(i) + m1;
01775           mvc[i1] = LLRin(i) + m0;
01776           LLRout(i) = mvc[i1]+m1;
01777           break;
01778         }
01779         case 3: {
01780           int i0=i;
01781           QLLR m0 = mcv[iind[i0]];
01782           int i1 = i0+nvar;
01783           QLLR m1 = mcv[iind[i1]];
01784           int i2 = i1+nvar;
01785           QLLR m2 = mcv[iind[i2]];
01786           LLRout(i) = LLRin(i)+m0+m1+m2;
01787           mvc[i0]=LLRout(i)-m0;
01788           mvc[i1]=LLRout(i)-m1;
01789           mvc[i2]=LLRout(i)-m2;
01790           break;
01791         }
01792         case 4: {
01793           int i0=i;
01794           QLLR m0 = mcv[iind[i0]];
01795           int i1 = i0+nvar;
01796           QLLR m1 = mcv[iind[i1]];
01797           int i2 = i1+nvar;
01798           QLLR m2 = mcv[iind[i2]];
01799           int i3 = i2+nvar;
01800           QLLR m3 = mcv[iind[i3]];
01801           LLRout(i)= LLRin(i)+m0+m1+m2+m3;
01802           mvc[i0]=LLRout(i)-m0;
01803           mvc[i1]=LLRout(i)-m1;
01804           mvc[i2]=LLRout(i)-m2;
01805           mvc[i3]=LLRout(i)-m3;
01806           break;
01807         }
01808         default:        { // differential update
01809           QLLR mvc_temp = LLRin(i);
01810           int index_iind = i; // tracks i+jp*nvar
01811           for (int jp=0; jp<sumX1(i); jp++) {
01812             mvc_temp +=  mcv[iind[index_iind]];
01813             index_iind += nvar;
01814           }
01815           LLRout(i) = mvc_temp;
01816           index_iind = i;  // tracks i+j*nvar
01817           for (int j=0; j<sumX1[i]; j++) {
01818             mvc[index_iind] = mvc_temp - mcv[iind[index_iind]];
01819             index_iind += nvar;
01820           }
01821         }
01822         }
01823       }
01824 
01825       if (psc && syndrome_check(LLRout)) {
01826         is_valid_codeword=true;
01827         break;
01828       }
01829     } while  (iter < max_iters);
01830 
01831     if (nvar>=100000) { it_info_debug(""); }
01832     return (is_valid_codeword ? iter : -iter);
01833   }
01834 
01835 
01836   bool LDPC_Code::syndrome_check(const bvec &x) const
01837   {
01838     QLLRvec llr=1-2*to_ivec(x);
01839     return syndrome_check(llr);
01840   }
01841 
01842   bool LDPC_Code::syndrome_check(const QLLRvec &LLR) const
01843   {
01844     // Please note the IT++ convention that a sure zero corresponds to
01845     // LLR=+infinity
01846     int i,j,synd,vi;
01847 
01848     for (j=0; j<ncheck; j++) {
01849       synd = 0;
01850       int vind = j; // tracks j+i*ncheck
01851       for (i=0; i<sumX2(j); i++) {
01852         vi = V(vind);
01853         if (LLR(vi)<0) {
01854           synd++;
01855         }
01856         vind += ncheck;
01857       }
01858       if ((synd&1)==1) {
01859         return false;  // codeword is invalid
01860       }
01861     }
01862     return true;   // codeword is valid
01863   };
01864 
01865 
01866   // ----------------------------------------------------------------------
01867   // LDPC_Code private methods
01868   // ----------------------------------------------------------------------
01869 
01870   void LDPC_Code::decoder_parameterization(const LDPC_Parity* const Hmat)
01871   {
01872     // copy basic parameters
01873     sumX1 = Hmat->sumX1;
01874     sumX2 = Hmat->sumX2;
01875     nvar = Hmat->nvar; //get_nvar();
01876     ncheck = Hmat->ncheck; //get_ncheck();
01877     int cmax = max(sumX1);
01878     int vmax = max(sumX2);
01879 
01880     // decoder parameterization
01881     V = zeros_i(ncheck*vmax);
01882     C = zeros_i(cmax*nvar);
01883     jind = zeros_i(ncheck*vmax);
01884     iind = zeros_i(nvar*cmax);
01885 
01886     it_info_debug("LDPC_Code::decoder_parameterization(): Computations "
01887                   "- phase 1");
01888     for (int i=0; i<nvar; i++) {
01889       ivec coli = Hmat->get_col(i).get_nz_indices();
01890       for (int j0=0; j0<length(coli); j0++) {
01891         C(j0+cmax*i) = coli(j0);
01892       }
01893     }
01894 
01895     it_info_debug("LDPC_Code::decoder_parameterization(): Computations "
01896                   "- phase 2");
01897     it_info_debug("Computing decoder parameterization. Phase 2");
01898     for (int j=0; j<ncheck; j++) {
01899       ivec rowj = Hmat->get_row(j).get_nz_indices();
01900       for (int i0=0; i0<length(rowj); i0++) {
01901         V(j+ncheck*i0) = rowj(i0);
01902       }
01903     }
01904 
01905     it_info_debug("LDPC_Code::decoder_parameterization(): Computations "
01906                   "- phase 3");
01907     it_info_debug("Computing decoder parameterization. Phase 3.");
01908     for (int j=0; j<ncheck; j++) {
01909       for (int ip=0; ip<sumX2(j); ip++) {
01910         int vip = V(j+ip*ncheck);
01911         int k=0;
01912         while (1==1)  {
01913           if (C(k+vip*cmax)==j) {
01914             break;
01915           }
01916           k++;
01917         }
01918         jind(j+ip*ncheck) = vip+k*nvar;
01919       }
01920     }
01921 
01922     it_info_debug("LDPC_Code::decoder_parameterization(): Computations "
01923                   "- phase 4");
01924     for (int i=0; i<nvar; i++) {
01925       for (int jp=0; jp<sumX1(i); jp++) {
01926         int cjp = C(jp+i*cmax);
01927         int k=0;
01928         while (1==1) {
01929           if (V(cjp+k*ncheck)==i) {break; }
01930           k++;
01931         }
01932         iind(i+jp*nvar) = cjp+k*ncheck;
01933       }
01934     }
01935 
01936     H_defined = true;
01937   }
01938 
01939 
01940   void LDPC_Code::setup_decoder()
01941   {
01942     if (H_defined) {
01943       mcv.set_size(max(sumX2) * ncheck);
01944       mvc.set_size(max(sumX1) * nvar);
01945     }
01946   }
01947 
01948 
01949   void LDPC_Code::integrity_check()
01950   {
01951     if (G_defined) {
01952       it_info_debug("LDPC_Code::integrity_check(): Checking integrity of "
01953                     "the LDPC_Parity and LDPC_Generator data");
01954       bvec bv(nvar-ncheck), cw;
01955       bv.clear();
01956       bv(0) = 1;
01957       for (int i = 0; i < nvar-ncheck; i++) {
01958         G->encode(bv, cw);
01959         it_assert(syndrome_check(cw),
01960                   "LDPC_Code::integrity_check(): Syndrome check failed");
01961         bv.shift_right(bv(nvar-ncheck-1));
01962       }
01963     }
01964     else {
01965       it_info_debug("LDPC_Code::integrity_check(): No generator defined "
01966                     "- no check performed");
01967     }
01968   }
01969 
01970   // ----------------------------------------------------------------------
01971   // Related functions
01972   // ----------------------------------------------------------------------
01973 
01974   std::ostream &operator<<(std::ostream &os, const LDPC_Code &C)
01975   {
01976     ivec rdeg = zeros_i(max(C.sumX2)+1);
01977     for (int i=0; i<C.ncheck; i++)     {
01978       rdeg(C.sumX2(i))++;
01979     }
01980 
01981     ivec cdeg = zeros_i(max(C.sumX1)+1);
01982     for (int j=0; j<C.nvar; j++)     {
01983       cdeg(C.sumX1(j))++;
01984     }
01985 
01986     os << "--- LDPC codec ----------------------------------\n"
01987        << "Nvar : " << C.get_nvar() << "\n"
01988        << "Ncheck : " << C.get_ncheck() << "\n"
01989        << "Rate : " << C.get_rate() << "\n"
01990        << "Column degrees (node perspective): " << cdeg << "\n"
01991        << "Row degrees (node perspective): " << rdeg << "\n"
01992        << "-------------------------------------------------\n"
01993        << "Decoder parameters:\n"
01994        << " - method : " << C.dec_method << "\n"
01995        << " - max. iterations : " << C.max_iters << "\n"
01996        << " - syndrome check at each iteration : " << C.psc << "\n"
01997        << " - syndrome check at start : " << C.pisc << "\n"
01998        << "-------------------------------------------------\n"
01999        << C.llrcalc << "\n";
02000     return os;
02001   }
02002 
02003 } // namespace itpp
SourceForge Logo

Generated on Mon Jan 7 22:28:57 2008 for IT++ by Doxygen 1.5.4