IT++ Logo
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  ivec perm;
932  GF2mat T, U;
933  for (int k = 0; k < nvar; k++) {
934  it_error_if(j1 >= nvar - ncheck, "LDPC_Generator_Systematic::construct(): "
935  "Unable to obtain enough independent columns.");
936 
937  bvec c = Hd.get_col(col_order(k));
938  if (j2 == 0) { // first column in P2 is number col_order(k)
939  P2 = GF2mat(c);
940  P2.T_fact(T, U, perm);/* return value not used */
941  actual_ordering(k) = nvar - ncheck;
942  j2++;
943  }
944  else {
945  if (j2 < ncheck) {
946  if (P2.T_fact_update_addcol(T, U, perm, c)) {
947  P2 = P2.concatenate_horizontal(c);
948  actual_ordering(k) = nvar - ncheck + j2;
949  j2++;
950  continue;
951  }
952  }
953  if (j1 == 0) {
954  P1 = GF2mat(c);
955  actual_ordering(k) = j1;
956  }
957  else {
958  P1 = P1.concatenate_horizontal(c);
959  actual_ordering(k) = j1;
960  }
961  j1++;
962  }
963  }
964 
965  it_info_debug("Rank of parity check matrix: " << j2);
966 
967  // -- Compute the systematic part of the generator matrix --
968  G = (P2.inverse() * P1).transpose();
969 
970  // -- Permute the columns of the parity check matrix --
971  GF2mat P = P1.concatenate_horizontal(P2);
972  *H = LDPC_Parity(ncheck, nvar);
973  // brute force copy from P to H
974  for (int i = 0; i < ncheck; i++) {
975  for (int j = 0; j < nvar; j++) {
976  if (P.get(i, j)) {
977  H->set(i, j, 1);
978  }
979  }
980  }
981 
982  // -- Check that the result was correct --
984  * (gf2dense_eye(nvar - ncheck).concatenate_horizontal(G)).transpose()).is_zero(),
985  "LDPC_Generator_Systematic::construct(): Incorrect generator matrix G");
986 
987  G = G.transpose(); // store the generator matrix in transposed form
988  it_info_debug("Systematic generator matrix computed.");
989 
991 
992  return actual_ordering;
993 }
994 
995 
996 void LDPC_Generator_Systematic::save(const std::string& filename) const
997 {
998  it_file f(filename);
999  int ver;
1000  f >> Name("Fileversion") >> ver;
1002  "LDPC_Generator_Systematic::save(): Unsupported file format");
1003  f << Name("G_type") << get_type();
1004  f << Name("G") << G;
1005  f.close();
1006 }
1007 
1008 
1009 void LDPC_Generator_Systematic::load(const std::string& filename)
1010 {
1011  it_ifile f(filename);
1012  int ver;
1013  f >> Name("Fileversion") >> ver;
1015  "LDPC_Generator_Systematic::load(): Unsupported file format");
1016  std::string gen_type;
1017  f >> Name("G_type") >> gen_type;
1018  it_assert(gen_type == get_type(),
1019  "LDPC_Generator_Systematic::load(): Wrong generator type");
1020  f >> Name("G") >> G;
1021  f.close();
1022 
1023  mark_initialized();
1024 }
1025 
1026 
1027 void LDPC_Generator_Systematic::encode(const bvec &input, bvec &output)
1028 {
1029  it_assert(is_initialized(), "LDPC_Generator_Systematic::encode(): Systematic "
1030  "generator not set up");
1031  it_assert(input.size() == G.cols(), "LDPC_Generator_Systematic::encode(): "
1032  "Improper input vector size (" << input.size() << " != "
1033  << G.cols() << ")");
1034 
1035  output = concat(input, G * input);
1036 }
1037 
1038 
1039 // ----------------------------------------------------------------------
1040 // BLDPC_Generator
1041 // ----------------------------------------------------------------------
1042 
1044  const std::string type):
1045  LDPC_Generator(type), H_enc(), N(0), M(0), K(0), Z(0)
1046 {
1047  construct(H);
1048 }
1049 
1050 
1051 void BLDPC_Generator::encode(const bvec &input, bvec &output)
1052 {
1053  it_assert(is_initialized(), "BLDPC_Generator::encode(): Cannot encode with not "
1054  "initialized generator");
1055  it_assert(input.size() == K, "BLDPC_Generator::encode(): Input vector "
1056  "length is not equal to K");
1057 
1058  // copy systematic bits first
1059  output = input;
1060  output.set_size(N, true);
1061 
1062  // backward substitution to obtain the first Z parity check bits
1063  for (int k = 0; k < Z; k++) {
1064  for (int j = 0; j < K; j++) {
1065  output(K + k) += H_enc(M - 1 - k, j) * input(j);
1066  }
1067  for (int j = 0; j < k; j++) {
1068  output(K + k) += H_enc(M - 1 - k, K + j) * output(K + j);
1069  }
1070  }
1071 
1072  // forward substitution to obtain the next M-Z parity check bits
1073  for (int k = 0; k < M - Z; k++) {
1074  for (int j = 0; j < K; j++) {
1075  output(K + Z + k) += H_enc(k, j) * input(j);
1076  }
1077  for (int j = K; j < K + Z; j++) {
1078  output(K + Z + k) += H_enc(k, j) * output(j);
1079  }
1080  for (int j = K + Z; j < K + Z + k; j++) {
1081  output(K + Z + k) += H_enc(k, j) * output(j);
1082  }
1083  }
1084 }
1085 
1086 
1088 {
1089  if (H != 0 && H->is_valid()) {
1090  H_enc = H->get_H(); // sparse to dense conversion
1091  Z = H->get_exp_factor();
1092  N = H_enc.cols();
1093  M = H_enc.rows();
1094  K = N - M;
1095 
1096  // ----------------------------------------------------------------------
1097  // STEP 1
1098  // ----------------------------------------------------------------------
1099  // loop over last M-Z columns of matrix H
1100  for (int i = 0; i < M - Z; i += Z) {
1101  // loop over last Z rows of matrix H
1102  for (int j = 0; j < Z; j++) {
1103  // Gaussian elimination by adding particular rows
1104  H_enc.add_rows(M - 1 - j, M - Z - 1 - j - i);
1105  }
1106  }
1107 
1108  // ----------------------------------------------------------------------
1109  // STEP 2
1110  // ----------------------------------------------------------------------
1111  // set first processed row index to M-Z
1112  int r1 = M - Z;
1113  // loop over columns with indexes K .. K+Z-1
1114  for (int c1 = K + Z - 1; c1 >= K; c1--) {
1115  int r2 = r1;
1116  // find the first '1' in column c1
1117  while (H_enc(r2, c1) == 0 && r2 < M - 1)
1118  r2++;
1119  // if necessary, swap rows r1 and r2
1120  if (r2 != r1)
1121  H_enc.swap_rows(r1, r2);
1122 
1123  // look for the other '1' in column c1 and get rid of them
1124  for (r2 = r1 + 1; r2 < M; r2++) {
1125  if (H_enc(r2, c1) == 1) {
1126  // Gaussian elimination by adding particular rows
1127  H_enc.add_rows(r2, r1);
1128  }
1129  }
1130 
1131  // increase first processed row index
1132  r1++;
1133  }
1134 
1135  mark_initialized();
1136  }
1137 }
1138 
1139 
1140 void BLDPC_Generator::save(const std::string& filename) const
1141 {
1143  "BLDPC_Generator::save(): Can not save not initialized generator");
1144  // Every Z-th row of H_enc until M-Z
1145  GF2mat H_T(M / Z - 1, N);
1146  for (int i = 0; i < M / Z - 1; i++) {
1147  H_T.set_row(i, H_enc.get_row(i*Z));
1148  }
1149  // Last Z preprocessed rows of H_enc
1150  GF2mat H_Z = H_enc.get_submatrix(M - Z, 0, M - 1, N - 1);
1151 
1152  it_file f(filename);
1153  int ver;
1154  f >> Name("Fileversion") >> ver;
1155  it_assert(ver == LDPC_binary_file_version, "BLDPC_Generator::save(): "
1156  "Unsupported file format");
1157  f << Name("G_type") << get_type();
1158  f << Name("H_T") << H_T;
1159  f << Name("H_Z") << H_Z;
1160  f << Name("Z") << Z;
1161  f.close();
1162 }
1163 
1164 void BLDPC_Generator::load(const std::string& filename)
1165 {
1166  GF2mat H_T, H_Z;
1167 
1168  it_ifile f(filename);
1169  int ver;
1170  f >> Name("Fileversion") >> ver;
1171  it_assert(ver == LDPC_binary_file_version, "BLDPC_Generator::load(): "
1172  "Unsupported file format");
1173  std::string gen_type;
1174  f >> Name("G_type") >> gen_type;
1175  it_assert(gen_type == get_type(),
1176  "BLDPC_Generator::load(): Wrong generator type");
1177  f >> Name("H_T") >> H_T;
1178  f >> Name("H_Z") >> H_Z;
1179  f >> Name("Z") >> Z;
1180  f.close();
1181 
1182  N = H_T.cols();
1183  M = (H_T.rows() + 1) * Z;
1184  K = N - M;
1185 
1186  H_enc = GF2mat(M - Z, N);
1187  for (int i = 0; i < H_T.rows(); i++) {
1188  for (int j = 0; j < Z; j++) {
1189  for (int k = 0; k < N; k++) {
1190  if (H_T(i, (k / Z)*Z + (k + Z - j) % Z)) {
1191  H_enc.set(i*Z + j, k, 1);
1192  }
1193  }
1194  }
1195  }
1197 
1198  mark_initialized();
1199 }
1200 
1201 
1202 // ----------------------------------------------------------------------
1203 // LDPC_Code
1204 // ----------------------------------------------------------------------
1205 
1206 LDPC_Code::LDPC_Code(): H_defined(false), G_defined(false), dec_method(new std::string),
1207  max_iters(50), psc(true), pisc(false),
1208  llrcalc(LLR_calc_unit()) { set_decoding_method("BP");}
1209 
1211  LDPC_Generator* const G_in,
1212  bool perform_integrity_check):
1213  H_defined(false), G_defined(false), dec_method(new std::string), max_iters(50),
1214  psc(true), pisc(false), llrcalc(LLR_calc_unit())
1215 {
1216  set_decoding_method("BP");
1217  set_code(H, G_in, perform_integrity_check);
1218 }
1219 
1220 LDPC_Code::LDPC_Code(const std::string& filename,
1221  LDPC_Generator* const G_in):
1222  H_defined(false), G_defined(false), dec_method(new std::string), max_iters(50),
1223  psc(true), pisc(false), llrcalc(LLR_calc_unit())
1224 {
1225  set_decoding_method("BP");
1226  load_code(filename, G_in);
1227 }
1228 
1229 
1230 void LDPC_Code::set_code(const LDPC_Parity* const H,
1231  LDPC_Generator* const G_in,
1232  bool perform_integrity_check)
1233 {
1235  setup_decoder();
1236  G = G_in;
1237  if (G != 0) {
1238  G_defined = true;
1239  if (perform_integrity_check) {
1240  integrity_check();
1241  } else {
1242  it_info_debug("LDPC_Code::set_code(): integrity check was not performed");
1243  }
1244  }
1245 }
1246 
1247 void LDPC_Code::load_code(const std::string& filename,
1248  LDPC_Generator* const G_in)
1249 {
1250  it_info_debug("LDPC_Code::load_code(): Loading LDPC codec from "
1251  << filename);
1252 
1253  it_ifile f(filename);
1254  int ver;
1255  f >> Name("Fileversion") >> ver;
1256  it_assert(ver == LDPC_binary_file_version, "LDPC_Code::load_code(): "
1257  "Unsupported file format");
1258  f >> Name("H_defined") >> H_defined;
1259  f >> Name("G_defined") >> G_defined;
1260  f >> Name("nvar") >> nvar;
1261  f >> Name("ncheck") >> ncheck;
1262  f >> Name("C") >> C;
1263  f >> Name("V") >> V;
1264  f >> Name("sumX1") >> sumX1;
1265  f >> Name("sumX2") >> sumX2;
1266  f >> Name("iind") >> iind;
1267  f >> Name("jind") >> jind;
1268  f.close();
1269 
1270  // load generator data
1271  if (G_defined) {
1272  it_assert(G_in != 0, "LDPC_Code::load_code(): Generator object is "
1273  "missing. Can not load the generator data from a file.");
1274  G = G_in;
1275  G->load(filename);
1276  }
1277  else {
1278  G = 0;
1279  it_info_debug("LDPC_Code::load_code(): Generator data not loaded. "
1280  "Generator object will not be used.");
1281  }
1282 
1283  it_info_debug("LDPC_Code::load_code(): Successfully loaded LDPC codec "
1284  "from " << filename);
1285 
1286  setup_decoder();
1287 }
1288 
1289 void LDPC_Code::save_code(const std::string& filename) const
1290 {
1291  it_assert(H_defined, "LDPC_Code::save_to_file(): There is no parity "
1292  "check matrix");
1293  it_info_debug("LDPC_Code::save_to_file(): Saving LDPC codec to "
1294  << filename);
1295 
1296  it_file f;
1297  f.open(filename, true);
1298  f << Name("Fileversion") << LDPC_binary_file_version;
1299  f << Name("H_defined") << H_defined;
1300  f << Name("G_defined") << G_defined;
1301  f << Name("nvar") << nvar;
1302  f << Name("ncheck") << ncheck;
1303  f << Name("C") << C;
1304  f << Name("V") << V;
1305  f << Name("sumX1") << sumX1;
1306  f << Name("sumX2") << sumX2;
1307  f << Name("iind") << iind;
1308  f << Name("jind") << jind;
1309  f.close();
1310 
1311  // save generator data;
1312  if (G_defined)
1313  G->save(filename);
1314  else
1315  it_info_debug("LDPC_Code::save_code(): Missing generator object - "
1316  "generator data not saved");
1317 
1318  it_info_debug("LDPC_Code::save_code(): Successfully saved LDPC codec to "
1319  << filename);
1320 }
1321 
1322 
1323 void LDPC_Code::set_decoding_method(const std::string& method_in)
1324 {
1325  it_assert((method_in == "bp") || (method_in == "BP"),
1326  "LDPC_Code::set_decoding_method(): Not implemented decoding method");
1327  *dec_method = method_in;
1328 }
1329 
1330 void LDPC_Code::set_exit_conditions(int max_iters_in,
1331  bool syndr_check_each_iter,
1332  bool syndr_check_at_start)
1333 {
1334  it_assert(max_iters >= 0, "LDPC_Code::set_nrof_iterations(): Maximum "
1335  "number of iterations can not be negative");
1336  max_iters = max_iters_in;
1337  psc = syndr_check_each_iter;
1338  pisc = syndr_check_at_start;
1339 }
1340 
1341 void LDPC_Code::set_llrcalc(const LLR_calc_unit& llrcalc_in)
1342 {
1343  llrcalc = llrcalc_in;
1344 }
1345 
1346 
1347 void LDPC_Code::encode(const bvec &input, bvec &output)
1348 {
1349  it_assert(G_defined, "LDPC_Code::encode(): LDPC Generator is required "
1350  "for encoding");
1351  G->encode(input, output);
1352  it_assert_debug(syndrome_check(output), "LDPC_Code::encode(): Syndrome "
1353  "check failed");
1354 }
1355 
1356 bvec LDPC_Code::encode(const bvec &input)
1357 {
1358  bvec result;
1359  encode(input, result);
1360  return result;
1361 }
1362 
1363 void LDPC_Code::decode(const vec &llr_in, bvec &syst_bits)
1364 {
1365  QLLRvec qllrin = llrcalc.to_qllr(llr_in);
1366  QLLRvec qllrout;
1367  bp_decode(qllrin, qllrout);
1368  syst_bits = (qllrout.left(nvar - ncheck) < 0);
1369 }
1370 
1371 bvec LDPC_Code::decode(const vec &llr_in)
1372 {
1373  bvec syst_bits;
1374  decode(llr_in, syst_bits);
1375  return syst_bits;
1376 }
1377 
1378 void LDPC_Code::decode_soft_out(const vec &llr_in, vec &llr_out)
1379 {
1380  QLLRvec qllrin = llrcalc.to_qllr(llr_in);
1381  QLLRvec qllrout;
1382  bp_decode(qllrin, qllrout);
1383  llr_out = llrcalc.to_double(qllrout);
1384 }
1385 
1386 vec LDPC_Code::decode_soft_out(const vec &llr_in)
1387 {
1388  vec llr_out;
1389  decode_soft_out(llr_in, llr_out);
1390  return llr_out;
1391 }
1392 
1393 int LDPC_Code::bp_decode(const QLLRvec &LLRin, QLLRvec &LLRout)
1394 {
1395  // Note the IT++ convention that a sure zero corresponds to
1396  // LLR=+infinity
1397  it_assert(H_defined, "LDPC_Code::bp_decode(): Parity check matrix not "
1398  "defined");
1399  it_assert((LLRin.size() == nvar) && (sumX1.size() == nvar)
1400  && (sumX2.size() == ncheck), "LDPC_Code::bp_decode(): Wrong "
1401  "input dimensions");
1402 
1403  if (pisc && syndrome_check(LLRin)) {
1404  LLRout = LLRin;
1405  return 0;
1406  }
1407 
1408  LLRout.set_size(LLRin.size());
1409 
1410  // allocate temporary variables used for the check node update
1411  ivec jj(max_cnd);
1412  QLLRvec m(max_cnd);
1413  QLLRvec ml(max_cnd);
1414  QLLRvec mr(max_cnd);
1415 
1416  // initial step
1417  for (int i = 0; i < nvar; i++) {
1418  int index = i;
1419  for (int j = 0; j < sumX1(i); j++) {
1420  mvc[index] = LLRin(i);
1421  index += nvar;
1422  }
1423  }
1424 
1425  bool is_valid_codeword = false;
1426  int iter = 0;
1427  do {
1428  iter++;
1429  if (nvar >= 100000) { it_info_no_endl_debug("."); }
1430  // --------- Step 1: check to variable nodes ----------
1431  for (int j = 0; j < ncheck; j++) {
1432  // The check node update calculations are hardcoded for degrees
1433  // up to 6. For larger degrees, a general algorithm is used.
1434  switch (sumX2(j)) {
1435  case 0:
1436  it_error("LDPC_Code::bp_decode(): sumX2(j)=0");
1437  case 1:
1438  it_error("LDPC_Code::bp_decode(): sumX2(j)=1");
1439  case 2: {
1440  mcv[j+ncheck] = mvc[jind[j]];
1441  mcv[j] = mvc[jind[j+ncheck]];
1442  break;
1443  }
1444  case 3: {
1445  int j0 = j;
1446  QLLR m0 = mvc[jind[j0]];
1447  int j1 = j0 + ncheck;
1448  QLLR m1 = mvc[jind[j1]];
1449  int j2 = j1 + ncheck;
1450  QLLR m2 = mvc[jind[j2]];
1451  mcv[j0] = llrcalc.Boxplus(m1, m2);
1452  mcv[j1] = llrcalc.Boxplus(m0, m2);
1453  mcv[j2] = llrcalc.Boxplus(m0, m1);
1454  break;
1455  }
1456  case 4: {
1457  int j0 = j;
1458  QLLR m0 = mvc[jind[j0]];
1459  int j1 = j0 + ncheck;
1460  QLLR m1 = mvc[jind[j1]];
1461  int j2 = j1 + ncheck;
1462  QLLR m2 = mvc[jind[j2]];
1463  int j3 = j2 + ncheck;
1464  QLLR m3 = mvc[jind[j3]];
1465  QLLR m01 = llrcalc.Boxplus(m0, m1);
1466  QLLR m23 = llrcalc.Boxplus(m2, m3);
1467  mcv[j0] = llrcalc.Boxplus(m1, m23);
1468  mcv[j1] = llrcalc.Boxplus(m0, m23);
1469  mcv[j2] = llrcalc.Boxplus(m01, m3);
1470  mcv[j3] = llrcalc.Boxplus(m01, m2);
1471  break;
1472  }
1473  case 5: {
1474  int j0 = j;
1475  QLLR m0 = mvc[jind[j0]];
1476  int j1 = j0 + ncheck;
1477  QLLR m1 = mvc[jind[j1]];
1478  int j2 = j1 + ncheck;
1479  QLLR m2 = mvc[jind[j2]];
1480  int j3 = j2 + ncheck;
1481  QLLR m3 = mvc[jind[j3]];
1482  int j4 = j3 + ncheck;
1483  QLLR m4 = mvc[jind[j4]];
1484  QLLR m01 = llrcalc.Boxplus(m0, m1);
1485  QLLR m02 = llrcalc.Boxplus(m01, m2);
1486  QLLR m34 = llrcalc.Boxplus(m3, m4);
1487  QLLR m24 = llrcalc.Boxplus(m2, m34);
1488  mcv[j0] = llrcalc.Boxplus(m1, m24);
1489  mcv[j1] = llrcalc.Boxplus(m0, m24);
1490  mcv[j2] = llrcalc.Boxplus(m01, m34);
1491  mcv[j3] = llrcalc.Boxplus(m02, m4);
1492  mcv[j4] = llrcalc.Boxplus(m02, m3);
1493  break;
1494  }
1495  case 6: {
1496  int j0 = j;
1497  QLLR m0 = mvc[jind[j0]];
1498  int j1 = j0 + ncheck;
1499  QLLR m1 = mvc[jind[j1]];
1500  int j2 = j1 + ncheck;
1501  QLLR m2 = mvc[jind[j2]];
1502  int j3 = j2 + ncheck;
1503  QLLR m3 = mvc[jind[j3]];
1504  int j4 = j3 + ncheck;
1505  QLLR m4 = mvc[jind[j4]];
1506  int j5 = j4 + ncheck;
1507  QLLR m5 = mvc[jind[j5]];
1508  QLLR m01 = llrcalc.Boxplus(m0, m1);
1509  QLLR m23 = llrcalc.Boxplus(m2, m3);
1510  QLLR m45 = llrcalc.Boxplus(m4, m5);
1511  QLLR m03 = llrcalc.Boxplus(m01, m23);
1512  QLLR m25 = llrcalc.Boxplus(m23, m45);
1513  QLLR m0145 = llrcalc.Boxplus(m01, m45);
1514  mcv[j0] = llrcalc.Boxplus(m1, m25);
1515  mcv[j1] = llrcalc.Boxplus(m0, m25);
1516  mcv[j2] = llrcalc.Boxplus(m0145, m3);
1517  mcv[j3] = llrcalc.Boxplus(m0145, m2);
1518  mcv[j4] = llrcalc.Boxplus(m03, m5);
1519  mcv[j5] = llrcalc.Boxplus(m03, m4);
1520  break;
1521  }
1522  default: {
1523  int nodes = sumX2(j);
1524  if( nodes > max_cnd ) {
1525  std::ostringstream m_sout;
1526  m_sout << "check node degrees >" << max_cnd << " not supported in this version";
1527  it_error( m_sout.str() );
1528  }
1529 
1530  nodes--;
1531  jj[0] = j;
1532  m[0] = mvc[jind[jj[0]]];
1533  for(int i = 1; i <= nodes; i++ ) {
1534  jj[i] = jj[i-1] + ncheck;
1535  m[i] = mvc[jind[jj[i]]];
1536  }
1537 
1538  // compute partial sums from the left and from the right
1539  ml[0] = m[0];
1540  mr[0] = m[nodes];
1541  for(int i = 1; i < nodes; i++ ) {
1542  ml[i] = llrcalc.Boxplus( ml[i-1], m[i] );
1543  mr[i] = llrcalc.Boxplus( mr[i-1], m[nodes-i] );
1544  }
1545 
1546  // merge partial sums
1547  mcv[jj[0]] = mr[nodes-1];
1548  mcv[jj[nodes]] = ml[nodes-1];
1549  for(int i = 1; i < nodes; i++ )
1550  mcv[jj[i]] = llrcalc.Boxplus( ml[i-1], mr[nodes-1-i] );
1551  }
1552  } // switch statement
1553  }
1554 
1555  // step 2: variable to check nodes
1556  for (int i = 0; i < nvar; i++) {
1557  switch (sumX1(i)) {
1558  case 0:
1559  it_error("LDPC_Code::bp_decode(): sumX1(i)=0");
1560  case 1: {
1561  /* This case is rare but apparently occurs for codes used in
1562  the DVB-T2 standard.
1563  */
1564  QLLR m0 = mcv[iind[i]];
1565  mvc[i] = LLRin(i);
1566  LLRout(i) = LLRin(i) + m0;
1567  break;
1568  }
1569  case 2: {
1570  QLLR m0 = mcv[iind[i]];
1571  int i1 = i + nvar;
1572  QLLR m1 = mcv[iind[i1]];
1573  mvc[i] = LLRin(i) + m1;
1574  mvc[i1] = LLRin(i) + m0;
1575  LLRout(i) = mvc[i1] + m1;
1576  break;
1577  }
1578  case 3: {
1579  int i0 = i;
1580  QLLR m0 = mcv[iind[i0]];
1581  int i1 = i0 + nvar;
1582  QLLR m1 = mcv[iind[i1]];
1583  int i2 = i1 + nvar;
1584  QLLR m2 = mcv[iind[i2]];
1585  LLRout(i) = LLRin(i) + m0 + m1 + m2;
1586  mvc[i0] = LLRout(i) - m0;
1587  mvc[i1] = LLRout(i) - m1;
1588  mvc[i2] = LLRout(i) - m2;
1589  break;
1590  }
1591  case 4: {
1592  int i0 = i;
1593  QLLR m0 = mcv[iind[i0]];
1594  int i1 = i0 + nvar;
1595  QLLR m1 = mcv[iind[i1]];
1596  int i2 = i1 + nvar;
1597  QLLR m2 = mcv[iind[i2]];
1598  int i3 = i2 + nvar;
1599  QLLR m3 = mcv[iind[i3]];
1600  LLRout(i) = LLRin(i) + m0 + m1 + m2 + m3;
1601  mvc[i0] = LLRout(i) - m0;
1602  mvc[i1] = LLRout(i) - m1;
1603  mvc[i2] = LLRout(i) - m2;
1604  mvc[i3] = LLRout(i) - m3;
1605  break;
1606  }
1607  default: { // differential update
1608  QLLR mvc_temp = LLRin(i);
1609  int index_iind = i; // tracks i+jp*nvar
1610  for (int jp = 0; jp < sumX1(i); jp++) {
1611  mvc_temp += mcv[iind[index_iind]];
1612  index_iind += nvar;
1613  }
1614  LLRout(i) = mvc_temp;
1615  index_iind = i; // tracks i+j*nvar
1616  for (int j = 0; j < sumX1[i]; j++) {
1617  mvc[index_iind] = mvc_temp - mcv[iind[index_iind]];
1618  index_iind += nvar;
1619  }
1620  }
1621  }
1622  }
1623 
1624  if (psc && syndrome_check(LLRout)) {
1625  is_valid_codeword = true;
1626  break;
1627  }
1628  }
1629  while (iter < max_iters);
1630 
1631  if (nvar >= 100000) { it_info_debug(""); }
1632  return (is_valid_codeword ? iter : -iter);
1633 }
1634 
1635 
1636 bool LDPC_Code::syndrome_check(const bvec &x) const
1637 {
1638  QLLRvec llr = 1 - 2 * to_ivec(x);
1639  return syndrome_check(llr);
1640 }
1641 
1642 bool LDPC_Code::syndrome_check(const QLLRvec &LLR) const
1643 {
1644  // Please note the IT++ convention that a sure zero corresponds to
1645  // LLR=+infinity
1646  int i, j, synd, vi;
1647 
1648  for (j = 0; j < ncheck; j++) {
1649  synd = 0;
1650  int vind = j; // tracks j+i*ncheck
1651  for (i = 0; i < sumX2(j); i++) {
1652  vi = V(vind);
1653  if (LLR(vi) < 0) {
1654  synd++;
1655  }
1656  vind += ncheck;
1657  }
1658  if ((synd&1) == 1) {
1659  return false; // codeword is invalid
1660  }
1661  }
1662  return true; // codeword is valid
1663 };
1664 
1665 QLLRvec LDPC_Code::soft_syndrome_check(const QLLRvec &LLR) const
1666 {
1667  QLLRvec result(ncheck);
1668  int i,j,vi,vind;
1669 
1670  for (j=0; j<ncheck; j++) {
1671  result(j)=-QLLR_MAX;
1672  vind=j;
1673  for (i=0; i<sumX2(j); i++) {
1674  vi=V(vind);
1675  result(j) = llrcalc.Boxplus(LLR(vi),result(j));
1676  vind += ncheck;
1677  }
1678  }
1679 
1680  return result;
1681 }
1682 
1683 // ----------------------------------------------------------------------
1684 // LDPC_Code private methods
1685 // ----------------------------------------------------------------------
1686 
1688 {
1689  // copy basic parameters
1690  sumX1 = Hmat->sumX1;
1691  sumX2 = Hmat->sumX2;
1692  nvar = Hmat->nvar; //get_nvar();
1693  ncheck = Hmat->ncheck; //get_ncheck();
1694  int cmax = max(sumX1);
1695  int vmax = max(sumX2);
1696 
1697  // decoder parameterization
1698  V = zeros_i(ncheck * vmax);
1699  C = zeros_i(cmax * nvar);
1700  jind = zeros_i(ncheck * vmax);
1701  iind = zeros_i(nvar * cmax);
1702 
1703  it_info_debug("LDPC_Code::decoder_parameterization(): Computations "
1704  "- phase 1");
1705  for (int i = 0; i < nvar; i++) {
1706  ivec coli = Hmat->get_col(i).get_nz_indices();
1707  for (int j0 = 0; j0 < length(coli); j0++) {
1708  C(j0 + cmax*i) = coli(j0);
1709  }
1710  }
1711 
1712  it_info_debug("LDPC_Code::decoder_parameterization(): Computations "
1713  "- phase 2");
1714  it_info_debug("Computing decoder parameterization. Phase 2");
1715  for (int j = 0; j < ncheck; j++) {
1716  ivec rowj = Hmat->get_row(j).get_nz_indices();
1717  for (int i0 = 0; i0 < length(rowj); i0++) {
1718  V(j + ncheck*i0) = rowj(i0);
1719  }
1720  }
1721 
1722  it_info_debug("LDPC_Code::decoder_parameterization(): Computations "
1723  "- phase 3");
1724  it_info_debug("Computing decoder parameterization. Phase 3.");
1725  for (int j = 0; j < ncheck; j++) {
1726  for (int ip = 0; ip < sumX2(j); ip++) {
1727  int vip = V(j + ip * ncheck);
1728  int k = 0;
1729  while (1 == 1) {
1730  if (C(k + vip*cmax) == j) {
1731  break;
1732  }
1733  k++;
1734  }
1735  jind(j + ip*ncheck) = vip + k * nvar;
1736  }
1737  }
1738 
1739  it_info_debug("LDPC_Code::decoder_parameterization(): Computations "
1740  "- phase 4");
1741  for (int i = 0; i < nvar; i++) {
1742  for (int jp = 0; jp < sumX1(i); jp++) {
1743  int cjp = C(jp + i * cmax);
1744  int k = 0;
1745  while (1 == 1) {
1746  if (V(cjp + k*ncheck) == i) {break; }
1747  k++;
1748  }
1749  iind(i + jp*nvar) = cjp + k * ncheck;
1750  }
1751  }
1752 
1753  H_defined = true;
1754 }
1755 
1756 
1758 {
1759  if (H_defined) {
1760  mcv.set_size(max(sumX2) * ncheck);
1761  mvc.set_size(max(sumX1) * nvar);
1762  }
1763 }
1764 
1765 
1767 {
1768  if (G_defined) {
1769  it_info_debug("LDPC_Code::integrity_check(): Checking integrity of "
1770  "the LDPC_Parity and LDPC_Generator data");
1771  bvec bv(nvar - ncheck), cw;
1772  bv.clear();
1773  bv(0) = 1;
1774  for (int i = 0; i < nvar - ncheck; i++) {
1775  G->encode(bv, cw);
1777  "LDPC_Code::integrity_check(): Syndrome check failed");
1778  bv.shift_right(bv(nvar - ncheck - 1));
1779  }
1780  }
1781  else {
1782  it_info_debug("LDPC_Code::integrity_check(): No generator defined "
1783  "- no check performed");
1784  }
1785 }
1786 
1787 // ----------------------------------------------------------------------
1788 // Related functions
1789 // ----------------------------------------------------------------------
1790 
1791 std::ostream &operator<<(std::ostream &os, const LDPC_Code &C)
1792 {
1793  ivec rdeg = zeros_i(max(C.sumX2) + 1);
1794  for (int i = 0; i < C.ncheck; i++) {
1795  rdeg(C.sumX2(i))++;
1796  }
1797 
1798  ivec cdeg = zeros_i(max(C.sumX1) + 1);
1799  for (int j = 0; j < C.nvar; j++) {
1800  cdeg(C.sumX1(j))++;
1801  }
1802 
1803  os << "--- LDPC codec ----------------------------------\n"
1804  << "Nvar : " << C.get_nvar() << "\n"
1805  << "Ncheck : " << C.get_ncheck() << "\n"
1806  << "Rate : " << C.get_rate() << "\n"
1807  << "Column degrees (node perspective): " << cdeg << "\n"
1808  << "Row degrees (node perspective): " << rdeg << "\n"
1809  << "-------------------------------------------------\n"
1810  << "Decoder parameters:\n"
1811  << " - method : " << C.get_decoding_method() << "\n"
1812  << " - max. iterations : " << C.max_iters << "\n"
1813  << " - syndrome check at each iteration : " << C.psc << "\n"
1814  << " - syndrome check at start : " << C.pisc << "\n"
1815  << "-------------------------------------------------\n"
1816  << C.llrcalc << "\n";
1817  return os;
1818 }
1819 
1820 } // namespace itpp
SourceForge Logo

Generated on Sat Jul 6 2013 10:54:22 for IT++ by Doxygen 1.8.2