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