40 inline double max(
double x,
double y) {
return std::max(x, y); }
45 int Rec_Syst_Conv_Code::calc_state_transition(
const int instate,
const int input, ivec &parity)
47 int i, j, in = 0, temp = (gen_pol_rev(0) & (instate << 1)), parity_temp, parity_bit;
49 for (i = 0; i < K; i++) {
55 parity.set_size(n - 1,
false);
56 for (j = 0; j < (n - 1); j++) {
57 parity_temp = ((instate << 1) + in) & gen_pol_rev(j + 1);
59 for (i = 0; i < K; i++) {
60 parity_bit = (parity_temp & 1) ^ parity_bit;
61 parity_temp = parity_temp >> 1;
63 parity(j) = parity_bit;
65 return in + ((instate << 1) & ((1 << m) - 1));
74 K = constraint_length;
78 gen_pol_rev.set_size(n,
false);
79 for (
int i = 0; i < n; i++) {
84 state_trans.set_size(Nstates, 2,
false);
85 rev_state_trans.set_size(Nstates, 2,
false);
86 output_parity.set_size(Nstates, 2*(n - 1),
false);
87 rev_output_parity.set_size(Nstates, 2*(n - 1),
false);
90 for (s_prim = 0; s_prim < Nstates; s_prim++) {
91 s0 = calc_state_transition(s_prim, 0, p0);
92 state_trans(s_prim, 0) = s0;
93 rev_state_trans(s0, 0) = s_prim;
94 for (j = 0; j < (n - 1); j++) {
95 output_parity(s_prim, 2*j + 0) = p0(j);
96 rev_output_parity(s0, 2*j + 0) = p0(j);
99 s1 = calc_state_transition(s_prim, 1, p1);
100 state_trans(s_prim, 1) = s1;
101 rev_state_trans(s1, 1) = s_prim;
102 for (j = 0; j < (n - 1); j++) {
103 output_parity(s_prim, 2*j + 1) = p1(j);
104 rev_output_parity(s1, 2*j + 1) = p1(j);
126 int i, j,
length = input.size(), target_state;
127 parity_bits.set_size(length + m, n - 1,
false);
128 tail.set_size(m,
false);
131 for (i = 0; i <
length; i++) {
132 for (j = 0; j < (n - 1); j++) {
133 parity_bits(i, j) = output_parity(encoder_state, 2 * j +
int(input(i)));
135 encoder_state = state_trans(encoder_state,
int(input(i)));
139 for (i = 0; i < m; i++) {
140 target_state = (encoder_state << 1) & ((1 << m) - 1);
141 if (state_trans(encoder_state, 0) == target_state) { tail(i) =
bin(0); }
142 else { tail(i) =
bin(1); }
143 for (j = 0; j < (n - 1); j++) {
144 parity_bits(length + i, j) = output_parity(encoder_state, 2 * j +
int(tail(i)));
146 encoder_state = target_state;
153 int i, j,
length = input.size();
154 parity_bits.set_size(length, n - 1,
false);
157 for (i = 0; i <
length; i++) {
158 for (j = 0; j < (n - 1); j++) {
159 parity_bits(i, j) = output_parity(encoder_state, 2 * j +
int(input(i)));
161 encoder_state = state_trans(encoder_state,
int(input(i)));
167 vec &extrinsic_output,
bool in_terminated)
169 double gamma_k_e, nom, den, temp0, temp1, exp_temp0, exp_temp1;
170 int j, s0, s1, k, kk, s, s_prim, s_prim0, s_prim1, block_length = rec_systematic.length();
173 alpha.set_size(Nstates, block_length + 1,
false);
174 beta.set_size(Nstates, block_length + 1,
false);
175 gamma.set_size(2*Nstates, block_length + 1,
false);
176 denom.set_size(block_length + 1,
false);
178 extrinsic_output.set_size(block_length,
false);
180 if (in_terminated) { terminated =
true; }
183 for (k = 1; k <= block_length; k++) {
185 for (s_prim = 0; s_prim < Nstates; s_prim++) {
188 for (j = 0; j < (n - 1); j++) {
189 exp_temp0 += 0.5 * Lc * rec_parity(kk, j) * double(1 - 2 * output_parity(s_prim, 2 * j + 0));
190 exp_temp1 += 0.5 * Lc * rec_parity(kk, j) * double(1 - 2 * output_parity(s_prim, 2 * j + 1));
213 gamma(2*s_prim + 0, k) =
trunc_exp(0.5 * (extrinsic_input(kk) + Lc * rec_systematic(kk)) + exp_temp0);
214 gamma(2*s_prim + 1, k) =
trunc_exp(-0.5 * (extrinsic_input(kk) + Lc * rec_systematic(kk)) + exp_temp1);
219 alpha.set_col(0,
zeros(Nstates));
223 for (k = 1; k <= block_length; k++) {
224 for (s = 0; s < Nstates; s++) {
225 s_prim0 = rev_state_trans(s, 0);
226 s_prim1 = rev_state_trans(s, 1);
227 temp0 = alpha(s_prim0, k - 1) * gamma(2 * s_prim0 + 0, k);
228 temp1 = alpha(s_prim1, k - 1) * gamma(2 * s_prim1 + 1, k);
229 alpha(s, k) = temp0 + temp1;
230 denom(k) += temp0 + temp1;
232 alpha.set_col(k, alpha.get_col(k) / denom(k));
237 beta.set_col(block_length,
zeros(Nstates));
238 beta(0, block_length) = 1.0;
241 beta.set_col(block_length, alpha.get_col(block_length));
245 for (k = block_length; k >= 2; k--) {
246 for (s_prim = 0; s_prim < Nstates; s_prim++) {
247 s0 = state_trans(s_prim, 0);
248 s1 = state_trans(s_prim, 1);
249 beta(s_prim, k - 1) = (beta(s0, k) * gamma(2 * s_prim + 0, k)) + (beta(s1, k) * gamma(2 * s_prim + 1, k));
251 beta.set_col(k - 1, beta.get_col(k - 1) / denom(k));
255 for (k = 1; k <= block_length; k++) {
259 for (s_prim = 0; s_prim < Nstates; s_prim++) {
260 s0 = state_trans(s_prim, 0);
261 s1 = state_trans(s_prim, 1);
264 for (j = 0; j < (n - 1); j++) {
265 exp_temp0 += 0.5 * Lc * rec_parity(kk, j) * double(1 - 2 * output_parity(s_prim, 2 * j + 0));
266 exp_temp1 += 0.5 * Lc * rec_parity(kk, j) * double(1 - 2 * output_parity(s_prim, 2 * j + 1));
270 nom += alpha(s_prim, k - 1) * gamma_k_e * beta(s0, k);
274 den += alpha(s_prim, k - 1) * gamma_k_e * beta(s1, k);
277 extrinsic_output(kk) =
trunc_log(nom / den);
283 const vec &extrinsic_input, vec &extrinsic_output,
bool in_terminated, std::string metric)
285 if (metric ==
"TABLE") {
289 QLLRvec rec_systematic_q = llrcalc.
to_qllr(rec_systematic);
290 QLLRmat rec_parity_q = llrcalc.
to_qllr(rec_parity);
291 QLLRvec extrinsic_input_q = llrcalc.
to_qllr(extrinsic_input);
292 QLLRvec extrinsic_output_q(
length(extrinsic_output));
293 log_decode(rec_systematic_q, rec_parity_q, extrinsic_input_q,
294 extrinsic_output_q, in_terminated);
295 extrinsic_output = llrcalc.
to_double(extrinsic_output_q);
299 double nom, den, exp_temp0, exp_temp1, rp, temp0, temp1;
300 int i, j, s0, s1, k, kk, l, s, s_prim, s_prim0, s_prim1, block_length = rec_systematic.length();
307 it_error(
"Rec_Syst_Conv_Code::log_decode: Illegal metric parameter");
310 alpha.set_size(Nstates, block_length + 1,
false);
311 beta.set_size(Nstates, block_length + 1,
false);
312 gamma.set_size(2*Nstates, block_length + 1,
false);
313 extrinsic_output.set_size(block_length,
false);
314 denom.set_size(block_length + 1,
false);
315 for (k = 0; k <= block_length; k++) { denom(k) = -infinity; }
317 if (in_terminated) { terminated =
true; }
321 "Rec_Syst_Conv_Code::log_decode: This function assumes that Lc = 1.0. Please use proper scaling of the input data");
324 for (k = 1; k <= block_length; k++) {
326 for (s_prim = 0; s_prim < Nstates; s_prim++) {
329 for (j = 0; j < (n - 1); j++) {
330 rp = rec_parity(kk, j);
331 if (output_parity(s_prim , 2*j + 0) == 0) { exp_temp0 += rp; }
332 else { exp_temp0 -= rp; }
333 if (output_parity(s_prim , 2*j + 1) == 0) { exp_temp1 += rp; }
334 else { exp_temp1 -= rp; }
336 gamma(2*s_prim + 0, k) = 0.5 * ((extrinsic_input(kk) + rec_systematic(kk)) + exp_temp0);
337 gamma(2*s_prim + 1, k) = -0.5 * ((extrinsic_input(kk) + rec_systematic(kk)) - exp_temp1);
342 for (j = 1; j < Nstates; j++) { alpha(j, 0) = -infinity; }
346 for (k = 1; k <= block_length; k++) {
347 for (s = 0; s < Nstates; s++) {
348 s_prim0 = rev_state_trans(s, 0);
349 s_prim1 = rev_state_trans(s, 1);
350 temp0 = alpha(s_prim0, k - 1) + gamma(2 * s_prim0 + 0, k);
351 temp1 = alpha(s_prim1, k - 1) + gamma(2 * s_prim1 + 1, k);
352 alpha(s, k) =
com_log(temp0, temp1);
353 denom(k) =
com_log(alpha(s, k), denom(k));
356 for (l = 0; l < Nstates; l++) { alpha(l, k) -= denom(k); }
361 for (i = 1; i < Nstates; i++) { beta(i, block_length) = -infinity; }
362 beta(0, block_length) = 0.0;
365 beta.set_col(block_length, alpha.get_col(block_length));
369 for (k = block_length; k >= 1; k--) {
370 for (s_prim = 0; s_prim < Nstates; s_prim++) {
371 s0 = state_trans(s_prim, 0);
372 s1 = state_trans(s_prim, 1);
373 beta(s_prim, k - 1) =
com_log(beta(s0, k) + gamma(2 * s_prim + 0, k) , beta(s1, k) + gamma(2 * s_prim + 1, k));
376 for (l = 0; l < Nstates; l++) { beta(l, k - 1) -= denom(k); }
380 for (k = 1; k <= block_length; k++) {
384 for (s_prim = 0; s_prim < Nstates; s_prim++) {
385 s0 = state_trans(s_prim, 0);
386 s1 = state_trans(s_prim, 1);
389 for (j = 0; j < (n - 1); j++) {
390 rp = rec_parity(kk, j);
391 if (output_parity(s_prim , 2*j + 0) == 0) { exp_temp0 += rp; }
392 else { exp_temp0 -= rp; }
393 if (output_parity(s_prim , 2*j + 1) == 0) { exp_temp1 += rp; }
394 else { exp_temp1 -= rp; }
396 nom =
com_log(nom, alpha(s_prim, kk) + 0.5 * exp_temp0 + beta(s0, k));
397 den =
com_log(den, alpha(s_prim, kk) + 0.5 * exp_temp1 + beta(s1, k));
399 extrinsic_output(kk) = nom - den;
405 const vec &extrinsic_input, vec &extrinsic_output,
bool in_terminated, std::string metric)
407 if (metric ==
"TABLE") {
408 QLLRvec rec_systematic_q = llrcalc.
to_qllr(rec_systematic);
409 QLLRvec rec_parity_q = llrcalc.
to_qllr(rec_parity);
410 QLLRvec extrinsic_input_q = llrcalc.
to_qllr(extrinsic_input);
411 QLLRvec extrinsic_output_q(
length(extrinsic_output));
412 log_decode_n2(rec_systematic_q, rec_parity_q, extrinsic_input_q,
413 extrinsic_output_q, in_terminated);
414 extrinsic_output = llrcalc.
to_double(extrinsic_output_q);
419 double nom, den, exp_temp0, exp_temp1, rp;
420 int k, kk, l, s, s_prim, s_prim0, s_prim1, block_length = rec_systematic.length();
421 int ext_info_length = extrinsic_input.length();
429 it_error(
"Rec_Syst_Conv_Code::log_decode_n2: Illegal metric parameter");
432 alpha.set_size(Nstates, block_length + 1,
false);
433 beta.set_size(Nstates, block_length + 1,
false);
434 gamma.set_size(2*Nstates, block_length + 1,
false);
435 extrinsic_output.set_size(ext_info_length,
false);
438 if (in_terminated) { terminated =
true; }
442 "Rec_Syst_Conv_Code::log_decode_n2: This function assumes that Lc = 1.0. Please use proper scaling of the input data");
445 for (s = 1; s < Nstates; s++) { alpha(s, 0) = -infinity; }
449 for (k = 1; k <= block_length; k++) {
451 if (kk < ext_info_length) {
452 ex = 0.5 * (extrinsic_input(kk) + rec_systematic(kk));
455 ex = 0.5 * rec_systematic(kk);
457 rp = 0.5 * rec_parity(kk);
458 for (s = 0; s < Nstates; s++) {
459 s_prim0 = rev_state_trans(s, 0);
460 s_prim1 = rev_state_trans(s, 1);
461 if (output_parity(s_prim0 , 0)) { exp_temp0 = -rp; }
462 else { exp_temp0 = rp; }
463 if (output_parity(s_prim1 , 1)) { exp_temp1 = -rp; }
464 else { exp_temp1 = rp; }
465 gamma(2*s_prim0 , k) = ex + exp_temp0;
466 gamma(2*s_prim1 + 1, k) = -ex + exp_temp1;
467 alpha(s, k) =
com_log(alpha(s_prim0, kk) + gamma(2 * s_prim0 , k),
468 alpha(s_prim1, kk) + gamma(2 * s_prim1 + 1, k));
472 for (l = 0; l < Nstates; l++) { alpha(l, k) -=
norm; }
477 for (s = 1; s < Nstates; s++) { beta(s, block_length) = -infinity; }
478 beta(0, block_length) = 0.0;
481 beta.set_col(block_length, alpha.get_col(block_length));
485 for (k = block_length; k >= 1; k--) {
487 for (s_prim = 0; s_prim < Nstates; s_prim++) {
488 beta(s_prim, kk) =
com_log(beta(state_trans(s_prim, 0), k) + gamma(2 * s_prim, k),
489 beta(state_trans(s_prim, 1), k) + gamma(2 * s_prim + 1, k));
492 for (l = 0; l < Nstates; l++) { beta(l, k) -=
norm; }
496 for (k = 1; k <= block_length; k++) {
498 if (kk < ext_info_length) {
501 rp = 0.5 * rec_parity(kk);
502 for (s_prim = 0; s_prim < Nstates; s_prim++) {
503 if (output_parity(s_prim , 0)) { exp_temp0 = -rp; }
504 else { exp_temp0 = rp; }
505 if (output_parity(s_prim , 1)) { exp_temp1 = -rp; }
506 else { exp_temp1 = rp; }
507 nom =
com_log(nom, alpha(s_prim, kk) + exp_temp0 + beta(state_trans(s_prim, 0), k));
508 den =
com_log(den, alpha(s_prim, kk) + exp_temp1 + beta(state_trans(s_prim, 1), k));
510 extrinsic_output(kk) = nom - den;
518 const QLLRvec &extrinsic_input,
519 QLLRvec &extrinsic_output,
bool in_terminated)
522 int nom, den, exp_temp0, exp_temp1, rp, temp0, temp1;
523 int i, j, s0, s1, k, kk, l, s, s_prim, s_prim0, s_prim1, block_length = rec_systematic.length();
526 alpha_q.set_size(Nstates, block_length + 1,
false);
527 beta_q.set_size(Nstates, block_length + 1,
false);
528 gamma_q.set_size(2*Nstates, block_length + 1,
false);
529 extrinsic_output.set_size(block_length,
false);
530 denom_q.set_size(block_length + 1,
false);
531 for (k = 0; k <= block_length; k++) { denom_q(k) = -QLLR_MAX; }
533 if (in_terminated) { terminated =
true; }
537 "Rec_Syst_Conv_Code::log_decode: This function assumes that Lc = 1.0. Please use proper scaling of the input data");
540 for (k = 1; k <= block_length; k++) {
542 for (s_prim = 0; s_prim < Nstates; s_prim++) {
545 for (j = 0; j < (n - 1); j++) {
546 rp = rec_parity(kk, j);
547 if (output_parity(s_prim , 2*j + 0) == 0) { exp_temp0 += rp; }
548 else { exp_temp0 -= rp; }
549 if (output_parity(s_prim , 2*j + 1) == 0) { exp_temp1 += rp; }
550 else { exp_temp1 -= rp; }
553 gamma_q(2*s_prim + 0, k) = ((extrinsic_input(kk) + rec_systematic(kk)) + exp_temp0) / 2;
554 gamma_q(2*s_prim + 1, k) = - ((extrinsic_input(kk) + rec_systematic(kk)) - exp_temp1) / 2;
559 for (j = 1; j < Nstates; j++) { alpha_q(j, 0) = -QLLR_MAX; }
563 for (k = 1; k <= block_length; k++) {
564 for (s = 0; s < Nstates; s++) {
565 s_prim0 = rev_state_trans(s, 0);
566 s_prim1 = rev_state_trans(s, 1);
567 temp0 = alpha_q(s_prim0, k - 1) + gamma_q(2 * s_prim0 + 0, k);
568 temp1 = alpha_q(s_prim1, k - 1) + gamma_q(2 * s_prim1 + 1, k);
571 alpha_q(s, k) = llrcalc.
jaclog(temp0, temp1);
572 denom_q(k) = llrcalc.
jaclog(alpha_q(s, k), denom_q(k));
575 for (l = 0; l < Nstates; l++) { alpha_q(l, k) -= denom_q(k); }
580 for (i = 1; i < Nstates; i++) { beta_q(i, block_length) = -QLLR_MAX; }
581 beta_q(0, block_length) = 0;
584 beta_q.set_col(block_length, alpha_q.get_col(block_length));
588 for (k = block_length; k >= 1; k--) {
589 for (s_prim = 0; s_prim < Nstates; s_prim++) {
590 s0 = state_trans(s_prim, 0);
591 s1 = state_trans(s_prim, 1);
593 beta_q(s_prim, k - 1) = llrcalc.
jaclog(beta_q(s0, k) + gamma_q(2 * s_prim + 0, k) , beta_q(s1, k) + gamma_q(2 * s_prim + 1, k));
596 for (l = 0; l < Nstates; l++) { beta_q(l, k - 1) -= denom_q(k); }
600 for (k = 1; k <= block_length; k++) {
604 for (s_prim = 0; s_prim < Nstates; s_prim++) {
605 s0 = state_trans(s_prim, 0);
606 s1 = state_trans(s_prim, 1);
609 for (j = 0; j < (n - 1); j++) {
610 rp = rec_parity(kk, j);
611 if (output_parity(s_prim , 2*j + 0) == 0) { exp_temp0 += rp; }
612 else { exp_temp0 -= rp; }
613 if (output_parity(s_prim , 2*j + 1) == 0) { exp_temp1 += rp; }
614 else { exp_temp1 -= rp; }
618 nom = llrcalc.
jaclog(nom, alpha_q(s_prim, kk) + exp_temp0 / 2 + beta_q(s0, k));
619 den = llrcalc.
jaclog(den, alpha_q(s_prim, kk) + exp_temp1 / 2 + beta_q(s1, k));
621 extrinsic_output(kk) = nom - den;
629 const QLLRvec &rec_parity,
630 const QLLRvec &extrinsic_input,
631 QLLRvec &extrinsic_output,
634 int nom, den, exp_temp0, exp_temp1, rp;
635 int k, kk, l, s, s_prim, s_prim0, s_prim1, block_length = rec_systematic.length();
636 int ext_info_length = extrinsic_input.length();
641 alpha_q.set_size(Nstates, block_length + 1,
false);
642 beta_q.set_size(Nstates, block_length + 1,
false);
643 gamma_q.set_size(2*Nstates, block_length + 1,
false);
644 extrinsic_output.set_size(ext_info_length,
false);
647 if (in_terminated) { terminated =
true; }
651 "Rec_Syst_Conv_Code::log_decode_n2: This function assumes that Lc = 1.0. Please use proper scaling of the input data");
654 for (s = 1; s < Nstates; s++) { alpha_q(s, 0) = -QLLR_MAX; }
658 for (k = 1; k <= block_length; k++) {
660 if (kk < ext_info_length) {
661 ex = (extrinsic_input(kk) + rec_systematic(kk)) / 2;
664 ex = rec_systematic(kk) / 2;
666 rp = rec_parity(kk) / 2;
667 for (s = 0; s < Nstates; s++) {
668 s_prim0 = rev_state_trans(s, 0);
669 s_prim1 = rev_state_trans(s, 1);
670 if (output_parity(s_prim0 , 0)) { exp_temp0 = -rp; }
671 else { exp_temp0 = rp; }
672 if (output_parity(s_prim1 , 1)) { exp_temp1 = -rp; }
673 else { exp_temp1 = rp; }
674 gamma_q(2*s_prim0 , k) = ex + exp_temp0;
675 gamma_q(2*s_prim1 + 1, k) = -ex + exp_temp1;
676 alpha_q(s, k) = llrcalc.
jaclog(alpha_q(s_prim0, kk) + gamma_q(2 * s_prim0 , k),
677 alpha_q(s_prim1, kk) + gamma_q(2 * s_prim1 + 1, k));
680 norm = alpha_q(0, k);
681 for (l = 0; l < Nstates; l++) { alpha_q(l, k) -=
norm; }
686 for (s = 1; s < Nstates; s++) { beta_q(s, block_length) = -QLLR_MAX; }
687 beta_q(0, block_length) = 0;
690 beta_q.set_col(block_length, alpha_q.get_col(block_length));
694 for (k = block_length; k >= 1; k--) {
696 for (s_prim = 0; s_prim < Nstates; s_prim++) {
697 beta_q(s_prim, kk) = llrcalc.
jaclog(beta_q(state_trans(s_prim, 0), k) + gamma_q(2 * s_prim, k),
698 beta_q(state_trans(s_prim, 1), k) + gamma_q(2 * s_prim + 1, k));
701 for (l = 0; l < Nstates; l++) { beta_q(l, k) -=
norm; }
705 for (k = 1; k <= block_length; k++) {
707 if (kk < ext_info_length) {
710 rp = rec_parity(kk) / 2;
711 for (s_prim = 0; s_prim < Nstates; s_prim++) {
712 if (output_parity(s_prim , 0)) { exp_temp0 = -rp; }
713 else { exp_temp0 = rp; }
714 if (output_parity(s_prim , 1)) { exp_temp1 = -rp; }
715 else { exp_temp1 = rp; }
716 nom = llrcalc.
jaclog(nom, alpha_q(s_prim, kk) + exp_temp0 + beta_q(state_trans(s_prim, 0), k));
717 den = llrcalc.
jaclog(den, alpha_q(s_prim, kk) + exp_temp1 + beta_q(state_trans(s_prim, 1), k));
719 extrinsic_output(kk) = nom - den;
726 llrcalc = in_llrcalc;