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