IT++ Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
ldpc.cpp
Go to the documentation of this file.
1 
29 #include <itpp/comm/ldpc.h>
30 #include <iomanip>
31 #include <sstream>
32 
33 namespace itpp
34 {
35 
42 static const int LDPC_binary_file_version = 2;
43 
44 // ---------------------------------------------------------------------------
45 // LDPC_Parity
46 // ---------------------------------------------------------------------------
47 
48 // public methods
49 
50 LDPC_Parity::LDPC_Parity(int nc, int nv): init_flag(false)
51 {
52  initialize(nc, nv);
53 }
54 
55 LDPC_Parity::LDPC_Parity(const std::string& filename,
56  const std::string& format): init_flag(false)
57 {
58  if (format == "alist") {
59  load_alist(filename);
60  }
61  else {
62  it_error("LDPC_Parity::LDPC_Parity(): Only 'alist' format is supported");
63  }
64 }
65 
67  init_flag(false)
68 {
69  import_alist(alist);
70 }
71 
72 void LDPC_Parity::initialize(int nc, int nv)
73 {
74  ncheck = nc;
75  nvar = nv;
78  sumX1 = zeros_i(nvar);
79  sumX2 = zeros_i(ncheck);
80  init_flag = true;
81 }
82 
83 void LDPC_Parity::set(int i, int j, bin x)
84 {
85  it_assert(init_flag, "LDPC_Parity::set(): Object not initialized");
86  it_assert_debug((i >= 0) && (i < ncheck),
87  "LDPC_Parity::set(): Wrong index i");
88  it_assert_debug((j >= 0) && (j < nvar),
89  "LDPC_Parity::set(): Wrong index j");
90  it_assert_debug(H(i, j) == Ht(j, i), "LDPC_Parity:set(): Internal error");
91 
92  int diff = static_cast<int>(x) - static_cast<int>(H(i, j));
93  sumX1(j) += diff;
94  sumX2(i) += diff;
95 
96  if (x == 1) {
97  H.set(i, j, 1);
98  Ht.set(j, i, 1);
99  }
100  else {
101  H.clear_elem(i, j);
102  Ht.clear_elem(j, i);
103  }
104 
105  it_assert_debug(H(i, j) == x, "LDPC_Parity::set(): Internal error");
106  it_assert_debug(Ht(j, i) == x, "LDPC_Parity::set(): Internal error");
107 }
108 
110 {
112  "LDPC_Parity::display_stats(): Object not initialized");
113  int cmax = max(sumX1);
114  int vmax = max(sumX2);
115  vec vdeg = zeros(cmax + 1); // number of variable nodes with n neighbours
116  vec cdeg = zeros(vmax + 1); // number of check nodes with n neighbours
117  for (int col = 0; col < nvar; col++) {
118  vdeg(length(get_col(col).get_nz_indices()))++;
119  it_assert(sumX1(col) == length(get_col(col).get_nz_indices()),
120  "LDPC_Parity::display_stats(): Internal error");
121  }
122  for (int row = 0; row < ncheck; row++) {
123  cdeg(length(get_row(row).get_nz_indices()))++;
124  it_assert(sumX2(row) == length(get_row(row).get_nz_indices()),
125  "LDPC_Parity::display_stats(): Internal error");
126  }
127 
128  // from edge perspective
129  // number of edges connected to vnodes of degree n
130  vec vdegedge = elem_mult(vdeg, linspace(0, vdeg.length() - 1,
131  vdeg.length()));
132  // number of edges connected to cnodes of degree n
133  vec cdegedge = elem_mult(cdeg, linspace(0, cdeg.length() - 1,
134  cdeg.length()));
135 
136  int edges = sum(elem_mult(to_ivec(linspace(0, vdeg.length() - 1,
137  vdeg.length())),
138  to_ivec(vdeg)));
139 
140  it_info("--- LDPC parity check matrix ---");
141  it_info("Dimension [ncheck x nvar]: " << ncheck << " x " << nvar);
142  it_info("Variable node degree distribution from node perspective:\n"
143  << vdeg / nvar);
144  it_info("Check node degree distribution from node perspective:\n"
145  << cdeg / ncheck);
146  it_info("Variable node degree distribution from edge perspective:\n"
147  << vdegedge / edges);
148  it_info("Check node degree distribution from edge perspective:\n"
149  << cdegedge / edges);
150  it_info("Rate: " << get_rate());
151  it_info("--------------------------------");
152 }
153 
154 
155 void LDPC_Parity::load_alist(const std::string& alist_file)
156 {
157  import_alist(GF2mat_sparse_alist(alist_file));
158 }
159 
160 void LDPC_Parity::save_alist(const std::string& alist_file) const
161 {
163  alist.write(alist_file);
164 }
165 
166 
168 {
169  GF2mat_sparse X = alist.to_sparse();
170 
171  initialize(X.rows(), X.cols());
172  // brute force copy from X to this
173  for (int i = 0; i < ncheck; i++) {
174  for (int j = 0; j < nvar; j++) {
175  if (X(i, j)) {
176  set(i, j, 1);
177  }
178  }
179  }
180 }
181 
183 {
185  "LDPC_Parity::export_alist(): Object not initialized");
186  GF2mat_sparse_alist alist;
187  alist.from_sparse(H);
188  return alist;
189 }
190 
191 
192 int LDPC_Parity::check_connectivity(int from_i, int from_j, int to_i,
193  int to_j, int godir, int L) const
194 {
196  "LDPC_Parity::check_connectivity(): Object not initialized");
197  int i, j, result;
198 
199  if (L < 0) { // unable to reach coordinate with given L
200  return (-3);
201  }
202 
203  // check if reached destination
204  if ((from_i == to_i) && (from_j == to_j) && (godir != 0)) {
205  return L;
206  }
207 
208  if (get(from_i, from_j) == 0) { // meaningless search
209  return (-2);
210  }
211 
212  if (L == 2) { // Treat this case separately for efficiency
213  if (godir == 2) { // go horizontally
214  if (get(from_i, to_j) == 1) { return 0; }
215  }
216  if (godir == 1) { // go vertically
217  if (get(to_i, from_j) == 1) { return 0; }
218  }
219  return (-3);
220  }
221 
222  if ((godir == 1) || (godir == 0)) { // go vertically
223  ivec cj = get_col(from_j).get_nz_indices();
224  for (i = 0; i < length(cj); i++) {
225  if (cj(i) != from_i) {
226  result = check_connectivity(cj(i), from_j, to_i, to_j, 2, L - 1);
227  if (result >= 0) {
228  return (result);
229  }
230  }
231  }
232  }
233 
234  if (godir == 2) { // go horizontally
235  ivec ri = get_row(from_i).get_nz_indices();
236  for (j = 0; j < length(ri); j++) {
237  if (ri(j) != from_j) {
238  result = check_connectivity(from_i, ri(j), to_i, to_j, 1, L - 1);
239  if (result >= 0) {
240  return (result);
241  }
242  }
243  }
244  }
245 
246  return (-1);
247 };
248 
250 {
252  "LDPC_Parity::check_for_cycles(): Object not initialized");
253  // looking for odd length cycles does not make sense
254  if ((L&1) == 1) { return (-1); }
255  if (L == 0) { return (-4); }
256 
257  int cycles = 0;
258  for (int i = 0; i < nvar; i++) {
259  ivec ri = get_col(i).get_nz_indices();
260  for (int j = 0; j < length(ri); j++) {
261  if (check_connectivity(ri(j), i, ri(j), i, 0, L) >= 0) {
262  cycles++;
263  }
264  }
265  }
266  return cycles;
267 };
268 
269 // ivec LDPC_Parity::get_rowdegree() const
270 // {
271 // ivec rdeg = zeros_i(Nmax);
272 // for (int i=0; i<ncheck; i++) {
273 // rdeg(sumX2(i))++;
274 // }
275 // return rdeg;
276 // };
277 
278 // ivec LDPC_Parity::get_coldegree() const
279 // {
280 // ivec cdeg = zeros_i(Nmax);
281 // for (int j=0; j<nvar; j++) {
282 // cdeg(sumX1(j))++;
283 // }
284 // return cdeg;
285 // };
286 
287 
288 // ----------------------------------------------------------------------
289 // LDPC_Parity_Unstructured
290 // ----------------------------------------------------------------------
291 
293 {
295  "LDPC_Parity::cycle_removal_MGW(): Object not initialized");
296  typedef Sparse_Mat<short> Ssmat;
297  typedef Sparse_Vec<short> Ssvec;
298 
299  Maxcyc -= 2;
300 
301  // Construct the adjacency matrix of the graph
302  Ssmat G(ncheck + nvar, ncheck + nvar, 5);
303  for (int j = 0; j < nvar; j++) {
304  GF2vec_sparse col = get_col(j);
305  for (int i = 0; i < col.nnz(); i++) {
306  if (get(col.get_nz_index(i), j) == 1) {
307  G.set(col.get_nz_index(i), j + ncheck, 1);
308  G.set(ncheck + j, col.get_nz_index(i), 1);
309  }
310  }
311  }
312 
313  Array<Ssmat> Gpow(Maxcyc);
314  Gpow(0).set_size(ncheck + nvar, ncheck + nvar, 1);
315  Gpow(0).clear();
316  for (int i = 0; i < ncheck + nvar; i++) {
317  Gpow(0).set(i, i, 1);
318  }
319  Gpow(1) = G;
320 
321  /* Main cycle elimination loop starts here. Note that G and all
322  powers of G are symmetric matrices. This fact is exploited in
323  the code.
324  */
325  int r;
326  int cycles_found = 0;
327  int scl = Maxcyc;
328  for (r = 4; r <= Maxcyc; r += 2) {
329  // compute the next power of the adjacency matrix
330  Gpow(r / 2) = Gpow(r / 2 - 1) * G;
331  bool traverse_again;
332  do {
333  traverse_again = false;
334  cycles_found = 0;
335  it_info_debug("Starting new pass of cycle elimination, target girth "
336  << (r + 2) << "...");
337  int pdone = 0;
338  for (int j = 0; j < ncheck + nvar; j++) { // loop over elements of G
339  for (int i = 0; i < ncheck + nvar; i++) {
340  int ptemp = floor_i(100.0 * (i + j * (ncheck + nvar)) /
341  ((nvar + ncheck) * (nvar + ncheck)));
342  if (ptemp > pdone + 10) {
343  it_info_debug(ptemp << "% done.");
344  pdone = ptemp;
345  }
346 
347  if (((Gpow(r / 2))(i, j) >= 2) && ((Gpow(r / 2 - 2))(i, j) == 0)) {
348  // Found a cycle.
349  cycles_found++;
350 
351  // choose k
352  ivec tmpi = (elem_mult(Gpow(r / 2 - 1).get_col(i),
353  G.get_col(j))).get_nz_indices();
354  // int k = tmpi(rand()%length(tmpi));
355  int k = tmpi(randi(0, length(tmpi) - 1));
356  it_assert_debug(G(j, k) == 1 && G(k, j) == 1,
357  "LDPC_Parity_Unstructured::cycle_removal_MGW(): "
358  "Internal error");
359 
360  // determine candidate edges for an edge swap
361  Ssvec rowjk = Gpow(r / 2) * (Gpow(r / 2 - 1).get_col(j)
362  + Gpow(r / 2 - 1).get_col(k));
363  int p, l;
364  ivec Ce_ind = sort_index(randu(nvar + ncheck)); // random order
365 
366  for (int s = 0; s < nvar + ncheck; s++) {
367  l = Ce_ind(s);
368  if (rowjk(l) != 0) { continue; }
369  ivec colcandi = G.get_col(l).get_nz_indices();
370  if (length(colcandi) > 0) {
371  // select a node p which is a member of Ce
372  for (int u = 0; u < length(colcandi); u++) {
373  p = colcandi(u);
374  if (p != l) {
375  if (rowjk(p) == 0) {
376  goto found_candidate_vector;
377  }
378  }
379  }
380  }
381  }
382  continue; // go to the next entry (i,j)
383 
384  found_candidate_vector:
385  // swap edges
386 
387  if (p >= ncheck) { int z = l; l = p; p = z; }
388  if (j >= ncheck) { int z = k; k = j; j = z; }
389 
390  // Swap endpoints of edges (p,l) and (j,k)
391  // cout << "(" << j << "," << k << ")<->("
392  // << p << "," << l << ") " ;
393  // cout << ".";
394  // cout.flush();
395 
396  // Update the matrix
397  it_assert_debug((get(j, k - ncheck) == 1) && (get(p, l - ncheck) == 1),
398  "LDPC_Parity_Unstructured::cycle_removal_MGW(): "
399  "Internal error");
400  set(j, k - ncheck, 0);
401  set(p, l - ncheck, 0);
402  it_assert_debug((get(j, l - ncheck) == 0) && (get(p, k - ncheck) == 0),
403  "LDPC_Parity_Unstructured::cycle_removal_MGW(): "
404  "Internal error");
405  set(j, l - ncheck, 1);
406  set(p, k - ncheck, 1);
407 
408  // Update adjacency matrix
409  it_assert_debug(G(p, l) == 1 && G(l, p) == 1 && G(j, k) == 1
410  && G(k, j) == 1, "G");
411  it_assert_debug(G(j, l) == 0 && G(l, j) == 0 && G(p, k) == 0
412  && G(k, p) == 0, "G");
413 
414  // Delta is the update term to G
415  Ssmat Delta(ncheck + nvar, ncheck + nvar, 2);
416  Delta.set(j, k, -1);
417  Delta.set(k, j, -1);
418  Delta.set(p, l, -1);
419  Delta.set(l, p, -1);
420  Delta.set(j, l, 1);
421  Delta.set(l, j, 1);
422  Delta.set(p, k, 1);
423  Delta.set(k, p, 1);
424 
425  // update G and its powers
426  G = G + Delta;
427  it_assert_debug(G(p, l) == 0 && G(l, p) == 0 && G(j, k) == 0
428  && G(k, j) == 0, "G");
429  it_assert_debug(G(j, l) == 1 && G(l, j) == 1 && G(p, k) == 1
430  && G(k, p) == 1, "G");
431 
432  Gpow(1) = G;
433  Gpow(2) = G * G;
434  for (int z = 3; z <= r / 2; z++) {
435  Gpow(z) = Gpow(z - 1) * G;
436  }
437 
438  traverse_again = true;
439  } // if G()...
440  } // loop over i
441  } // loop over j
442  if ((!traverse_again) && (cycles_found > 0)) { // no point continue
443  scl = r - 2;
444  goto finished;
445  }
446  }
447  while (cycles_found != 0);
448  scl = r; // there were no cycles of length r; move on to next r
449  it_info_debug("Achieved girth " << (scl + 2)
450  << ". Proceeding to next level.");
451  } // loop over r
452 
453 finished:
454  int girth = scl + 2; // scl=length of smallest cycle
455  it_info_debug("Cycle removal (MGW algoritm) finished. Graph girth: "
456  << girth << ". Cycles remaining on next girth level: "
457  << cycles_found);
458 
459  return girth;
460 }
461 
463  const ivec& R,
464  const ivec& cycopt)
465 {
466  // Method based on random permutation. Attempts to avoid placing new
467  // edges so that cycles are created. More aggressive cycle avoidance
468  // for degree-2 nodes. EGL January 2007.
469 
470  initialize(sum(R), sum(C));
471  // C, R: Target number of columns/rows with certain number of ones
472 
473  // compute number of edges
474  int Ne = 0;
475  for (int i = 0;i < C.length();i++) {
476  for (int j = 0; j < C(i); j++) {
477  for (int m = 0; m < i; m++) Ne++;
478  }
479  }
480 
481  // compute connectivity matrix
482  ivec vcon(Ne);
483  ivec ccon(Ne);
484  ivec vd(nvar);
485  ivec cd(ncheck);
486  int k = 0;
487  int l = 0;
488  for (int i = 0;i < C.length();i++) {
489  for (int j = 0; j < C(i); j++) {
490  for (int m = 0; m < i; m++) {
491  vcon(k) = l;
492  vd(l) = i;
493  k++;
494  }
495  l++;
496  }
497  }
498  k = 0;
499  l = 0;
500  for (int i = 0;i < R.length();i++) {
501  for (int j = 0; j < R(i); j++) {
502  for (int m = 0; m < i; m++) {
503  ccon(k) = l;
504  cd(l) = i;
505  k++;
506  }
507  l++;
508  }
509  }
510  it_assert(k == Ne, "C/R mismatch");
511 
512  // compute random permutations
513  ivec ind = sort_index(randu(Ne));
514  ivec cp = sort_index(randu(nvar));
515  ivec rp = sort_index(randu(ncheck));
516 
517  // set the girth goal for various variable node degrees
518  ivec Laim = zeros_i(Nmax);
519  for (int i = 0; i < length(cycopt); i++) {
520  Laim(i + 2) = cycopt(i);
521  }
522  for (int i = length(cycopt); i < Nmax - 2; i++) {
523  Laim(i + 2) = cycopt(length(cycopt) - 1);
524  }
525  it_info_debug("Running with Laim=" << Laim.left(25));
526 
527  int failures = 0;
528  const int Max_attempts = 100;
529  const int apcl = 10; // attempts before reducing girth target
530  for (int k = 0; k < Ne; k++) {
531  const int el = Ne - k - 2;
532  if (k % 250 == 0) {
533  it_info_debug("Processing edge: " << k << " out of " << Ne
534  << ". Variable node degree: " << vd(vcon(k))
535  << ". Girth target: " << Laim(vd(vcon(k)))
536  << ". Accumulated failures: " << failures);
537  }
538  const int c = cp(vcon(k));
539  int L = Laim(vd(vcon(k)));
540  int attempt = 0;
541  while (true) {
542  if (attempt > 0 && attempt % apcl == 0 && L >= 6) { L -= 2; };
543  int r = rp(ccon(ind(k)));
544  if (get(r, c)) { // double edge
545  // set(r,c,0);
546  if (el > 0) {
547  // int t=k+1+rand()%el;
548  int t = k + 1 + randi(0, el - 1);
549  int x = ind(t);
550  ind(t) = ind(k);
551  ind(k) = x;
552  attempt++;
553  if (attempt == Max_attempts) {
554  failures++;
555  break;
556  }
557  }
558  else { // almost at the last edge
559  break;
560  }
561  }
562  else {
563  set(r, c, 1);
564  if (L > 0) { // attempt to avoid cycles
565  if (check_connectivity(r, c, r, c, 0, L) >= 0) { // detected cycle
566  set(r, c, 0);
567  if (el > 0) {
568  // make a swap in the index permutation
569  // int t=k+1+rand()%el;
570  int t = k + 1 + randi(0, el - 1);
571  int x = ind(t);
572  ind(t) = ind(k);
573  ind(k) = x;
574  attempt++;
575  if (attempt == Max_attempts) { // give up
576  failures++;
577  set(r, c, 1);
578  break;
579  }
580  }
581  else { // no edges left
582  set(r, c, 1);
583  break;
584  }
585  }
586  else {
587  break;
588  }
589  }
590  else {
591  break;
592  }
593  }
594  }
595  }
596 }
597 
598 void LDPC_Parity_Unstructured::compute_CR(const vec& var_deg, const vec& chk_deg, const int Nvar,
599  ivec &C, ivec &R)
600 {
601  // compute the degree distributions from a node perspective
602  vec Vi = linspace(1, length(var_deg), length(var_deg));
603  vec Ci = linspace(1, length(chk_deg), length(chk_deg));
604  // Compute number of cols with n 1's
605  // C, R: Target number of columns/rows with certain number of ones
606  C = to_ivec(round(Nvar * elem_div(var_deg, Vi)
607  / sum(elem_div(var_deg, Vi))));
608  C = concat(0, C);
609  int edges = sum(elem_mult(to_ivec(linspace(0, C.length() - 1,
610  C.length())), C));
611  R = to_ivec(round(edges * elem_div(chk_deg, Ci)));
612  R = concat(0, R);
613  vec Ri = linspace(0, length(R) - 1, length(R));
614  vec Coli = linspace(0, length(C) - 1, length(C));
615 
616  // trim to deal with inconsistencies due to rounding errors
617  if (sum(C) != Nvar) {
618  ivec ind = find(C == max(C));
619  C(ind(0)) = C(ind(0)) - (sum(C) - Nvar);
620  }
621 
622  //the number of edges calculated from R must match the number of
623  //edges calculated from C
624  while (sum(elem_mult(to_vec(R), Ri)) !=
625  sum(elem_mult(to_vec(C), Coli))) {
626  //we're only changing R, this is probably(?) better for irac codes
627  if (sum(elem_mult(to_vec(R), Ri)) > sum(elem_mult(to_vec(C), Coli))) {
628  //remove an edge from R
629  ivec ind = find(R == max(R));
630  int old = R(ind(0));
631  R.set(ind(0), old - 1);
632  old = R(ind(0) - 1);
633  R.set(ind(0) - 1, old + 1);
634  }
635  else {
636  ivec ind = find(R == max(R));
637  if (ind(0) == R.length() - 1) {
638  R = concat(R, 0);
639  Ri = linspace(0, length(R) - 1, length(R));
640  }
641  int old = R(ind(0));
642  R.set(ind(0), old - 1);
643  old = R(ind(0) + 1);
644  R.set(ind(0) + 1, old + 1);
645  }
646  }
647 
648  C = concat(C, zeros_i(Nmax - length(C)));
649  R = concat(R, zeros_i(Nmax - length(R)));
650 
651  it_info_debug("C=" << C << std::endl);
652  it_info_debug("R=" << R << std::endl);
653 
654 }
655 
656 // ----------------------------------------------------------------------
657 // LDPC_Parity_Regular
658 // ----------------------------------------------------------------------
659 
661  const std::string& method,
662  const ivec& options)
663 {
664  generate(Nvar, k, l, method, options);
665 }
666 
667 void LDPC_Parity_Regular::generate(int Nvar, int k, int l,
668  const std::string& method,
669  const ivec& options)
670 {
671  vec var_deg = zeros(k);
672  vec chk_deg = zeros(l);
673  var_deg(k - 1) = 1;
674  chk_deg(l - 1) = 1;
675 
676  ivec C, R;
677  compute_CR(var_deg, chk_deg, Nvar, C, R);
678  it_info_debug("sum(C)=" << sum(C) << " Nvar=" << Nvar);
679  it_info_debug("sum(R)=" << sum(R) << " approximate target=" << round_i(Nvar * k / static_cast<double>(l)));
680 
681  if (method == "rand") {
682  generate_random_H(C, R, options);
683  }
684  else {
685  it_error("not implemented");
686  };
687 }
688 
689 
690 // ----------------------------------------------------------------------
691 // LDPC_Parity_Irregular
692 // ----------------------------------------------------------------------
693 
695  const vec& var_deg,
696  const vec& chk_deg,
697  const std::string& method,
698  const ivec& options)
699 {
700  generate(Nvar, var_deg, chk_deg, method, options);
701 }
702 
703 void LDPC_Parity_Irregular::generate(int Nvar, const vec& var_deg,
704  const vec& chk_deg,
705  const std::string& method,
706  const ivec& options)
707 {
708  ivec C, R;
709  compute_CR(var_deg, chk_deg, Nvar, C, R);
710 
711  if (method == "rand") {
712  generate_random_H(C, R, options);
713  }
714  else {
715  it_error("not implemented");
716  };
717 }
718 
719 // ----------------------------------------------------------------------
720 // BLDPC_Parity
721 // ----------------------------------------------------------------------
722 
723 BLDPC_Parity::BLDPC_Parity(const imat& base_matrix, int exp_factor)
724 {
725  expand_base(base_matrix, exp_factor);
726 }
727 
728 BLDPC_Parity::BLDPC_Parity(const std::string& filename, int exp_factor)
729 {
730  load_base_matrix(filename);
731  expand_base(H_b, exp_factor);
732 }
733 
734 void BLDPC_Parity::expand_base(const imat& base_matrix, int exp_factor)
735 {
736  Z = exp_factor;
737  H_b = base_matrix;
738  H_b_valid = true;
739 
740  initialize(H_b.rows() * Z, H_b.cols() * Z);
741  for (int r = 0; r < H_b.rows(); r++) {
742  for (int c = 0; c < H_b.cols(); c++) {
743  int rz = r * Z;
744  int cz = c * Z;
745  switch (H_b(r, c)) {
746  case -1:
747  break;
748  case 0:
749  for (int i = 0; i < Z; ++i)
750  set(rz + i, cz + i, 1);
751  break;
752  default:
753  for (int i = 0; i < Z; ++i)
754  set(rz + i, cz + (i + H_b(r, c)) % Z, 1);
755  break;
756  }
757  }
758  }
759 }
760 
761 
763 {
764  return Z;
765 }
766 
768 {
769  return H_b;
770 }
771 
772 void BLDPC_Parity::set_exp_factor(int exp_factor)
773 {
774  Z = exp_factor;
775  if (H_b_valid) {
776  expand_base(H_b, exp_factor);
777  }
778  else {
779  calculate_base_matrix();
780  }
781 }
782 
783 
784 void BLDPC_Parity::load_base_matrix(const std::string& filename)
785 {
786  std::ifstream bm_file(filename.c_str());
787  it_assert(bm_file.is_open(), "BLDPC_Parity::load_base_matrix(): Could not "
788  "open file \"" << filename << "\" for reading");
789  // delete old base matrix content
790  H_b.set_size(0, 0);
791  // parse text file content, row by row
792  std::string line;
793  int line_counter = 0;
794  getline(bm_file, line);
795  while (!bm_file.eof()) {
796  line_counter++;
797  std::stringstream ss(line);
798  ivec row(0);
799  while (ss.good()) {
800  int val;
801  ss >> val;
802  row = concat(row, val);
803  }
804  if ((H_b.rows() == 0) || (row.size() == H_b.cols()))
805  H_b.append_row(row);
806  else
807  it_warning("BLDPC_Parity::load_base_matrix(): Wrong size of "
808  "a parsed row number " << line_counter);
809  getline(bm_file, line);
810  }
811  bm_file.close();
812 
813  // transpose parsed base matrix if necessary
814  if (H_b.rows() > H_b.cols())
815  H_b = H_b.transpose();
816 
817  H_b_valid = true;
818  init_flag = false;
819 }
820 
821 void BLDPC_Parity::save_base_matrix(const std::string& filename) const
822 {
823  it_assert(H_b_valid, "BLDPC_Parity::save_base_matrix(): Base matrix is "
824  "not valid");
825  std::ofstream bm_file(filename.c_str());
826  it_assert(bm_file.is_open(), "BLDPC_Parity::save_base_matrix(): Could not "
827  "open file \"" << filename << "\" for writing");
828 
829  for (int r = 0; r < H_b.rows(); r++) {
830  for (int c = 0; c < H_b.cols(); c++) {
831  bm_file << std::setw(3) << H_b(r, c);
832  }
833  bm_file << "\n";
834  }
835 
836  bm_file.close();
837 }
838 
839 
840 void BLDPC_Parity::calculate_base_matrix()
841 {
842  std::string error_str = "BLDPC_Parity::calculate_base_matrix(): "
843  "Invalid BLDPC matrix. Cannot calculate base matrix from it.";
844  int rows = H.rows() / Z;
845  int cols = H.cols() / Z;
846  it_assert((rows * Z == H.rows()) && (cols * Z == H.cols()), error_str);
847  H_b.set_size(rows, cols);
848 
849  for (int r = 0; r < rows; ++r) {
850  int rz = r * Z;
851  for (int c = 0; c < cols; ++c) {
852  int cz = c * Z;
853  GF2mat_sparse H_Z = H.get_submatrix(rz, rz + Z - 1, cz, cz + Z - 1);
854 
855  if (H_Z.nnz() == 0) {
856  H_b(r, c) = -1;
857  }
858  else if (H_Z.nnz() == Z) {
859  // check for cyclic-shifted ZxZ matrix
860  int shift = 0;
861  while ((shift < Z) && (H_Z(0, shift) != 1))
862  ++shift;
863  it_assert(shift < Z, error_str);
864  for (int i = 1; i < Z; ++i)
865  it_assert(H_Z(0, shift) == H_Z(i, (i + shift) % Z), error_str);
866  H_b(r, c) = shift;
867  }
868  else {
869  it_error(error_str);
870  }
871  } // for (int c = 0; c < cols; ++c)
872  } // for (int r = 0; r < rows; ++r)
873 
874  it_info("Base matrix calculated");
875  H_b_valid = true;
876 }
877 
878 
879 
880 // ----------------------------------------------------------------------
881 // LDPC_Generator_Systematic
882 // ----------------------------------------------------------------------
883 
885  bool natural_ordering,
886  const ivec& ind):
887  LDPC_Generator("systematic"), G()
888 {
889  ivec tmp;
890  tmp = construct(H, natural_ordering, ind);
891 }
892 
893 
895  bool natural_ordering,
896  const ivec& avoid_cols)
897 {
898  int nvar = H->get_nvar();
899  int ncheck = H->get_ncheck();
900 
901  // create dense representation of parity check matrix
902  GF2mat Hd(H->get_H());
903 
904  // -- Determine initial column ordering --
905  ivec col_order(nvar);
906  if (natural_ordering) {
907  for (int i = 0; i < nvar; i++) {
908  col_order(i) = i;
909  }
910  }
911  else {
912  // take the columns in random order, but the ones to avoid at last
913  vec col_importance = randu(nvar);
914  for (int i = 0; i < length(avoid_cols); i++) {
915  col_importance(avoid_cols(i)) = (-col_importance(avoid_cols(i)));
916  }
917  col_order = sort_index(-col_importance);
918  }
919 
920  ivec actual_ordering(nvar);
921 
922  // Now partition P as P=[P1 P2]. Then find G so [P1 P2][I G]'=0.
923 
924  // -- Create P1 and P2 --
925  GF2mat P1; //(ncheck,nvar-ncheck); // non-invertible part
926  GF2mat P2; //(ncheck,ncheck); // invertible part
927 
928  it_info_debug("Computing a systematic generator matrix...");
929 
930  int j1 = 0, j2 = 0;
931  int rank;
932  ivec perm;
933  GF2mat T, U;
934  for (int k = 0; k < nvar; k++) {
935  it_error_if(j1 >= nvar - ncheck, "LDPC_Generator_Systematic::construct(): "
936  "Unable to obtain enough independent columns.");
937 
938  bvec c = Hd.get_col(col_order(k));
939  if (j2 == 0) { // first column in P2 is number col_order(k)
940  P2 = GF2mat(c);
941  rank = P2.T_fact(T, U, perm);
942  actual_ordering(k) = nvar - ncheck;
943  j2++;
944  }
945  else {
946  if (j2 < ncheck) {
947  if (P2.T_fact_update_addcol(T, U, perm, c)) {
948  P2 = P2.concatenate_horizontal(c);
949  actual_ordering(k) = nvar - ncheck + j2;
950  j2++;
951  continue;
952  }
953  }
954  if (j1 == 0) {
955  P1 = GF2mat(c);
956  actual_ordering(k) = j1;
957  }
958  else {
959  P1 = P1.concatenate_horizontal(c);
960  actual_ordering(k) = j1;
961  }
962  j1++;
963  }
964  }
965 
966  it_info_debug("Rank of parity check matrix: " << j2);
967 
968  // -- Compute the systematic part of the generator matrix --
969  G = (P2.inverse() * P1).transpose();
970 
971  // -- Permute the columns of the parity check matrix --
972  GF2mat P = P1.concatenate_horizontal(P2);
973  *H = LDPC_Parity(ncheck, nvar);
974  // brute force copy from P to H
975  for (int i = 0; i < ncheck; i++) {
976  for (int j = 0; j < nvar; j++) {
977  if (P.get(i, j)) {
978  H->set(i, j, 1);
979  }
980  }
981  }
982 
983  // -- Check that the result was correct --
985  * (gf2dense_eye(nvar - ncheck).concatenate_horizontal(G)).transpose()).is_zero(),
986  "LDPC_Generator_Systematic::construct(): Incorrect generator matrix G");
987 
988  G = G.transpose(); // store the generator matrix in transposed form
989  it_info_debug("Systematic generator matrix computed.");
990 
991  init_flag = true;
992 
993  return actual_ordering;
994 }
995 
996 
997 void LDPC_Generator_Systematic::save(const std::string& filename) const
998 {
999  it_file f(filename);
1000  int ver;
1001  f >> Name("Fileversion") >> ver;
1003  "LDPC_Generator_Systematic::save(): Unsupported file format");
1004  f << Name("G_type") << type;
1005  f << Name("G") << G;
1006  f.close();
1007 }
1008 
1009 
1010 void LDPC_Generator_Systematic::load(const std::string& filename)
1011 {
1012  it_ifile f(filename);
1013  int ver;
1014  f >> Name("Fileversion") >> ver;
1016  "LDPC_Generator_Systematic::load(): Unsupported file format");
1017  std::string gen_type;
1018  f >> Name("G_type") >> gen_type;
1019  it_assert(gen_type == type,
1020  "LDPC_Generator_Systematic::load(): Wrong generator type");
1021  f >> Name("G") >> G;
1022  f.close();
1023 
1024  init_flag = true;
1025 }
1026 
1027 
1028 void LDPC_Generator_Systematic::encode(const bvec &input, bvec &output)
1029 {
1030  it_assert(init_flag, "LDPC_Generator_Systematic::encode(): Systematic "
1031  "generator not set up");
1032  it_assert(input.size() == G.cols(), "LDPC_Generator_Systematic::encode(): "
1033  "Improper input vector size (" << input.size() << " != "
1034  << G.cols() << ")");
1035 
1036  output = concat(input, G * input);
1037 }
1038 
1039 
1040 // ----------------------------------------------------------------------
1041 // BLDPC_Generator
1042 // ----------------------------------------------------------------------
1043 
1045  const std::string type):
1046  LDPC_Generator(type), H_enc(), N(0), M(0), K(0), Z(0)
1047 {
1048  construct(H);
1049 }
1050 
1051 
1052 void BLDPC_Generator::encode(const bvec &input, bvec &output)
1053 {
1054  it_assert(init_flag, "BLDPC_Generator::encode(): Cannot encode with not "
1055  "initialized generator");
1056  it_assert(input.size() == K, "BLDPC_Generator::encode(): Input vector "
1057  "length is not equal to K");
1058 
1059  // copy systematic bits first
1060  output = input;
1061  output.set_size(N, true);
1062 
1063  // backward substitution to obtain the first Z parity check bits
1064  for (int k = 0; k < Z; k++) {
1065  for (int j = 0; j < K; j++) {
1066  output(K + k) += H_enc(M - 1 - k, j) * input(j);
1067  }
1068  for (int j = 0; j < k; j++) {
1069  output(K + k) += H_enc(M - 1 - k, K + j) * output(K + j);
1070  }
1071  }
1072 
1073  // forward substitution to obtain the next M-Z parity check bits
1074  for (int k = 0; k < M - Z; k++) {
1075  for (int j = 0; j < K; j++) {
1076  output(K + Z + k) += H_enc(k, j) * input(j);
1077  }
1078  for (int j = K; j < K + Z; j++) {
1079  output(K + Z + k) += H_enc(k, j) * output(j);
1080  }
1081  for (int j = K + Z; j < K + Z + k; j++) {
1082  output(K + Z + k) += H_enc(k, j) * output(j);
1083  }
1084  }
1085 }
1086 
1087 
1089 {
1090  if (H != 0 && H->is_valid()) {
1091  H_enc = H->get_H(); // sparse to dense conversion
1092  Z = H->get_exp_factor();
1093  N = H_enc.cols();
1094  M = H_enc.rows();
1095  K = N - M;
1096 
1097  // ----------------------------------------------------------------------
1098  // STEP 1
1099  // ----------------------------------------------------------------------
1100  // loop over last M-Z columns of matrix H
1101  for (int i = 0; i < M - Z; i += Z) {
1102  // loop over last Z rows of matrix H
1103  for (int j = 0; j < Z; j++) {
1104  // Gaussian elimination by adding particular rows
1105  H_enc.add_rows(M - 1 - j, M - Z - 1 - j - i);
1106  }
1107  }
1108 
1109  // ----------------------------------------------------------------------
1110  // STEP 2
1111  // ----------------------------------------------------------------------
1112  // set first processed row index to M-Z
1113  int r1 = M - Z;
1114  // loop over columns with indexes K .. K+Z-1
1115  for (int c1 = K + Z - 1; c1 >= K; c1--) {
1116  int r2 = r1;
1117  // find the first '1' in column c1
1118  while (H_enc(r2, c1) == 0 && r2 < M - 1)
1119  r2++;
1120  // if necessary, swap rows r1 and r2
1121  if (r2 != r1)
1122  H_enc.swap_rows(r1, r2);
1123 
1124  // look for the other '1' in column c1 and get rid of them
1125  for (r2 = r1 + 1; r2 < M; r2++) {
1126  if (H_enc(r2, c1) == 1) {
1127  // Gaussian elimination by adding particular rows
1128  H_enc.add_rows(r2, r1);
1129  }
1130  }
1131 
1132  // increase first processed row index
1133  r1++;
1134  }
1135 
1136  init_flag = true;
1137  }
1138 }
1139 
1140 
1141 void BLDPC_Generator::save(const std::string& filename) const
1142 {
1144  "BLDPC_Generator::save(): Can not save not initialized generator");
1145  // Every Z-th row of H_enc until M-Z
1146  GF2mat H_T(M / Z - 1, N);
1147  for (int i = 0; i < M / Z - 1; i++) {
1148  H_T.set_row(i, H_enc.get_row(i*Z));
1149  }
1150  // Last Z preprocessed rows of H_enc
1151  GF2mat H_Z = H_enc.get_submatrix(M - Z, 0, M - 1, N - 1);
1152 
1153  it_file f(filename);
1154  int ver;
1155  f >> Name("Fileversion") >> ver;
1156  it_assert(ver == LDPC_binary_file_version, "BLDPC_Generator::save(): "
1157  "Unsupported file format");
1158  f << Name("G_type") << type;
1159  f << Name("H_T") << H_T;
1160  f << Name("H_Z") << H_Z;
1161  f << Name("Z") << Z;
1162  f.close();
1163 }
1164 
1165 void BLDPC_Generator::load(const std::string& filename)
1166 {
1167  GF2mat H_T, H_Z;
1168 
1169  it_ifile f(filename);
1170  int ver;
1171  f >> Name("Fileversion") >> ver;
1172  it_assert(ver == LDPC_binary_file_version, "BLDPC_Generator::load(): "
1173  "Unsupported file format");
1174  std::string gen_type;
1175  f >> Name("G_type") >> gen_type;
1176  it_assert(gen_type == type,
1177  "BLDPC_Generator::load(): Wrong generator type");
1178  f >> Name("H_T") >> H_T;
1179  f >> Name("H_Z") >> H_Z;
1180  f >> Name("Z") >> Z;
1181  f.close();
1182 
1183  N = H_T.cols();
1184  M = (H_T.rows() + 1) * Z;
1185  K = N - M;
1186 
1187  H_enc = GF2mat(M - Z, N);
1188  for (int i = 0; i < H_T.rows(); i++) {
1189  for (int j = 0; j < Z; j++) {
1190  for (int k = 0; k < N; k++) {
1191  if (H_T(i, (k / Z)*Z + (k + Z - j) % Z)) {
1192  H_enc.set(i*Z + j, k, 1);
1193  }
1194  }
1195  }
1196  }
1198 
1199  init_flag = true;
1200 }
1201 
1202 
1203 // ----------------------------------------------------------------------
1204 // LDPC_Code
1205 // ----------------------------------------------------------------------
1206 
1207 LDPC_Code::LDPC_Code(): H_defined(false), G_defined(false), dec_method("BP"),
1208  max_iters(50), psc(true), pisc(false),
1209  llrcalc(LLR_calc_unit()) { }
1210 
1212  LDPC_Generator* const G_in):
1213  H_defined(false), G_defined(false), dec_method("BP"), max_iters(50),
1214  psc(true), pisc(false), llrcalc(LLR_calc_unit())
1215 {
1216  set_code(H, G_in);
1217 }
1218 
1219 LDPC_Code::LDPC_Code(const std::string& filename,
1220  LDPC_Generator* const G_in):
1221  H_defined(false), G_defined(false), dec_method("BP"), max_iters(50),
1222  psc(true), pisc(false), llrcalc(LLR_calc_unit())
1223 {
1224  load_code(filename, G_in);
1225 }
1226 
1227 
1228 void LDPC_Code::set_code(const LDPC_Parity* const H,
1229  LDPC_Generator* const G_in)
1230 {
1232  setup_decoder();
1233  G = G_in;
1234  if (G != 0) {
1235  G_defined = true;
1236  integrity_check();
1237  }
1238 }
1239 
1240 void LDPC_Code::load_code(const std::string& filename,
1241  LDPC_Generator* const G_in)
1242 {
1243  it_info_debug("LDPC_Code::load_code(): Loading LDPC codec from "
1244  << filename);
1245 
1246  it_ifile f(filename);
1247  int ver;
1248  f >> Name("Fileversion") >> ver;
1249  it_assert(ver == LDPC_binary_file_version, "LDPC_Code::load_code(): "
1250  "Unsupported file format");
1251  f >> Name("H_defined") >> H_defined;
1252  f >> Name("G_defined") >> G_defined;
1253  f >> Name("nvar") >> nvar;
1254  f >> Name("ncheck") >> ncheck;
1255  f >> Name("C") >> C;
1256  f >> Name("V") >> V;
1257  f >> Name("sumX1") >> sumX1;
1258  f >> Name("sumX2") >> sumX2;
1259  f >> Name("iind") >> iind;
1260  f >> Name("jind") >> jind;
1261  f.close();
1262 
1263  // load generator data
1264  if (G_defined) {
1265  it_assert(G_in != 0, "LDPC_Code::load_code(): Generator object is "
1266  "missing. Can not load the generator data from a file.");
1267  G = G_in;
1268  G->load(filename);
1269  }
1270  else {
1271  G = 0;
1272  it_info_debug("LDPC_Code::load_code(): Generator data not loaded. "
1273  "Generator object will not be used.");
1274  }
1275 
1276  it_info_debug("LDPC_Code::load_code(): Successfully loaded LDPC codec "
1277  "from " << filename);
1278 
1279  setup_decoder();
1280 }
1281 
1282 void LDPC_Code::save_code(const std::string& filename) const
1283 {
1284  it_assert(H_defined, "LDPC_Code::save_to_file(): There is no parity "
1285  "check matrix");
1286  it_info_debug("LDPC_Code::save_to_file(): Saving LDPC codec to "
1287  << filename);
1288 
1289  it_file f;
1290  f.open(filename, true);
1291  f << Name("Fileversion") << LDPC_binary_file_version;
1292  f << Name("H_defined") << H_defined;
1293  f << Name("G_defined") << G_defined;
1294  f << Name("nvar") << nvar;
1295  f << Name("ncheck") << ncheck;
1296  f << Name("C") << C;
1297  f << Name("V") << V;
1298  f << Name("sumX1") << sumX1;
1299  f << Name("sumX2") << sumX2;
1300  f << Name("iind") << iind;
1301  f << Name("jind") << jind;
1302  f.close();
1303 
1304  // save generator data;
1305  if (G_defined)
1306  G->save(filename);
1307  else
1308  it_info_debug("LDPC_Code::save_code(): Missing generator object - "
1309  "generator data not saved");
1310 
1311  it_info_debug("LDPC_Code::save_code(): Successfully saved LDPC codec to "
1312  << filename);
1313 }
1314 
1315 
1316 void LDPC_Code::set_decoding_method(const std::string& method_in)
1317 {
1318  it_assert((method_in == "bp") || (method_in == "BP"),
1319  "LDPC_Code::set_decoding_method(): Not implemented decoding method");
1320  dec_method = method_in;
1321 }
1322 
1323 void LDPC_Code::set_exit_conditions(int max_iters_in,
1324  bool syndr_check_each_iter,
1325  bool syndr_check_at_start)
1326 {
1327  it_assert(max_iters >= 0, "LDPC_Code::set_nrof_iterations(): Maximum "
1328  "number of iterations can not be negative");
1329  max_iters = max_iters_in;
1330  psc = syndr_check_each_iter;
1331  pisc = syndr_check_at_start;
1332 }
1333 
1334 void LDPC_Code::set_llrcalc(const LLR_calc_unit& llrcalc_in)
1335 {
1336  llrcalc = llrcalc_in;
1337 }
1338 
1339 
1340 void LDPC_Code::encode(const bvec &input, bvec &output)
1341 {
1342  it_assert(G_defined, "LDPC_Code::encode(): LDPC Generator is required "
1343  "for encoding");
1344  G->encode(input, output);
1345  it_assert_debug(syndrome_check(output), "LDPC_Code::encode(): Syndrome "
1346  "check failed");
1347 }
1348 
1349 bvec LDPC_Code::encode(const bvec &input)
1350 {
1351  bvec result;
1352  encode(input, result);
1353  return result;
1354 }
1355 
1356 void LDPC_Code::decode(const vec &llr_in, bvec &syst_bits)
1357 {
1358  QLLRvec qllrin = llrcalc.to_qllr(llr_in);
1359  QLLRvec qllrout;
1360  bp_decode(qllrin, qllrout);
1361  syst_bits = (qllrout.left(nvar - ncheck) < 0);
1362 }
1363 
1364 bvec LDPC_Code::decode(const vec &llr_in)
1365 {
1366  bvec syst_bits;
1367  decode(llr_in, syst_bits);
1368  return syst_bits;
1369 }
1370 
1371 void LDPC_Code::decode_soft_out(const vec &llr_in, vec &llr_out)
1372 {
1373  QLLRvec qllrin = llrcalc.to_qllr(llr_in);
1374  QLLRvec qllrout;
1375  bp_decode(qllrin, qllrout);
1376  llr_out = llrcalc.to_double(qllrout);
1377 }
1378 
1379 vec LDPC_Code::decode_soft_out(const vec &llr_in)
1380 {
1381  vec llr_out;
1382  decode_soft_out(llr_in, llr_out);
1383  return llr_out;
1384 }
1385 
1386 int LDPC_Code::bp_decode(const QLLRvec &LLRin, QLLRvec &LLRout)
1387 {
1388  // Note the IT++ convention that a sure zero corresponds to
1389  // LLR=+infinity
1390  it_assert(H_defined, "LDPC_Code::bp_decode(): Parity check matrix not "
1391  "defined");
1392  it_assert((LLRin.size() == nvar) && (sumX1.size() == nvar)
1393  && (sumX2.size() == ncheck), "LDPC_Code::bp_decode(): Wrong "
1394  "input dimensions");
1395 
1396  if (pisc && syndrome_check(LLRin)) {
1397  LLRout = LLRin;
1398  return 0;
1399  }
1400 
1401  LLRout.set_size(LLRin.size());
1402 
1403  // allocate temporary variables used for the check node update
1404  ivec jj(max_cnd);
1405  QLLRvec m(max_cnd);
1406  QLLRvec ml(max_cnd);
1407  QLLRvec mr(max_cnd);
1408 
1409  // initial step
1410  for (int i = 0; i < nvar; i++) {
1411  int index = i;
1412  for (int j = 0; j < sumX1(i); j++) {
1413  mvc[index] = LLRin(i);
1414  index += nvar;
1415  }
1416  }
1417 
1418  bool is_valid_codeword = false;
1419  int iter = 0;
1420  do {
1421  iter++;
1422  if (nvar >= 100000) { it_info_no_endl_debug("."); }
1423  // --------- Step 1: check to variable nodes ----------
1424  for (int j = 0; j < ncheck; j++) {
1425  // The check node update calculations are hardcoded for degrees
1426  // up to 6. For larger degrees, a general algorithm is used.
1427  switch (sumX2(j)) {
1428  case 0:
1429  it_error("LDPC_Code::bp_decode(): sumX2(j)=0");
1430  case 1:
1431  it_error("LDPC_Code::bp_decode(): sumX2(j)=1");
1432  case 2: {
1433  mcv[j+ncheck] = mvc[jind[j]];
1434  mcv[j] = mvc[jind[j+ncheck]];
1435  break;
1436  }
1437  case 3: {
1438  int j0 = j;
1439  QLLR m0 = mvc[jind[j0]];
1440  int j1 = j0 + ncheck;
1441  QLLR m1 = mvc[jind[j1]];
1442  int j2 = j1 + ncheck;
1443  QLLR m2 = mvc[jind[j2]];
1444  mcv[j0] = llrcalc.Boxplus(m1, m2);
1445  mcv[j1] = llrcalc.Boxplus(m0, m2);
1446  mcv[j2] = llrcalc.Boxplus(m0, m1);
1447  break;
1448  }
1449  case 4: {
1450  int j0 = j;
1451  QLLR m0 = mvc[jind[j0]];
1452  int j1 = j0 + ncheck;
1453  QLLR m1 = mvc[jind[j1]];
1454  int j2 = j1 + ncheck;
1455  QLLR m2 = mvc[jind[j2]];
1456  int j3 = j2 + ncheck;
1457  QLLR m3 = mvc[jind[j3]];
1458  QLLR m01 = llrcalc.Boxplus(m0, m1);
1459  QLLR m23 = llrcalc.Boxplus(m2, m3);
1460  mcv[j0] = llrcalc.Boxplus(m1, m23);
1461  mcv[j1] = llrcalc.Boxplus(m0, m23);
1462  mcv[j2] = llrcalc.Boxplus(m01, m3);
1463  mcv[j3] = llrcalc.Boxplus(m01, m2);
1464  break;
1465  }
1466  case 5: {
1467  int j0 = j;
1468  QLLR m0 = mvc[jind[j0]];
1469  int j1 = j0 + ncheck;
1470  QLLR m1 = mvc[jind[j1]];
1471  int j2 = j1 + ncheck;
1472  QLLR m2 = mvc[jind[j2]];
1473  int j3 = j2 + ncheck;
1474  QLLR m3 = mvc[jind[j3]];
1475  int j4 = j3 + ncheck;
1476  QLLR m4 = mvc[jind[j4]];
1477  QLLR m01 = llrcalc.Boxplus(m0, m1);
1478  QLLR m02 = llrcalc.Boxplus(m01, m2);
1479  QLLR m34 = llrcalc.Boxplus(m3, m4);
1480  QLLR m24 = llrcalc.Boxplus(m2, m34);
1481  mcv[j0] = llrcalc.Boxplus(m1, m24);
1482  mcv[j1] = llrcalc.Boxplus(m0, m24);
1483  mcv[j2] = llrcalc.Boxplus(m01, m34);
1484  mcv[j3] = llrcalc.Boxplus(m02, m4);
1485  mcv[j4] = llrcalc.Boxplus(m02, m3);
1486  break;
1487  }
1488  case 6: {
1489  int j0 = j;
1490  QLLR m0 = mvc[jind[j0]];
1491  int j1 = j0 + ncheck;
1492  QLLR m1 = mvc[jind[j1]];
1493  int j2 = j1 + ncheck;
1494  QLLR m2 = mvc[jind[j2]];
1495  int j3 = j2 + ncheck;
1496  QLLR m3 = mvc[jind[j3]];
1497  int j4 = j3 + ncheck;
1498  QLLR m4 = mvc[jind[j4]];
1499  int j5 = j4 + ncheck;
1500  QLLR m5 = mvc[jind[j5]];
1501  QLLR m01 = llrcalc.Boxplus(m0, m1);
1502  QLLR m23 = llrcalc.Boxplus(m2, m3);
1503  QLLR m45 = llrcalc.Boxplus(m4, m5);
1504  QLLR m03 = llrcalc.Boxplus(m01, m23);
1505  QLLR m25 = llrcalc.Boxplus(m23, m45);
1506  QLLR m0145 = llrcalc.Boxplus(m01, m45);
1507  mcv[j0] = llrcalc.Boxplus(m1, m25);
1508  mcv[j1] = llrcalc.Boxplus(m0, m25);
1509  mcv[j2] = llrcalc.Boxplus(m0145, m3);
1510  mcv[j3] = llrcalc.Boxplus(m0145, m2);
1511  mcv[j4] = llrcalc.Boxplus(m03, m5);
1512  mcv[j5] = llrcalc.Boxplus(m03, m4);
1513  break;
1514  }
1515  default: {
1516  int nodes = sumX2(j);
1517  if( nodes > max_cnd ) {
1518  std::ostringstream m_sout;
1519  m_sout << "check node degrees >" << max_cnd << " not supported in this version";
1520  it_error( m_sout.str() );
1521  }
1522 
1523  nodes--;
1524  jj[0] = j;
1525  m[0] = mvc[jind[jj[0]]];
1526  for(int i = 1; i <= nodes; i++ ) {
1527  jj[i] = jj[i-1] + ncheck;
1528  m[i] = mvc[jind[jj[i]]];
1529  }
1530 
1531  // compute partial sums from the left and from the right
1532  ml[0] = m[0];
1533  mr[0] = m[nodes];
1534  for(int i = 1; i < nodes; i++ ) {
1535  ml[i] = llrcalc.Boxplus( ml[i-1], m[i] );
1536  mr[i] = llrcalc.Boxplus( mr[i-1], m[nodes-i] );
1537  }
1538 
1539  // merge partial sums
1540  mcv[jj[0]] = mr[nodes-1];
1541  mcv[jj[nodes]] = ml[nodes-1];
1542  for(int i = 1; i < nodes; i++ )
1543  mcv[jj[i]] = llrcalc.Boxplus( ml[i-1], mr[nodes-1-i] );
1544  }
1545  } // switch statement
1546  }
1547 
1548  // step 2: variable to check nodes
1549  for (int i = 0; i < nvar; i++) {
1550  switch (sumX1(i)) {
1551  case 0:
1552  it_error("LDPC_Code::bp_decode(): sumX1(i)=0");
1553  case 1: {
1554  // Degenerate case-should not occur for good coded. A lonely
1555  // variable node only sends its incoming message
1556  mvc[i] = LLRin(i);
1557  LLRout(i) = LLRin(i);
1558  break;
1559  }
1560  case 2: {
1561  QLLR m0 = mcv[iind[i]];
1562  int i1 = i + nvar;
1563  QLLR m1 = mcv[iind[i1]];
1564  mvc[i] = LLRin(i) + m1;
1565  mvc[i1] = LLRin(i) + m0;
1566  LLRout(i) = mvc[i1] + m1;
1567  break;
1568  }
1569  case 3: {
1570  int i0 = i;
1571  QLLR m0 = mcv[iind[i0]];
1572  int i1 = i0 + nvar;
1573  QLLR m1 = mcv[iind[i1]];
1574  int i2 = i1 + nvar;
1575  QLLR m2 = mcv[iind[i2]];
1576  LLRout(i) = LLRin(i) + m0 + m1 + m2;
1577  mvc[i0] = LLRout(i) - m0;
1578  mvc[i1] = LLRout(i) - m1;
1579  mvc[i2] = LLRout(i) - m2;
1580  break;
1581  }
1582  case 4: {
1583  int i0 = i;
1584  QLLR m0 = mcv[iind[i0]];
1585  int i1 = i0 + nvar;
1586  QLLR m1 = mcv[iind[i1]];
1587  int i2 = i1 + nvar;
1588  QLLR m2 = mcv[iind[i2]];
1589  int i3 = i2 + nvar;
1590  QLLR m3 = mcv[iind[i3]];
1591  LLRout(i) = LLRin(i) + m0 + m1 + m2 + m3;
1592  mvc[i0] = LLRout(i) - m0;
1593  mvc[i1] = LLRout(i) - m1;
1594  mvc[i2] = LLRout(i) - m2;
1595  mvc[i3] = LLRout(i) - m3;
1596  break;
1597  }
1598  default: { // differential update
1599  QLLR mvc_temp = LLRin(i);
1600  int index_iind = i; // tracks i+jp*nvar
1601  for (int jp = 0; jp < sumX1(i); jp++) {
1602  mvc_temp += mcv[iind[index_iind]];
1603  index_iind += nvar;
1604  }
1605  LLRout(i) = mvc_temp;
1606  index_iind = i; // tracks i+j*nvar
1607  for (int j = 0; j < sumX1[i]; j++) {
1608  mvc[index_iind] = mvc_temp - mcv[iind[index_iind]];
1609  index_iind += nvar;
1610  }
1611  }
1612  }
1613  }
1614 
1615  if (psc && syndrome_check(LLRout)) {
1616  is_valid_codeword = true;
1617  break;
1618  }
1619  }
1620  while (iter < max_iters);
1621 
1622  if (nvar >= 100000) { it_info_debug(""); }
1623  return (is_valid_codeword ? iter : -iter);
1624 }
1625 
1626 
1627 bool LDPC_Code::syndrome_check(const bvec &x) const
1628 {
1629  QLLRvec llr = 1 - 2 * to_ivec(x);
1630  return syndrome_check(llr);
1631 }
1632 
1633 bool LDPC_Code::syndrome_check(const QLLRvec &LLR) const
1634 {
1635  // Please note the IT++ convention that a sure zero corresponds to
1636  // LLR=+infinity
1637  int i, j, synd, vi;
1638 
1639  for (j = 0; j < ncheck; j++) {
1640  synd = 0;
1641  int vind = j; // tracks j+i*ncheck
1642  for (i = 0; i < sumX2(j); i++) {
1643  vi = V(vind);
1644  if (LLR(vi) < 0) {
1645  synd++;
1646  }
1647  vind += ncheck;
1648  }
1649  if ((synd&1) == 1) {
1650  return false; // codeword is invalid
1651  }
1652  }
1653  return true; // codeword is valid
1654 };
1655 
1656 
1657 // ----------------------------------------------------------------------
1658 // LDPC_Code private methods
1659 // ----------------------------------------------------------------------
1660 
1662 {
1663  // copy basic parameters
1664  sumX1 = Hmat->sumX1;
1665  sumX2 = Hmat->sumX2;
1666  nvar = Hmat->nvar; //get_nvar();
1667  ncheck = Hmat->ncheck; //get_ncheck();
1668  int cmax = max(sumX1);
1669  int vmax = max(sumX2);
1670 
1671  // decoder parameterization
1672  V = zeros_i(ncheck * vmax);
1673  C = zeros_i(cmax * nvar);
1674  jind = zeros_i(ncheck * vmax);
1675  iind = zeros_i(nvar * cmax);
1676 
1677  it_info_debug("LDPC_Code::decoder_parameterization(): Computations "
1678  "- phase 1");
1679  for (int i = 0; i < nvar; i++) {
1680  ivec coli = Hmat->get_col(i).get_nz_indices();
1681  for (int j0 = 0; j0 < length(coli); j0++) {
1682  C(j0 + cmax*i) = coli(j0);
1683  }
1684  }
1685 
1686  it_info_debug("LDPC_Code::decoder_parameterization(): Computations "
1687  "- phase 2");
1688  it_info_debug("Computing decoder parameterization. Phase 2");
1689  for (int j = 0; j < ncheck; j++) {
1690  ivec rowj = Hmat->get_row(j).get_nz_indices();
1691  for (int i0 = 0; i0 < length(rowj); i0++) {
1692  V(j + ncheck*i0) = rowj(i0);
1693  }
1694  }
1695 
1696  it_info_debug("LDPC_Code::decoder_parameterization(): Computations "
1697  "- phase 3");
1698  it_info_debug("Computing decoder parameterization. Phase 3.");
1699  for (int j = 0; j < ncheck; j++) {
1700  for (int ip = 0; ip < sumX2(j); ip++) {
1701  int vip = V(j + ip * ncheck);
1702  int k = 0;
1703  while (1 == 1) {
1704  if (C(k + vip*cmax) == j) {
1705  break;
1706  }
1707  k++;
1708  }
1709  jind(j + ip*ncheck) = vip + k * nvar;
1710  }
1711  }
1712 
1713  it_info_debug("LDPC_Code::decoder_parameterization(): Computations "
1714  "- phase 4");
1715  for (int i = 0; i < nvar; i++) {
1716  for (int jp = 0; jp < sumX1(i); jp++) {
1717  int cjp = C(jp + i * cmax);
1718  int k = 0;
1719  while (1 == 1) {
1720  if (V(cjp + k*ncheck) == i) {break; }
1721  k++;
1722  }
1723  iind(i + jp*nvar) = cjp + k * ncheck;
1724  }
1725  }
1726 
1727  H_defined = true;
1728 }
1729 
1730 
1732 {
1733  if (H_defined) {
1734  mcv.set_size(max(sumX2) * ncheck);
1735  mvc.set_size(max(sumX1) * nvar);
1736  }
1737 }
1738 
1739 
1741 {
1742  if (G_defined) {
1743  it_info_debug("LDPC_Code::integrity_check(): Checking integrity of "
1744  "the LDPC_Parity and LDPC_Generator data");
1745  bvec bv(nvar - ncheck), cw;
1746  bv.clear();
1747  bv(0) = 1;
1748  for (int i = 0; i < nvar - ncheck; i++) {
1749  G->encode(bv, cw);
1751  "LDPC_Code::integrity_check(): Syndrome check failed");
1752  bv.shift_right(bv(nvar - ncheck - 1));
1753  }
1754  }
1755  else {
1756  it_info_debug("LDPC_Code::integrity_check(): No generator defined "
1757  "- no check performed");
1758  }
1759 }
1760 
1761 // ----------------------------------------------------------------------
1762 // Related functions
1763 // ----------------------------------------------------------------------
1764 
1765 std::ostream &operator<<(std::ostream &os, const LDPC_Code &C)
1766 {
1767  ivec rdeg = zeros_i(max(C.sumX2) + 1);
1768  for (int i = 0; i < C.ncheck; i++) {
1769  rdeg(C.sumX2(i))++;
1770  }
1771 
1772  ivec cdeg = zeros_i(max(C.sumX1) + 1);
1773  for (int j = 0; j < C.nvar; j++) {
1774  cdeg(C.sumX1(j))++;
1775  }
1776 
1777  os << "--- LDPC codec ----------------------------------\n"
1778  << "Nvar : " << C.get_nvar() << "\n"
1779  << "Ncheck : " << C.get_ncheck() << "\n"
1780  << "Rate : " << C.get_rate() << "\n"
1781  << "Column degrees (node perspective): " << cdeg << "\n"
1782  << "Row degrees (node perspective): " << rdeg << "\n"
1783  << "-------------------------------------------------\n"
1784  << "Decoder parameters:\n"
1785  << " - method : " << C.dec_method << "\n"
1786  << " - max. iterations : " << C.max_iters << "\n"
1787  << " - syndrome check at each iteration : " << C.psc << "\n"
1788  << " - syndrome check at start : " << C.pisc << "\n"
1789  << "-------------------------------------------------\n"
1790  << C.llrcalc << "\n";
1791  return os;
1792 }
1793 
1794 } // namespace itpp
SourceForge Logo

Generated on Fri Mar 21 2014 17:14:13 for IT++ by Doxygen 1.8.1.2