58 result(1) = -
llrcalc.jaclog(0, -l);
59 result(0) = result(1) - l;
62 result(0) = -
llrcalc.jaclog(0, 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) {
83 p1(b) =
llrcalc.jaclog(p1(b), scaled_norm
84 + log_apriori_prob_const_point);
87 p0(b) =
llrcalc.jaclog(p0(b), scaled_norm
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) {
112 p1(b) =
llrcalc.jaclog(p1(b), scaled_norm
113 + log_apriori_prob_const_point);
116 p0(b) =
llrcalc.jaclog(p0(b), scaled_norm
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;
209 logsum0 =
llrcalc.jaclog(*(Qptr++), logsum0);
210 logsum1 =
llrcalc.jaclog(*(Qptr++), logsum1);
211 logsum1 =
llrcalc.jaclog(*(Qptr++), logsum1);
212 while(Qptr < addr3) {
213 logsum0 =
llrcalc.jaclog(*(Qptr++), logsum0);
214 logsum0 =
llrcalc.jaclog(*(Qptr++), logsum0);
215 logsum1 =
llrcalc.jaclog(*(Qptr++), logsum1);
216 logsum1 =
llrcalc.jaclog(*(Qptr++), logsum1);
218 logsum0 =
llrcalc.jaclog(*(Qptr++), logsum0);
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;
229 logsum0 =
llrcalc.jaclog(*(Qptr++), logsum0);
230 logsum0 =
llrcalc.jaclog(*(Qptr++), logsum0);
231 logsum1 =
llrcalc.jaclog(*(Qptr++), logsum1);
232 logsum1 =
llrcalc.jaclog(*(Qptr++), logsum1);
233 logsum1 =
llrcalc.jaclog(*(Qptr++), logsum1);
234 logsum1 =
llrcalc.jaclog(*(Qptr++), logsum1);
235 while(Qptr < addr3) {
236 logsum0 =
llrcalc.jaclog(*(Qptr++), logsum0);
237 logsum0 =
llrcalc.jaclog(*(Qptr++), logsum0);
238 logsum0 =
llrcalc.jaclog(*(Qptr++), logsum0);
239 logsum0 =
llrcalc.jaclog(*(Qptr++), logsum0);
240 logsum1 =
llrcalc.jaclog(*(Qptr++), logsum1);
241 logsum1 =
llrcalc.jaclog(*(Qptr++), logsum1);
242 logsum1 =
llrcalc.jaclog(*(Qptr++), logsum1);
243 logsum1 =
llrcalc.jaclog(*(Qptr++), logsum1);
245 logsum0 =
llrcalc.jaclog(*(Qptr++), logsum0);
246 logsum0 =
llrcalc.jaclog(*(Qptr++), logsum0);
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;
257 logsum0 =
llrcalc.jaclog(*(Qptr++), logsum0);
258 logsum0 =
llrcalc.jaclog(*(Qptr++), logsum0);
259 logsum0 =
llrcalc.jaclog(*(Qptr++), logsum0);
260 logsum0 =
llrcalc.jaclog(*(Qptr++), logsum0);
261 logsum1 =
llrcalc.jaclog(*(Qptr++), logsum1);
262 logsum1 =
llrcalc.jaclog(*(Qptr++), logsum1);
263 logsum1 =
llrcalc.jaclog(*(Qptr++), logsum1);
264 logsum1 =
llrcalc.jaclog(*(Qptr++), logsum1);
265 logsum1 =
llrcalc.jaclog(*(Qptr++), logsum1);
266 logsum1 =
llrcalc.jaclog(*(Qptr++), logsum1);
267 logsum1 =
llrcalc.jaclog(*(Qptr++), logsum1);
268 logsum1 =
llrcalc.jaclog(*(Qptr++), logsum1);
269 while(Qptr < addr3) {
270 logsum0 =
llrcalc.jaclog(*(Qptr++), logsum0);
271 logsum0 =
llrcalc.jaclog(*(Qptr++), logsum0);
272 logsum0 =
llrcalc.jaclog(*(Qptr++), logsum0);
273 logsum0 =
llrcalc.jaclog(*(Qptr++), logsum0);
274 logsum0 =
llrcalc.jaclog(*(Qptr++), logsum0);
275 logsum0 =
llrcalc.jaclog(*(Qptr++), logsum0);
276 logsum0 =
llrcalc.jaclog(*(Qptr++), logsum0);
277 logsum0 =
llrcalc.jaclog(*(Qptr++), logsum0);
278 logsum1 =
llrcalc.jaclog(*(Qptr++), logsum1);
279 logsum1 =
llrcalc.jaclog(*(Qptr++), logsum1);
280 logsum1 =
llrcalc.jaclog(*(Qptr++), logsum1);
281 logsum1 =
llrcalc.jaclog(*(Qptr++), logsum1);
282 logsum1 =
llrcalc.jaclog(*(Qptr++), logsum1);
283 logsum1 =
llrcalc.jaclog(*(Qptr++), logsum1);
284 logsum1 =
llrcalc.jaclog(*(Qptr++), logsum1);
285 logsum1 =
llrcalc.jaclog(*(Qptr++), logsum1);
287 logsum0 =
llrcalc.jaclog(*(Qptr++), logsum0);
288 logsum0 =
llrcalc.jaclog(*(Qptr++), logsum0);
289 logsum0 =
llrcalc.jaclog(*(Qptr++), logsum0);
290 logsum0 =
llrcalc.jaclog(*(Qptr++), logsum0);
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) {
457 "The length of M vector is different than length of the symbols vector.");
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)));
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");
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;
511 for(
int si = 0; si <
M[ci] - 1; si++) {
518 unsigned bitstring = 0, ind = 0;
526 const itpp::QLLRvec& llr_apr,
528 Soft_Demod_Method method)
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++) {
552 yspacings(ci)[si] = 2*(ytil(ci) * xspacing);
554 unsigned bitstring = 0, ind = 0;
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);
593 const QLLRvec &LLR_apriori,
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));
624 QLLR scaled_norm =
llrcalc.to_qllr(norm2);
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)];
647 bitstring ^= 1 << bit;
655 using namespace itpp;
660 unsigned oldi =
gray2dec(col)[bitstring & (
M[col] - 1)];
662 unsigned newi =
gray2dec(col)[bitstring & (
M[col] - 1)];
671 bitstring ^= 1 << bit;
674 lapr += ((bitstring >> bit) & 1) ? -
llrapr[bit] :
llrapr[bit];
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) {
703 "The length of M vector is different than length of the symbols vector.");
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)));
733 using namespace itpp;
735 it_assert_debug(
sum(
k) < 32,
"Number of total bits per transmission can not be larger than 32.\n");
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;
756 for(
int si = 0; si <
M[ci] - 1; si++) {
763 unsigned bitstring = 0, ind = 0;
764 hxnormupdate(Hx, bitstring, ind,
nb - 1);
769 const itpp::QLLRvec& llr_apr,
771 Soft_Demod_Method method)
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++) {
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) {
817 double sigma_zf = std::sqrt(
real(inv_HhtH(i, i)) * sigma2);
833 const QLLRvec &LLR_apriori,
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));
865 QLLR scaled_norm =
llrcalc.to_qllr(norm2);
866 update_LLR(logP_apriori, j, scaled_norm, i, bnum, bdenom);
868 LLR_aposteriori.set_subvector(b, bnum - bdenom);
874void 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)];
887 hxnormupdate(Hx, bitstring, ind, bit - 1);
889 bitstring ^= 1 << bit;
892 hxnormupdate(Hx, bitstring, ind, bit - 1);
895void 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)];
913 yxnormupdate(yx, lapr, bitstring, ind, bit - 1);
915 bitstring ^= 1 << bit;
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;
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.");
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;
987 spacing(i) = 2.0 / scaling_factor;
991int 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);
1045 z(
k) = std::max(z(
k), zrange(
k, 0));
1046 z(
k) = std::min(z(
k), zrange(
k, 1));
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;
1175 for(
int i = 0; i <
nt; ++i) {
1178 "ND_UQAM::set_M(): M is not a power of 2");
1180 L(i) =
round_i(std::sqrt(
static_cast<double>(
M(i))));
1181 it_assert(
L(i)*
L(i) ==
M(i),
"ND_UQAM: constellation M must be square");
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) {
1193 std::complex<double>(((
L(i) - 1) - j2 * 2.0) / scaling_factor,
1194 ((
L(i) - 1) - j1 * 2.0) / scaling_factor);
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) {
1256 for(
int i = 0; i <
nt; ++i) {
1259 "ND_UPSK::set_M(): M is not a power of 2");
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));
1272 if(
std::abs(std::real(symb)) < epsilon) {
1273 symbols(i)(j) = std::complex<double>(0.0, std::imag(symb));
1275 else if(
std::abs(std::imag(symb)) < epsilon) {
1276 symbols(i)(j) = std::complex<double>(std::real(symb), 0.0);
Definitions of Cholesky factorisation functions.
Array< T > left(int n) const
Get n left elements of the array.
void set_size(int n, bool copy=false)
Resizing an Array<T>.
QLLR jaclog(QLLR a, QLLR b) const
Jacobian logarithm.
Modulator_NCD()
Constructor.
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...
itpp::Array< itpp::Array< itpp::cvec > > hspacings
The spacing between different constellation points multiplied by the different H columns.
void init_soft_demodulator(const itpp::cmat &H, const double &sigma2)
Soft MAP demodulation for multidimensional channel, by "brute-force" enumeration of all constellation...
itpp::Array< itpp::vec > yspacings
The spacing between different constellation points scaled by different y elements.
itpp::cmat H
Complex-valued channel matrix.
Array< cvec > get_symbols() const
Get modulation symbols per dimension.
void modulate_bits(const bvec &bits, cvec &symbols) const
Modulate bits into symbols.
Array< cvec > symbols
Vectors of modulation symbols (along each dimension)
itpp::QLLRvec Qnorms
Norms part depending on both H and y.
double gaussnorm
The normalization factor in the exponent (in front of the square norm) in the Gaussian distribution.
Array< ivec > bits2symbols
Bit pattern in decimal form ordered and the corresponding symbols (one pattern per dimension)
int nt
Number of dimensions.
void demodllrbit2(itpp::QLLR &llr) const
Hardcoded implementation of 3:rd bit demodulation.
itpp::vec hnorms
Norms part dependent on H.
void demodmaxbit1(itpp::QLLR &maxllr) const
Hardcoded implementation of 2:nd bit demodulation.
ivec M
Number of modulation symbols along each dimension.
void update_LLR(const Array< QLLRvec > &logP_apriori, const ivec &s, QLLR scaled_norm, QLLRvec &num, QLLRvec &denom)
Update LLR (for internal use)
itpp::ivec bitcumsum
The cumulative sum of bits in the symbol vector.
ivec k
Number of bits per modulation symbol.
itpp::QLLRvec llrapr
A prioi information.
Soft_Demod_Method
Soft demodulation method.
@ FULL_ENUM_MAXLOG
Max-Log demodulation by "brute-force" enumeration of all points.
@ FULL_ENUM_LOGMAP
Log-MAP demodulation by "brute-force" enumeration of all points.
@ ZF_LOGMAP
Zero-Forcing Log-MAP approximated demodulation.
void demodllrbit1(itpp::QLLR &llr) const
Hardcoded implementation of 2:nd bit demodulation.
void marginalize_bits(itpp::QLLRvec &llr, Soft_Demod_Method method) const
Marginalize (sum) over the bits.
LLR_calc_unit llrcalc
LLR calculation unit.
itpp::Array< itpp::Vec< unsigned > > gray2dec
The Gray to decimal mapping.
int nb
Number of bits in the symbol vector.
QLLRvec probabilities(QLLR l)
Convert LLR to log-probabilities.
void demodmaxbit0(itpp::QLLR &maxllr) const
Hardcoded implementation of 1:st bit demodulation.
Array< bmat > bitmap
Bit mapping table (one table per dimension)
void demodllrbit0(itpp::QLLR &llr) const
Hardcoded implementation of 1:st bit demodulation.
void demodmaxbit2(itpp::QLLR &maxllr) const
Hardcoded implementation of 3:rd bit demodulation.
bool demod_initialized
Flag indicating whether the demodulator has been initialized.
itpp::ivec bpos2cpos
The bit to column mapping.
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...
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.
void hxnormupdate(itpp::vec &Hx, unsigned &bitstring, unsigned &ind, unsigned bit)
Calculation of the part of the norms that depends on H.
itpp::Array< itpp::vec > yspacings
The spacing between different constellation points scaled by different y elements.
itpp::mat H
Real channel matrix.
Array< vec > symbols
Vectors of modulation symbols (along each dimension)
void modulate_bits(const bvec &bits, vec &symbols) const
Modulate bits into symbols.
Array< vec > get_symbols() const
Get modulation symbols per dimension.
Modulator_NRD()
Constructor.
itpp::Array< itpp::Array< itpp::vec > > hspacings
The spacing between different constellation points multiplied by the different H columns.
void init_soft_demodulator(const itpp::mat &H, const double &sigma2)
Soft MAP demodulation for multidimensional channel, by "brute-force" enumeration of all constellation...
ND_UPAM(int nt=1, int Mary=2)
Constructor.
int sphere_decoding(const vec &y, const mat &H, double rmin, double rmax, double stepup, QLLRvec &detected_bits)
Sphere decoding.
void set_M(int nt=1, int Mary=2)
Set component modulators to M-PAM with Gray mapping.
void set_M(int nt=1, int Mary=4)
Set component modulators to M-QAM with Gray mapping.
ND_UPSK(int nt=1, int Mary=4)
Constructor.
ivec L
the square root of M
void set_constellation_points(const int nth, const cvec &inConstellation, const ivec &in_bit2symbols)
Set the constellation points.
ND_UQAM(int nt=1, int Mary=4)
Constructor.
void set_M(int nt=1, int Mary=4)
Set component modulators to M-QAM with Gray mapping.
Definitions of some specific functions useful in communications.
Definitions of converters between different vector and matrix types.
Elementary mathematical functions - header file.
#define it_error(s)
Abort unconditionally.
#define it_assert_debug(t, s)
Abort if t is not true and NDEBUG is not defined.
#define it_assert(t, s)
Abort if t is not true.
bool inv(const mat &X, mat &Y)
Inverse of real square matrix.
int pow2i(int x)
Calculate two to the power of x (2^x); x is integer.
vec log2(const vec &x)
log-2 of the elements
int levels2bits(int n)
Calculate the number of bits needed to represent n different values (levels).
Vec< T > cumsum(const Vec< T > &v)
Cumulative sum of all elements in the vector.
T sum(const Vec< T > &v)
Sum of all elements in the vector.
int length(const Vec< T > &v)
Length of vector.
T prod(const Vec< T > &v)
The product of all elements in the vector.
bool chol(const mat &X, mat &F)
Cholesky factorisation of real symmetric and positive definite matrix.
bmat graycode(int m)
Generate Gray code of blocklength m.
int mod(int k, int n)
Calculates the modulus, i.e. the signed reminder after division.
bool is_even(int x)
Return true if x is an even integer.
vec real(const cvec &data)
Real part of complex values.
vec sqr(const cvec &data)
Absolute square of elements.
cvec conj(const cvec &x)
Conjugate of complex value.
Vec< T > reverse(const Vec< T > &in)
Reverse the input vector.
ITPP_EXPORT vec zeros(int size)
A Double vector of zeros.
ITPP_EXPORT ivec ones_i(int size)
A Int vector of ones.
ITPP_EXPORT cvec ones_c(int size)
A float Complex vector of ones.
ITPP_EXPORT vec ones(int size)
A float vector of ones.
ITPP_EXPORT ivec zeros_i(int size)
A Int vector of zeros.
double norm(const cvec &v)
Calculate the 2-norm: norm(v)=sqrt(sum(abs(v).^2))
Definitions of matrix inversion routines.
IT++ compatibility types and functions.
Logarithmic and exponenential functions - header file.
Mat< bin > bmat
bin matrix
Various functions on vectors and matrices - header file.
Miscellaneous statistics functions and classes - header file.
Definition of vector (MIMO) modulator classes.
ITPP_EXPORT int round_i(double x)
Round to nearest integer.
std::ostream & operator<<(std::ostream &output, const bin &inbin)
Output stream of bin.
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.
const Array< T > concat(const Array< T > &a, const T &e)
Append element e to the end of the Array a.
int gray_code(int x)
Convert to Gray Code.
ITPP_EXPORT bvec dec2bin(int length, int index)
Convert a decimal int index to bvec using length bits in the representation.
const double pi
Constant Pi.
int floor_i(double x)
The nearest smaller integer.
int abs(const itpp::bin &inbin)
absolute value of bin