50 result(0) = result(1) - l;
54 result(1) = result(0) + l;
62 for (
int i = 0; i <
length(l); i++) {
69 QLLR scaled_norm,
int j, QLLRvec &p1,
72 QLLR log_apriori_prob_const_point = 0;
74 for (
int i = 0; i <
k(j); i++) {
75 log_apriori_prob_const_point +=
76 ((
bitmap(j)(s, i) == 0) ? logP_apriori(b)(1) : logP_apriori(b)(0));
81 for (
int i = 0; i <
k(j); i++) {
82 if (
bitmap(j)(s, i) == 0) {
84 + log_apriori_prob_const_point);
88 + log_apriori_prob_const_point);
95 const ivec &s, QLLR scaled_norm,
96 QLLRvec &p1, QLLRvec &p0)
98 QLLR log_apriori_prob_const_point = 0;
100 for (
int j = 0; j <
nt; j++) {
101 for (
int i = 0; i <
k(j); i++) {
102 log_apriori_prob_const_point +=
103 ((
bitmap(j)(s[j], i) == 0) ? logP_apriori(b)(1) : logP_apriori(b)(0));
109 for (
int j = 0; j <
nt; j++) {
110 for (
int i = 0; i <
k(j); i++) {
111 if (
bitmap(j)(s[j], i) == 0) {
113 + log_apriori_prob_const_point);
117 + log_apriori_prob_const_point);
131 "The number of input bits does not match.");
133 out_symbols.set_size(
nt);
136 for (
int i = 0; i <
nt; ++i) {
137 int symb =
bin2dec(bits.mid(b,
k(i)));
153 const QLLRvec &LLR_apriori,
154 QLLRvec &LLR_aposteriori,
162 it_assert(H.rows() >= H.cols(),
"Modulator_NRD::demodulate_soft_bits():"
163 " ZF demodulation impossible for undetermined systems");
166 mat inv_HtH =
inv(Ht * H);
167 vec shat = inv_HtH * Ht * y;
168 vec h =
ones(shat.size());
169 for (
int i = 0; i < shat.size(); ++i) {
171 double sigma_zf =
std::sqrt(inv_HtH(i, i) * sigma2);
179 it_error(
"Modulator_NRD::demodulate_soft_bits(): Improper soft "
180 "demodulation method");
186 const QLLRvec &LLR_apriori,
196 const QLLRvec &LLR_apriori,
197 QLLRvec &LLR_aposteriori)
200 "Modulator_NRD::demodulate_soft_bits(): Wrong sizes");
202 "Modulator_NRD::demodulate_soft_bits(): Wrong sizes");
205 LLR_aposteriori.set_size(LLR_apriori.size());
208 double moo2s2 = -1.0 / (2.0 * sigma2);
211 for (
int i = 0; i <
nt; ++i) {
212 QLLRvec bnum = -QLLR_MAX *
ones_i(
k(i));
213 QLLRvec bdenom = bnum;
215 for (
int j = 0; j <
M(i); ++j) {
216 double norm2 = moo2s2 *
sqr(y(i) - h(i) *
symbols(i)(j));
218 update_LLR(logP_apriori, j, scaled_norm, i, bnum, bdenom);
220 LLR_aposteriori.set_subvector(b, bnum - bdenom);
227 const QLLRvec &LLR_apriori,
228 QLLRvec &LLR_aposteriori)
233 "Modulator_NRD::demodulate_soft_bits(): Wrong sizes");
235 "Modulator_NRD::demodulate_soft_bits(): Wrong sizes");
238 LLR_aposteriori.set_size(LLR_apriori.size());
241 double moo2s2 = -1.0 / (2.0 * sigma2);
243 bool diff_update =
false;
244 for (
int i = 0; i <
length(
M); ++i) {
258 QLLRvec bnum = -QLLR_MAX *
ones_i(np);
259 QLLRvec bdenom = bnum;
273 if (s(r) >
M(r) - 1) {
283 for (
int p = 0; p < nr; ++p) {
285 for (
int i = 0; i <
nt; ++i) {
286 d -= H(p, i) *
symbols(i)[s[i]];
292 update_LLR(logP_apriori, s, scaled_norm, bnum, bdenom);
307 LLR_aposteriori = bnum - bdenom;
313 const vec &ytH,
const mat &HtH,
320 norm +=
sqr(cdiff) * HtH(k, k);
322 norm -= cdiff * ytH[
k];
323 for (
int i = 0; i < m; ++i) {
324 norm += cdiff * HtH(i, k) *
symbols(k)[s[i]];
331 os <<
"--- REAL MIMO (NRD) CHANNEL ---------" << std::endl;
332 os <<
"Dimension (nt): " << mod.
nt << std::endl;
333 os <<
"Bits per dimension (k): " << mod.
k << std::endl;
334 os <<
"Symbols per dimension (M):" << mod.
M << std::endl;
335 for (
int i = 0; i < mod.
nt; i++) {
336 os <<
"Bitmap for dimension " << i <<
": " << mod.
bitmap(i) << std::endl;
338 os <<
"Symbol coordinates for dimension " << i <<
": " << mod.
symbols(i).
left(mod.
M(i)) << std::endl;
351 "The number of input bits does not match.");
353 out_symbols.set_size(
nt);
356 for (
int i = 0; i <
nt; ++i) {
357 int symb =
bin2dec(bits.mid(b,
k(i)));
373 const QLLRvec &LLR_apriori,
374 QLLRvec &LLR_aposteriori,
382 it_assert(H.rows() >= H.cols(),
"Modulator_NCD::demodulate_soft_bits():"
383 " ZF demodulation impossible for undetermined systems");
386 cmat inv_HhtH =
inv(Hht * H);
387 cvec shat = inv_HhtH * Hht * y;
388 cvec h =
ones_c(shat.size());
389 for (
int i = 0; i < shat.size(); ++i) {
398 it_error(
"Modulator_NCD::demodulate_soft_bits(): Improper soft "
399 "demodulation method");
405 const QLLRvec &LLR_apriori,
416 const QLLRvec &LLR_apriori,
417 QLLRvec &LLR_aposteriori)
420 "Modulator_NCD::demodulate_soft_bits(): Wrong sizes");
422 "Modulator_NCD::demodulate_soft_bits(): Wrong sizes");
425 LLR_aposteriori.set_size(LLR_apriori.size());
428 double moos2 = -1.0 / sigma2;
431 for (
int i = 0; i <
nt; ++i) {
432 QLLRvec bnum = -QLLR_MAX *
ones_i(
k(i));
433 QLLRvec bdenom = -QLLR_MAX *
ones_i(
k(i));
435 for (
int j = 0; j <
M(i); ++j) {
436 double norm2 = moos2 *
sqr(y(i) - h(i) *
symbols(i)(j));
438 update_LLR(logP_apriori, j, scaled_norm, i, bnum, bdenom);
440 LLR_aposteriori.set_subvector(b, bnum - bdenom);
447 const QLLRvec &LLR_apriori,
448 QLLRvec &LLR_aposteriori)
453 "Modulator_NRD::demodulate_soft_bits(): Wrong sizes");
455 "Modulator_NRD::demodulate_soft_bits(): Wrong sizes");
458 LLR_aposteriori.set_size(LLR_apriori.size());
461 double moos2 = -1.0 / sigma2;
463 bool diff_update =
false;
464 for (
int i = 0; i <
length(
M); ++i) {
473 cmat HtH = H.hermitian_transpose() * H;
474 cvec ytH =
conj(H.hermitian_transpose() * y);
476 QLLRvec bnum = -QLLR_MAX *
ones_i(np);
477 QLLRvec bdenom = -QLLR_MAX *
ones_i(np);
492 if (s(r) >
M(r) - 1) {
502 for (
int p = 0; p < nr; ++p) {
503 std::complex<double> d = y(p);
504 for (
int i = 0; i <
nt; ++i) {
505 d -= H(p, i) *
symbols(i)[s[i]];
511 update_LLR(logP_apriori, s, scaled_norm, bnum, bdenom);
526 LLR_aposteriori = bnum - bdenom;
531 const cvec &ytH,
const cmat &HtH,
537 norm +=
sqr(cdiff) * (HtH(k, k).real());
539 norm -= (cdiff.real() * ytH[
k].real() - cdiff.imag() * ytH[
k].imag());
540 for (
int i = 0; i < m; i++) {
548 os <<
"--- COMPLEX MIMO (NCD) CHANNEL --------" << std::endl;
549 os <<
"Dimension (nt): " << mod.
nt << std::endl;
550 os <<
"Bits per dimension (k): " << mod.
k << std::endl;
551 os <<
"Symbols per dimension (M):" << mod.
M << std::endl;
552 for (
int i = 0; i < mod.
nt; i++) {
553 os <<
"Bitmap for dimension " << i <<
": "
554 << mod.
bitmap(i) << std::endl;
555 os <<
"Symbol coordinates for dimension " << i <<
": "
588 spacing.set_size(
nt);
590 for (
int i = 0; i <
nt; i++) {
593 "ND_UPAM::set_M(): M is not a power of 2.");
598 double average_energy = (
M(i) *
M(i) - 1) / 3.0;
599 double scaling_factor =
std::sqrt(average_energy);
601 for (
int j = 0; j <
M(i); ++j) {
602 symbols(i)(j) = ((
M(i) - 1) - j * 2) / scaling_factor;
610 spacing(i) = 2.0 / scaling_factor;
614 int ND_UPAM::sphere_search_SE(
const vec &y_in,
const mat &H,
615 const imat &zrange,
double r, ivec &zhat)
626 mat R =
chol(H.transpose() * H);
629 vec y = Q.transpose() * y_in;
630 mat Vi = Ri.transpose();
635 double bestdist = r * r;
639 for (
int i = 0; i < n; i++) {
640 for (
int j = 0; j < n; j++) {
641 E(i*n + n - 1) += y(j) * Vi(j + n * i);
647 z(n - 1) =
floor_i(0.5 + E(n * n - 1));
648 z(n - 1) =
std::max(z(n - 1), zrange(n - 1, 0));
649 z(n - 1) =
std::min(z(n - 1), zrange(n - 1, 1));
650 double p = (E(n * n - 1) - z(n - 1)) / Vi(n * n - 1);
652 step(n - 1) = sign_nozero_i(p);
658 double newdist = dist(k) + p * p;
660 if ((newdist < bestdist) && (k != 0)) {
661 for (
int i = 0; i <
k; i++) {
662 E(k - 1 + i*n) = E(k + i * n) - p * Vi(k + i * n);
667 z(k) =
floor_i(0.5 + E(k + k * n));
668 z(k) =
std::max(z(k), zrange(k, 0));
669 z(k) =
std::min(z(k), zrange(k, 1));
670 p = (E(k + k * n) - z(k)) / Vi(k + k * n);
672 step(k) = sign_nozero_i(p);
676 if (newdist < bestdist) {
681 else if (k == n - 1) {
690 if ((z(k) < zrange(k, 0)) || (z(k) > zrange(k, 1))) {
691 step(k) = (-step(k) - sign_nozero_i(step(k)));
695 if ((z(k) >= zrange(k, 0)) && (z(k) <= zrange(k, 1))) {
700 p = (E(k + k * n) - z(k)) / Vi(k + k * n);
701 step(k) = (-step(k) - sign_nozero_i(step(k)));
711 double rmax,
double stepup,
712 QLLRvec &detected_bits)
715 "ND_UPAM::sphere_decoding(): dimension mismatch");
717 "ND_UPAM::sphere_decoding(): dimension mismatch");
718 it_assert(rstart > 0,
"ND_UPAM::sphere_decoding(): radius error");
719 it_assert(rmax > rstart,
"ND_UPAM::sphere_decoding(): radius error");
724 mat Htemp(H.rows(), H.cols());
725 for (
int i = 0; i < H.cols(); i++) {
726 Htemp.set_col(i, H.get_col(i)*spacing(i));
727 ytemp += Htemp.get_col(i) * 0.5 * (
M(i) - 1.0);
731 for (
int i = 0; i <
nt; i++) {
733 crange(i, 1) =
M(i) - 1;
740 status = sphere_search_SE(ytemp, Htemp, crange, r, s);
743 detected_bits.set_size(
sum(k));
745 for (
int j = 0; j <
nt; j++) {
746 for (
int i = 0; i <
k(j); i++) {
747 if (
bitmap(j)((
M(j) - 1 - s[j]), i) == 0) {
748 detected_bits(b) = 1000;
751 detected_bits(b) = -1000;
798 for (
int i = 0; i <
nt; ++i) {
801 "ND_UQAM::set_M(): M is not a power of 2");
804 it_assert(
L(i)*
L(i) ==
M(i),
"ND_UQAM: constellation M must be square");
809 double average_energy = (
M(i) - 1) * 2.0 / 3.0;
810 double scaling_factor =
std::sqrt(average_energy);
813 for (
int j1 = 0; j1 <
L(i); ++j1) {
814 for (
int j2 = 0; j2 <
L(i); ++j2) {
816 std::complex<double>(((
L(i) - 1) - j2 * 2.0) / scaling_factor,
817 ((
L(i) - 1) - j1 * 2.0) / scaling_factor);
818 bitmap(i).set_row(j1*
L(i) + j2,
concat(gray_code.get_row(j1),
819 gray_code.get_row(j2)));
858 for (
int i = 0; i <
nt; ++i) {
861 "ND_UPSK::set_M(): M is not a power of 2");
867 double delta = 2.0 *
pi /
M(i);
868 double epsilon = delta / 10000.0;
870 for (
int j = 0; j <
M(i); ++j) {
871 std::complex<double> symb
872 = std::complex<double>(std::polar(1.0, delta * j));