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