59 result(0) = result(1) - l;
63 result(1) = result(0) + l;
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);
133 QLLR logsum0, logsum1;
134 const QLLR *
const addrfirst =
Qnorms._data();
135 const QLLR *
const addrsemilast = addrfirst + (1 << (
nb - 1)), *
const addrlast = addrfirst + (1 <<
nb);
137 for(
int bi = 3; bi <
nb - 1 ; bi++) {
140 const int forhalfdiff = 1 << bi, fordiff = 2 * forhalfdiff, fordbldiff = 2 * fordiff;
142 const QLLR *
const addr1 = addrfirst + forhalfdiff, *
const addr2 = addr1 + fordiff, *
const addr3 = addrlast - fordiff;
143 while(Qptr < addr1) logsum0 =
llrcalc.
jaclog(*(Qptr++), logsum0);
144 while(Qptr < addr2) logsum1 =
llrcalc.
jaclog(*(Qptr++), logsum1);
145 const QLLR *addrdyn0, *addrdyn1;
146 while(Qptr < addr3) {
147 addrdyn0 = Qptr + fordiff;
148 addrdyn1 = Qptr + fordbldiff;
149 while(Qptr < addrdyn0) logsum0 =
llrcalc.
jaclog(*(Qptr++), logsum0);
150 while(Qptr < addrdyn1) logsum1 =
llrcalc.
jaclog(*(Qptr++), logsum1);
152 while(Qptr < addrlast) logsum0 =
llrcalc.
jaclog(*(Qptr++), logsum0);
153 llr[bi] = logsum0 - logsum1;
159 while(Qptr < addrsemilast) logsum0 =
llrcalc.
jaclog(*(Qptr++), logsum0);
160 while(Qptr < addrlast) logsum1 =
llrcalc.
jaclog(*(Qptr++), logsum1);
161 llr[
nb - 1] = logsum0 - logsum1;
170 QLLR logmax0, logmax1;
171 const QLLR *
const addrfirst =
Qnorms._data();
172 const QLLR *
const addrsemilast = addrfirst + (1 << (
nb - 1)), *
const addrlast = addrfirst + (1 <<
nb);
174 for(
int bi = 3; bi <
nb - 1; bi++) {
177 const int forhalfdiff = 1 << bi, fordiff = 2 * forhalfdiff, fordbldiff = 2 * fordiff;
179 const QLLR *
const addr1 = addrfirst + forhalfdiff, *
const addr2 = addr1 + fordiff, *
const addr3 = addrlast - fordiff;
180 for(; Qptr < addr1; Qptr++) logmax0 = *Qptr > logmax0 ? *Qptr : logmax0;
181 for(; Qptr < addr2; Qptr++) logmax1 = *Qptr > logmax1 ? *Qptr : logmax1;
182 const QLLR *addrdyn0, *addrdyn1;
183 while(Qptr < addr3) {
184 addrdyn0 = Qptr + fordiff;
185 addrdyn1 = Qptr + fordbldiff;
186 for(; Qptr < addrdyn0; Qptr++) logmax0 = *Qptr > logmax0 ? *Qptr : logmax0;
187 for(; Qptr < addrdyn1; Qptr++) logmax1 = *Qptr > logmax1 ? *Qptr : logmax1;
189 for(; Qptr < addrlast; Qptr++) logmax0 = *Qptr > logmax0 ? *Qptr : logmax0;
190 llr[bi] = logmax0 - logmax1;
196 for(; Qptr < addrsemilast; Qptr++) logmax0 = *Qptr > logmax0 ? *Qptr : logmax0;
197 for(; Qptr < addrlast; Qptr++) logmax1 = *Qptr > logmax1 ? *Qptr : logmax1;
198 llr[
nb - 1] = logmax0 - logmax1;
200 else it_error(
"Improper soft demodulation method\n.");
205 using namespace itpp;
206 QLLR logsum0 = -QLLR_MAX, logsum1 = -QLLR_MAX;
207 const QLLR *
const addrfirst =
Qnorms._data(), *
const addr3 = addrfirst + (1 <<
nb) - 1;
208 const QLLR *Qptr = addrfirst;
212 while(Qptr < addr3) {
219 llr = logsum0 - logsum1;
224 using namespace itpp;
225 QLLR logsum0 = -QLLR_MAX, logsum1 = -QLLR_MAX;
226 const QLLR *
const addrfirst =
Qnorms._data(), *
const addr3 = addrfirst + (1 <<
nb) - 2;
227 const QLLR *Qptr = addrfirst;
235 while(Qptr < addr3) {
247 llr = logsum0 - logsum1;
252 using namespace itpp;
253 QLLR logsum0 = -QLLR_MAX, logsum1 = -QLLR_MAX;
254 const QLLR *
const addrfirst =
Qnorms._data(), *
const addr3 = addrfirst + (1 <<
nb) - 4;
255 const QLLR *Qptr = addrfirst;
269 while(Qptr < addr3) {
291 llr = logsum0 - logsum1;
296 using namespace itpp;
297 QLLR logmax0 = -QLLR_MAX, logmax1 = -QLLR_MAX;
298 const QLLR *
const addrfirst =
Qnorms._data(), *
const addr3 = addrfirst + (1 <<
nb) - 1;
299 const QLLR *Qptr = addrfirst;
300 logmax0 = *Qptr > logmax0 ? *Qptr : logmax0;
302 logmax1 = *Qptr > logmax1 ? *Qptr : logmax1;
304 logmax1 = *Qptr > logmax1 ? *Qptr : logmax1;
306 while(Qptr < addr3) {
307 logmax0 = *Qptr > logmax0 ? *Qptr : logmax0;
309 logmax0 = *Qptr > logmax0 ? *Qptr : logmax0;
311 logmax1 = *Qptr > logmax1 ? *Qptr : logmax1;
313 logmax1 = *Qptr > logmax1 ? *Qptr : logmax1;
316 logmax0 = *Qptr > logmax0 ? *Qptr : logmax0;
317 maxllr = logmax0 - logmax1;
322 using namespace itpp;
323 QLLR logmax0 = -QLLR_MAX, logmax1 = -QLLR_MAX;
324 const QLLR *
const addrfirst =
Qnorms._data(), *
const addr3 = addrfirst + (1 <<
nb) - 2;
325 const QLLR *Qptr = addrfirst;
326 logmax0 = *Qptr > logmax0 ? *Qptr : logmax0;
328 logmax0 = *Qptr > logmax0 ? *Qptr : logmax0;
330 logmax1 = *Qptr > logmax1 ? *Qptr : logmax1;
332 logmax1 = *Qptr > logmax1 ? *Qptr : logmax1;
334 logmax1 = *Qptr > logmax1 ? *Qptr : logmax1;
336 logmax1 = *Qptr > logmax1 ? *Qptr : logmax1;
338 while(Qptr < addr3) {
339 logmax0 = *Qptr > logmax0 ? *Qptr : logmax0;
341 logmax0 = *Qptr > logmax0 ? *Qptr : logmax0;
343 logmax0 = *Qptr > logmax0 ? *Qptr : logmax0;
345 logmax0 = *Qptr > logmax0 ? *Qptr : logmax0;
347 logmax1 = *Qptr > logmax1 ? *Qptr : logmax1;
349 logmax1 = *Qptr > logmax1 ? *Qptr : logmax1;
351 logmax1 = *Qptr > logmax1 ? *Qptr : logmax1;
353 logmax1 = *Qptr > logmax1 ? *Qptr : logmax1;
356 logmax0 = *Qptr > logmax0 ? *Qptr : logmax0;
358 logmax0 = *Qptr > logmax0 ? *Qptr : logmax0;
359 maxllr = logmax0 - logmax1;
364 using namespace itpp;
365 QLLR logmax0 = -QLLR_MAX, logmax1 = -QLLR_MAX;
366 const QLLR *
const addrfirst =
Qnorms._data(), *
const addr3 = addrfirst + (1 <<
nb) - 4;
367 const QLLR *Qptr = addrfirst;
368 logmax0 = *Qptr > logmax0 ? *Qptr : logmax0;
370 logmax0 = *Qptr > logmax0 ? *Qptr : logmax0;
372 logmax0 = *Qptr > logmax0 ? *Qptr : logmax0;
374 logmax0 = *Qptr > logmax0 ? *Qptr : logmax0;
376 logmax1 = *Qptr > logmax1 ? *Qptr : logmax1;
378 logmax1 = *Qptr > logmax1 ? *Qptr : logmax1;
380 logmax1 = *Qptr > logmax1 ? *Qptr : logmax1;
382 logmax1 = *Qptr > logmax1 ? *Qptr : logmax1;
384 logmax1 = *Qptr > logmax1 ? *Qptr : logmax1;
386 logmax1 = *Qptr > logmax1 ? *Qptr : logmax1;
388 logmax1 = *Qptr > logmax1 ? *Qptr : logmax1;
390 logmax1 = *Qptr > logmax1 ? *Qptr : logmax1;
392 while(Qptr < addr3) {
393 logmax0 = *Qptr > logmax0 ? *Qptr : logmax0;
395 logmax0 = *Qptr > logmax0 ? *Qptr : logmax0;
397 logmax0 = *Qptr > logmax0 ? *Qptr : logmax0;
399 logmax0 = *Qptr > logmax0 ? *Qptr : logmax0;
401 logmax0 = *Qptr > logmax0 ? *Qptr : logmax0;
403 logmax0 = *Qptr > logmax0 ? *Qptr : logmax0;
405 logmax0 = *Qptr > logmax0 ? *Qptr : logmax0;
407 logmax0 = *Qptr > logmax0 ? *Qptr : logmax0;
409 logmax1 = *Qptr > logmax1 ? *Qptr : logmax1;
411 logmax1 = *Qptr > logmax1 ? *Qptr : logmax1;
413 logmax1 = *Qptr > logmax1 ? *Qptr : logmax1;
415 logmax1 = *Qptr > logmax1 ? *Qptr : logmax1;
417 logmax1 = *Qptr > logmax1 ? *Qptr : logmax1;
419 logmax1 = *Qptr > logmax1 ? *Qptr : logmax1;
421 logmax1 = *Qptr > logmax1 ? *Qptr : logmax1;
423 logmax1 = *Qptr > logmax1 ? *Qptr : logmax1;
426 logmax0 = *Qptr > logmax0 ? *Qptr : logmax0;
428 logmax0 = *Qptr > logmax0 ? *Qptr : logmax0;
430 logmax0 = *Qptr > logmax0 ? *Qptr : logmax0;
432 logmax0 = *Qptr > logmax0 ? *Qptr : logmax0;
433 maxllr = logmax0 - logmax1;
440 for(
int i = 0; i <
length(l); i++) {
455 for(
int i = 0; i <
nt; ++i) {
456 it_assert(
M.length() == symbols.length(),
"Modulator_NRD::get_symbols(): " 457 "The length of M vector is different than length of the symbols vector.");
458 retvec(i) = symbols(i).
left(
M(i));
466 "The number of input bits does not match.");
468 out_symbols.set_size(
nt);
471 for(
int i = 0; i <
nt; ++i) {
472 int symb =
bin2dec(bits.mid(b,
k(i)));
481 modulate_bits(bits, result);
488 using namespace itpp;
489 it_assert(H_in.cols() ==
nt,
"Number of Tx antennas is wrong.\n");
490 it_assert(
sum(
k) < 32,
"Number of total bits per transmission can not be larger than 32.\n");
491 it_assert(
pow2i(
sum(
k)) ==
prod(
M),
"Modulator must use exhaustive constellations, i.e., #bits=log2(#symbs).\n");
497 hspacings.set_size(
nt);
498 yspacings.set_size(
nt);
502 vec startsymbvec(
nt);
503 for(
int ci = 0; ci <
nt; ci++) startsymbvec[ci] = symbols(ci)[0];
504 itpp::vec Hx = H * startsymbvec;
505 for(
int ci = 0, bcs = 0; ci <
nt; bcs +=
k[ci++]) {
506 for(
int bi = 0; bi <
k[ci]; bi++)
bpos2cpos[bcs + bi] = ci;
508 for(
int si = 0; si < M[ci]; si++) gray2dec(ci)[si ^(si >> 1)] = si;
509 yspacings(ci).set_size(
M[ci] - 1);
510 hspacings(ci).set_size(
M[ci] - 1);
511 for(
int si = 0; si <
M[ci] - 1; si++) {
512 double xspacing = symbols(ci)[
bits2symbols(ci)[(si + 1) ^((si + 1) >> 1)]];
513 xspacing -= symbols(ci)[
bits2symbols(ci)[si ^(si >> 1)]];
514 hspacings(ci)(si) = H.get_col(ci) * xspacing;
518 unsigned bitstring = 0, ind = 0;
519 hxnormupdate(Hx, bitstring, ind,
nb - 1);
526 const itpp::QLLRvec& llr_apr,
530 using namespace itpp;
533 it_assert_debug(H.rows() == y.length(),
"The dimensions are not correct.\n");
542 vec ytil = H.T() * y;
543 vec startsymbvec(
nt);
544 for(
int ci = 0; ci <
nt; ci++) startsymbvec[ci] = symbols(ci)[0];
545 double yx = 2*(ytil * startsymbvec);
549 for(
int ci = 0; ci <
nt; ci++)
for(
int si = 0; si <
M[ci] - 1; si++) {
550 double xspacing = symbols(ci)[
bits2symbols(ci)[(si + 1) ^((si + 1) >> 1)]];
551 xspacing -= symbols(ci)[
bits2symbols(ci)[si ^(si >> 1)]];
552 yspacings(ci)[si] = 2*(ytil(ci) * xspacing);
554 unsigned bitstring = 0, ind = 0;
555 yxnormupdate(yx, lapr, bitstring, ind,
nb - 1);
562 const QLLRvec &LLR_apriori,
563 QLLRvec &LLR_aposteriori,
568 it_assert(H.rows() >= H.cols(),
"Modulator_NRD::demodulate_soft_bits():" 569 " ZF demodulation impossible for undetermined systems");
572 mat inv_HtH =
inv(Ht * H);
573 vec shat = inv_HtH * Ht * y;
574 vec h =
ones(shat.size());
575 for(
int i = 0; i < shat.size(); ++i) {
577 double sigma_zf =
std::sqrt(inv_HtH(i, i) * sigma2);
581 demodulate_soft_bits(shat, h, 1.0,
zeros_i(
sum(
k)), LLR_aposteriori);
585 init_soft_demodulator(H, sigma2);
586 demodulate_soft_bits(y, LLR_apriori, LLR_aposteriori, method);
593 const QLLRvec &LLR_apriori,
597 demodulate_soft_bits(y, H, sigma2, LLR_apriori, result, method);
603 const QLLRvec &LLR_apriori,
604 QLLRvec &LLR_aposteriori)
607 "Modulator_NRD::demodulate_soft_bits(): Wrong sizes");
609 "Modulator_NRD::demodulate_soft_bits(): Wrong sizes");
612 LLR_aposteriori.set_size(LLR_apriori.size());
615 double moo2s2 = -1.0 / (2.0 * sigma2);
618 for(
int i = 0; i <
nt; ++i) {
619 QLLRvec bnum = -QLLR_MAX *
ones_i(
k(i));
620 QLLRvec bdenom = bnum;
622 for(
int j = 0; j <
M(i); ++j) {
623 double norm2 = moo2s2 *
sqr(y(i) - h(i) * symbols(i)(j));
625 update_LLR(logP_apriori, j, scaled_norm, i, bnum, bdenom);
627 LLR_aposteriori.set_subvector(b, bnum - bdenom);
634 using namespace itpp;
638 unsigned oldi =
gray2dec(col)[bitstring & (
M[col] - 1)];
640 unsigned newi =
gray2dec(col)[bitstring & (
M[col] - 1)];
641 Hx += oldi > newi ? -hspacings(col)(newi) : hspacings(col)(oldi);
645 hxnormupdate(Hx, bitstring, ind, bit - 1);
647 bitstring ^= 1 << bit;
649 Hx += oldi > newi ? -hspacings(col)(newi) : hspacings(col)(oldi);
650 hxnormupdate(Hx, bitstring, ind, bit - 1);
655 using namespace itpp;
660 unsigned oldi =
gray2dec(col)[bitstring & (
M[col] - 1)];
662 unsigned newi =
gray2dec(col)[bitstring & (
M[col] - 1)];
663 yx += oldi > newi ? -yspacings(col)[newi] : yspacings(col)[oldi];
669 yxnormupdate(yx, lapr, bitstring, ind, bit - 1);
671 bitstring ^= 1 << bit;
673 yx += oldi > newi ? -yspacings(col)[newi] : yspacings(col)[oldi];
674 lapr += ((bitstring >> bit) & 1) ? -
llrapr[bit] :
llrapr[bit];
675 yxnormupdate(yx, lapr, bitstring, ind, bit - 1);
681 os <<
"--- REAL MIMO (NRD) CHANNEL ---------" << std::endl;
682 os <<
"Dimension (nt): " <<
mod.nt << std::endl;
683 os <<
"Bits per dimension (k): " <<
mod.k << std::endl;
684 os <<
"Symbols per dimension (M):" <<
mod.M << std::endl;
685 for(
int i = 0; i <
mod.nt; i++) {
686 os <<
"Bitmap for dimension " << i <<
": " <<
mod.bitmap(i) << std::endl;
688 os <<
"Symbol coordinates for dimension " << i <<
": " <<
mod.symbols(i).left(
mod.M(i)) << std::endl;
690 os <<
mod.get_llrcalc() << std::endl;
701 for(
int i = 0; i <
nt; ++i) {
702 it_assert(
M.length() == symbols.length(),
"Modulator_NRD::get_symbols(): " 703 "The length of M vector is different than length of the symbols vector.");
704 retvec(i) = symbols(i).
left(
M(i));
712 "The number of input bits does not match.");
714 out_symbols.set_size(
nt);
717 for(
int i = 0; i <
nt; ++i) {
718 int symb =
bin2dec(bits.mid(b,
k(i)));
727 modulate_bits(bits, result);
733 using namespace itpp;
735 it_assert_debug(
sum(
k) < 32,
"Number of total bits per transmission can not be larger than 32.\n");
742 hspacings.set_size(
nt);
743 yspacings.set_size(
nt);
747 cvec startsymbvec(
nt);
748 for(
int ci = 0; ci <
nt; ci++) startsymbvec[ci] = symbols(ci)[0];
749 cvec Hx = H * startsymbvec;
750 for(
int ci = 0, bcs = 0; ci <
nt; bcs +=
k[ci++]) {
751 for(
int bi = 0; bi <
k[ci]; bi++)
bpos2cpos[bcs + bi] = ci;
753 for(
int si = 0; si < M[ci]; si++) gray2dec(ci)[si ^(si >> 1)] = si;
754 yspacings(ci).set_size(
M[ci] - 1);
755 hspacings(ci).set_size(
M[ci] - 1);
756 for(
int si = 0; si <
M[ci] - 1; si++) {
757 std::complex<double> xspacing = symbols(ci)[
bits2symbols(ci)[(si + 1) ^((si + 1) >> 1)]];
758 xspacing -= symbols(ci)[
bits2symbols(ci)[si ^(si >> 1)]];
759 hspacings(ci)(si) = H.get_col(ci) * xspacing;
763 unsigned bitstring = 0, ind = 0;
764 hxnormupdate(Hx, bitstring, ind,
nb - 1);
769 const itpp::QLLRvec& llr_apr,
773 using namespace itpp;
776 it_assert_debug(H.rows() == y.length(),
"The dimensions are not correct.\n");
784 cvec ytil =
conj(H.H() * y);
785 cvec startsymbvec(
nt);
786 for(
int ci = 0; ci <
nt; ci++) startsymbvec[ci] = symbols(ci)[0];
787 double yx = 2*(ytil * startsymbvec).
real();
790 for(
int ci = 0; ci <
nt; ci++)
for(
int si = 0; si <
M[ci] - 1; si++) {
791 std::complex<double> xspacing = symbols(ci)[
bits2symbols(ci)[(si + 1) ^((si + 1) >> 1)]];
792 xspacing -= symbols(ci)[
bits2symbols(ci)[si ^(si >> 1)]];
793 yspacings(ci)[si] = 2*(ytil[ci] * xspacing).
real();
795 unsigned bitstring = 0, ind = 0;
796 yxnormupdate(yx, lapr, bitstring, ind,
nb - 1);
803 const QLLRvec &LLR_apriori,
804 QLLRvec &LLR_aposteriori,
809 it_assert(H.rows() >= H.cols(),
"Modulator_NCD::demodulate_soft_bits():" 810 " ZF demodulation impossible for undetermined systems");
813 cmat inv_HhtH =
inv(Hht * H);
814 cvec shat = inv_HhtH * Hht * y;
815 cvec h =
ones_c(shat.size());
816 for(
int i = 0; i < shat.size(); ++i) {
821 demodulate_soft_bits(shat, h, 1.0,
zeros_i(
sum(
k)), LLR_aposteriori);
825 init_soft_demodulator(H, sigma2);
826 demodulate_soft_bits(y, LLR_apriori, LLR_aposteriori, method);
833 const QLLRvec &LLR_apriori,
837 demodulate_soft_bits(y, H, sigma2, LLR_apriori, result, method);
844 const QLLRvec &LLR_apriori,
845 QLLRvec &LLR_aposteriori)
848 "Modulator_NCD::demodulate_soft_bits(): Wrong sizes");
850 "Modulator_NCD::demodulate_soft_bits(): Wrong sizes");
853 LLR_aposteriori.set_size(LLR_apriori.size());
856 double moos2 = -1.0 / sigma2;
859 for(
int i = 0; i <
nt; ++i) {
860 QLLRvec bnum = -QLLR_MAX *
ones_i(
k(i));
861 QLLRvec bdenom = -QLLR_MAX *
ones_i(
k(i));
863 for(
int j = 0; j <
M(i); ++j) {
864 double norm2 = moos2 *
sqr(y(i) - h(i) * symbols(i)(j));
866 update_LLR(logP_apriori, j, scaled_norm, i, bnum, bdenom);
868 LLR_aposteriori.set_subvector(b, bnum - bdenom);
874 void Modulator_NCD::hxnormupdate(itpp::cvec& Hx,
unsigned& bitstring,
unsigned& ind,
unsigned bit)
876 using namespace itpp;
880 unsigned oldi =
gray2dec(col)[bitstring & (
M[col] - 1)];
882 unsigned newi =
gray2dec(col)[bitstring & (
M[col] - 1)];
883 Hx += oldi > newi ? -hspacings(col)(newi) : hspacings(col)(oldi);
887 hxnormupdate(Hx, bitstring, ind, bit - 1);
889 bitstring ^= 1 << bit;
891 Hx += oldi > newi ? -hspacings(col)(newi) : hspacings(col)(oldi);
892 hxnormupdate(Hx, bitstring, ind, bit - 1);
895 void Modulator_NCD::yxnormupdate(
double& yx, itpp::QLLR& lapr,
unsigned& bitstring,
unsigned& ind,
unsigned bit)
897 using namespace itpp;
904 unsigned oldi =
gray2dec(col)[bitstring & (
M[col] - 1)];
906 unsigned newi =
gray2dec(col)[bitstring & (
M[col] - 1)];
907 yx += oldi > newi ? -yspacings(col)[newi] : yspacings(col)[oldi];
913 yxnormupdate(yx, lapr, bitstring, ind, bit - 1);
915 bitstring ^= 1 << bit;
917 yx += oldi > newi ? -yspacings(col)[newi] : yspacings(col)[oldi];
918 lapr += ((bitstring >> bit) & 1) ? -
llrapr[bit] :
llrapr[bit];
919 yxnormupdate(yx, lapr, bitstring, ind, bit - 1);
925 os <<
"--- COMPLEX MIMO (NCD) CHANNEL --------" << std::endl;
926 os <<
"Dimension (nt): " <<
mod.nt << std::endl;
927 os <<
"Bits per dimension (k): " <<
mod.k << std::endl;
928 os <<
"Symbols per dimension (M):" <<
mod.M << std::endl;
929 for(
int i = 0; i <
mod.nt; i++) {
930 os <<
"Bitmap for dimension " << i <<
": " 931 <<
mod.bitmap(i) << std::endl;
932 os <<
"Symbol coordinates for dimension " << i <<
": " 933 <<
mod.symbols(i).left(
mod.M(i)) << std::endl;
935 os <<
mod.get_llrcalc() << std::endl;
953 set_M(
nt, Mary_temp);
963 symbols.set_size(
nt);
965 spacing.set_size(
nt);
967 for(
int i = 0; i <
nt; i++) {
970 "ND_UPAM::set_M(): M is not a power of 2.");
972 symbols(i).set_size(
M(i) + 1);
975 double average_energy = (
M(i) *
M(i) - 1) / 3.0;
976 double scaling_factor =
std::sqrt(average_energy);
978 for(
int j = 0; j <
M(i); ++j) {
979 symbols(i)(j) = ((
M(i) - 1) - j * 2) / scaling_factor;
985 symbols(i)(
M(i)) = 0.0;
987 spacing(i) = 2.0 / scaling_factor;
991 int ND_UPAM::sphere_search_SE(
const vec &y_in,
const mat &H,
992 const imat &zrange,
double r, ivec &zhat)
1003 mat R =
chol(H.transpose() * H);
1006 vec y = Q.transpose() * y_in;
1007 mat Vi = Ri.transpose();
1012 double bestdist = r * r;
1015 mat E =
zeros(n, n);
1016 for(
int i = 0; i < n; i++) {
1017 for(
int j = 0; j < n; j++) {
1018 E(i * n + n - 1) += y(j) * Vi(j + n * i);
1024 z(n - 1) =
floor_i(0.5 + E(n * n - 1));
1025 z(n - 1) =
std::max(z(n - 1), zrange(n - 1, 0));
1026 z(n - 1) =
std::min(z(n - 1), zrange(n - 1, 1));
1027 double p = (E(n * n - 1) - z(n - 1)) / Vi(n * n - 1);
1029 step(n - 1) = sign_nozero_i(p);
1035 double newdist = dist(
k) + p * p;
1037 if((newdist < bestdist) && (
k != 0)) {
1038 for(
int i = 0; i <
k; i++) {
1039 E(
k - 1 + i * n) = E(
k + i * n) - p * Vi(
k + i * n);
1047 p = (E(
k +
k * n) - z(
k)) / Vi(
k +
k * n);
1049 step(
k) = sign_nozero_i(p);
1053 if(newdist < bestdist) {
1058 else if(
k == n - 1) {
1067 if((z(
k) < zrange(
k, 0)) || (z(
k) > zrange(
k, 1))) {
1068 step(
k) = (-step(
k) - sign_nozero_i(step(
k)));
1072 if((z(
k) >= zrange(
k, 0)) && (z(
k) <= zrange(
k, 1))) {
1077 p = (E(
k +
k * n) - z(
k)) / Vi(
k +
k * n);
1078 step(
k) = (-step(
k) - sign_nozero_i(step(
k)));
1088 double rmax,
double stepup,
1089 QLLRvec &detected_bits)
1092 "ND_UPAM::sphere_decoding(): dimension mismatch");
1094 "ND_UPAM::sphere_decoding(): dimension mismatch");
1095 it_assert(rstart > 0,
"ND_UPAM::sphere_decoding(): radius error");
1096 it_assert(rmax > rstart,
"ND_UPAM::sphere_decoding(): radius error");
1101 mat Htemp(H.rows(), H.cols());
1102 for(
int i = 0; i < H.cols(); i++) {
1103 Htemp.set_col(i, H.get_col(i)*spacing(i));
1104 ytemp += Htemp.get_col(i) * 0.5 * (
M(i) - 1.0);
1108 for(
int i = 0; i <
nt; i++) {
1110 crange(i, 1) =
M(i) - 1;
1117 status = sphere_search_SE(ytemp, Htemp, crange, r, s);
1120 detected_bits.set_size(
sum(
k));
1122 for(
int j = 0; j <
nt; j++) {
1123 for(
int i = 0; i <
k(j); i++) {
1124 if(
bitmap(j)((
M(j) - 1 - s[j]), i) == 0) {
1125 detected_bits(b) = 1000;
1128 detected_bits(b) = -1000;
1161 set_M(
nt, Mary_temp);
1172 symbols.set_size(
nt);
1175 for(
int i = 0; i <
nt; ++i) {
1178 "ND_UQAM::set_M(): M is not a power of 2");
1181 it_assert(L(i)*L(i) ==
M(i),
"ND_UQAM: constellation M must be square");
1183 symbols(i).set_size(
M(i) + 1);
1186 double average_energy = (
M(i) - 1) * 2.0 / 3.0;
1187 double scaling_factor =
std::sqrt(average_energy);
1190 for(
int j1 = 0; j1 < L(i); ++j1) {
1191 for(
int j2 = 0; j2 < L(i); ++j2) {
1192 symbols(i)(j1 * L(i) + j2) =
1193 std::complex<double>(((L(i) - 1) - j2 * 2.0) / scaling_factor,
1194 ((L(i) - 1) - j1 * 2.0) / scaling_factor);
1204 symbols(i)(
M(i)) = 0.0;
1210 it_assert(
nt > nth,
"ND_UQAM::set_constellation_points(): Number of input to change is out of the size");
1211 it_assert(inConstellation.size() == in_bit2symbols.size(),
1212 "ND_UQAM::set_constellation_points(): Number of constellation and bits2symbols does not match");
1214 "ND_UQAM::set_constellation_points(): Number of symbols needs to be even and non-zero");
1216 symbols(nth).replace_mid(0, inConstellation);
1220 for(
int m = 0; m <
M(nth); ++m) {
1226 symbols(nth)(
M(nth)) = 0.0;
1243 set_M(
nt, Mary_temp);
1253 symbols.set_size(
nt);
1256 for(
int i = 0; i <
nt; ++i) {
1259 "ND_UPSK::set_M(): M is not a power of 2");
1261 symbols(i).set_size(
M(i) + 1);
1265 double delta = 2.0 *
pi /
M(i);
1266 double epsilon = delta / 10000.0;
1268 for(
int j = 0; j <
M(i); ++j) {
1269 std::complex<double> symb
1270 = std::complex<double>(std::polar(1.0, delta * j));
1273 symbols(i)(j) = std::complex<double>(0.0,
std::imag(symb));
1276 symbols(i)(j) = std::complex<double>(
std::real(symb), 0.0);
1279 symbols(i)(j) = symb;
1287 symbols(i)(
M(i)) = 0.0;
Various functions on vectors and matrices - header file.
ITPP_EXPORT int round_i(double x)
Round to nearest integer.
itpp::ivec bitcumsum
The cumulative sum of bits in the symbol vector.
ITPP_EXPORT ivec zeros_i(int size)
A Int vector of zeros.
ND_UPSK(int nt=1, int Mary=4)
Constructor.
Array< cvec > get_symbols() const
Get modulation symbols per dimension.
QLLRvec probabilities(QLLR l)
Convert LLR to log-probabilities.
ITPP_EXPORT ivec ones_i(int size)
A Int vector of ones.
Vec< T > reverse(const Vec< T > &in)
Reverse the input vector.
int nt
Number of dimensions.
vec imag(const cvec &data)
Imaginary part of complex values.
void yxnormupdate(double &yx, itpp::QLLR &lapr, unsigned &bitstring, unsigned &ind, unsigned bit)
Calculation of the remaining part of the norms that depends both on H and y.
Base class for N-dimensional vector (MIMO) channel modulators/demodulators with real-valued component...
LLR_calc_unit llrcalc
LLR calculation unit.
Array< ivec > bits2symbols
Bit pattern in decimal form ordered and the corresponding symbols (one pattern per dimension) ...
bool is_even(int x)
Return true if x is an even integer.
void set_M(int nt=1, int Mary=4)
Set component modulators to M-QAM with Gray mapping.
ND_UQAM(int nt=1, int Mary=4)
Constructor.
bool demod_initialized
Flag indicating whether the demodulator has been initialized.
double norm(const cvec &v)
Calculate the 2-norm: norm(v)=sqrt(sum(abs(v).^2))
cvec conj(const cvec &x)
Conjugate of complex value.
int mod(int k, int n)
Calculates the modulus, i.e. the signed reminder after division.
std::ostream & operator<<(std::ostream &output, const bin &inbin)
Output stream of bin.
QLLR to_qllr(double l) const
Convert a "real" LLR value to an LLR type.
void hxnormupdate(itpp::vec &Hx, unsigned &bitstring, unsigned &ind, unsigned bit)
Calculation of the part of the norms that depends on H.
void demodulate_soft_bits(const vec &y, const QLLRvec &LLR_apriori, QLLRvec &LLR_aposteriori, Soft_Demod_Method method=FULL_ENUM_LOGMAP)
Soft MAP demodulation for multidimensional channel, by "brute-force" enumeration of all constellation...
T sum(const Vec< T > &v)
Sum of all elements in the vector.
ivec M
Number of modulation symbols along each dimension.
void modulate_bits(const bvec &bits, cvec &symbols) const
Modulate bits into symbols.
Definitions of matrix inversion routines.
#define it_assert(t, s)
Abort if t is not true.
Logarithmic and exponenential functions - header file.
int length(const Vec< T > &v)
Length of vector.
#define it_assert_debug(t, s)
Abort if t is not true and NDEBUG is not defined.
int nb
Number of bits in the symbol vector.
bool chol(const mat &X, mat &F)
Cholesky factorisation of real symmetric and positive definite matrix.
void set_constellation_points(const int nth, const cvec &inConstellation, const ivec &in_bit2symbols)
Set the constellation points.
void set_size(int n, bool copy=false)
Resizing an Array<T>.
int gray_code(int x)
Convert to Gray Code.
void demodllrbit2(itpp::QLLR &llr) const
Hardcoded implementation of 3:rd bit demodulation.
T prod(const Vec< T > &v)
The product of all elements in the vector.
ITPP_EXPORT vec zeros(int size)
A Double vector of zeros.
Definitions of converters between different vector and matrix types.
int sphere_decoding(const vec &y, const mat &H, double rmin, double rmax, double stepup, QLLRvec &detected_bits)
Sphere decoding.
T min(const Vec< T > &in)
Minimum value of vector.
const double pi
Constant Pi.
ivec k
Number of bits per modulation symbol.
vec log2(const vec &x)
log-2 of the elements
ND_UPAM(int nt=1, int Mary=2)
Constructor.
void update_LLR(const Array< QLLRvec > &logP_apriori, const ivec &s, QLLR scaled_norm, QLLRvec &num, QLLRvec &denom)
Update LLR (for internal use)
T max(const Vec< T > &v)
Maximum value of vector.
Definitions of some specific functions useful in communications.
int pow2i(int x)
Calculate two to the power of x (2^x); x is integer.
Zero-Forcing Log-MAP approximated demodulation.
itpp::QLLRvec Qnorms
Norms part depending on both H and y.
Vec< T > cumsum(const Vec< T > &v)
Cumulative sum of all elements in the vector.
void set_M(int nt=1, int Mary=4)
Set component modulators to M-QAM with Gray mapping.
itpp::QLLRvec llrapr
A prioi information.
int floor_i(double x)
The nearest smaller integer.
void init_soft_demodulator(const itpp::mat &H, const double &sigma2)
Soft MAP demodulation for multidimensional channel, by "brute-force" enumeration of all constellation...
ITPP_EXPORT vec ones(int size)
A float vector of ones.
Miscellaneous statistics functions and classes - header file.
itpp::vec hnorms
Norms part dependent on H.
ITPP_EXPORT bvec dec2bin(int length, int index)
Convert a decimal int index to bvec using length bits in the representation.
bool inv(const mat &X, mat &Y)
Inverse of real square matrix.Calculate the inverse of the real matrix .
void marginalize_bits(itpp::QLLRvec &llr, Soft_Demod_Method method) const
Marginalize (sum) over the bits.
bmat graycode(int m)
Generate Gray code of blocklength m.The codes are contained as binary codewords {0,1} in the rows of the returned matrix. See also the gray() function in math/scalfunc.h.
void init_soft_demodulator(const itpp::cmat &H, const double &sigma2)
Soft MAP demodulation for multidimensional channel, by "brute-force" enumeration of all constellation...
IT++ compatibility types and functions.
vec sqr(const cvec &data)
Absolute square of elements.
Definition of vector (MIMO) modulator classes.
Base class for vector (MIMO) channel modulator/demodulators with complex valued components.
void demodulate_soft_bits(const cvec &y, const QLLRvec &LLR_apriori, QLLRvec &LLR_aposteriori, Soft_Demod_Method method=FULL_ENUM_LOGMAP)
Soft MAP demodulation for multidimensional channel, by "brute-force" enumeration of all constellation...
void modulate_bits(const bvec &bits, vec &symbols) const
Modulate bits into symbols.
Array< vec > get_symbols() const
Get modulation symbols per dimension.
Soft_Demod_Method
Soft demodulation method.
void demodmaxbit2(itpp::QLLR &maxllr) const
Hardcoded implementation of 3:rd bit demodulation.
#define it_error(s)
Abort unconditionally.
vec sqrt(const vec &x)
Square root of the elements.
void set_M(int nt=1, int Mary=2)
Set component modulators to M-PAM with Gray mapping.
void demodllrbit1(itpp::QLLR &llr) const
Hardcoded implementation of 2:nd bit demodulation.
Max-Log demodulation by "brute-force" enumeration of all points.
itpp::Array< itpp::Vec< unsigned > > gray2dec
The Gray to decimal mapping.
bin abs(const bin &inbin)
absolute value of bin
Array< T > left(int n) const
Get n left elements of the array.
Array< bmat > bitmap
Bit mapping table (one table per dimension)
ITPP_EXPORT cvec ones_c(int size)
A float Complex vector of ones.
QLLR jaclog(QLLR a, QLLR b) const
Jacobian logarithm.
void demodllrbit0(itpp::QLLR &llr) const
Hardcoded implementation of 1:st bit demodulation.
ITPP_EXPORT int bin2dec(const bvec &inbvec, bool msb_first=true)
Convert a bvec to decimal int with the first bit as MSB if msb_first == true.
Elementary mathematical functions - header file.
vec real(const cvec &data)
Real part of complex values.
itpp::ivec bpos2cpos
The bit to column mapping.
void demodmaxbit0(itpp::QLLR &maxllr) const
Hardcoded implementation of 1:st bit demodulation.
double gaussnorm
The normalization factor in the exponent (in front of the square norm) in the Gaussian distribution...
int levels2bits(int n)
Calculate the number of bits needed to represent n different values (levels).
void demodmaxbit1(itpp::QLLR &maxllr) const
Hardcoded implementation of 2:nd bit demodulation.
Log-MAP demodulation by "brute-force" enumeration of all points.
Definitions of Cholesky factorisation functions.
const Array< T > concat(const Array< T > &a, const T &e)
Append element e to the end of the Array a.