IT++ Logo
modulator_nd.cpp
Go to the documentation of this file.
1
30#include <itpp/comm/commfunc.h>
33#include <itpp/base/matfunc.h>
37#include <itpp/base/itcompat.h>
38#include <itpp/base/sort.h>
39#include <itpp/stat/misc_stat.h>
40
41#include <cmath>
42#include <iostream>
43#include <iomanip>
44
45namespace itpp
46{
47
48// ----------------------------------------------------------------------
49// Modulator_ND
50// ----------------------------------------------------------------------
51
52
54{
55 QLLRvec result(2);
56
57 if(l < 0) { // this can be done more efficiently
58 result(1) = -llrcalc.jaclog(0, -l);
59 result(0) = result(1) - l;
60 }
61 else {
62 result(0) = -llrcalc.jaclog(0, l);
63 result(1) = result(0) + l;
64 }
65 return result;
66}
67
68 void Modulator_ND::update_LLR(const Array<QLLRvec> &logP_apriori, int s,
69 QLLR scaled_norm, int j, QLLRvec &p1,
70 QLLRvec &p0)
71{
72 QLLR log_apriori_prob_const_point = 0;
73 int b = 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));
77 b++;
78 }
79
80 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);
85 }
86 else {
87 p0(b) = llrcalc.jaclog(p0(b), scaled_norm
88 + log_apriori_prob_const_point);
89 }
90 b++;
91 }
92}
93
95 const ivec &s, QLLR scaled_norm,
96 QLLRvec &p1, QLLRvec &p0)
97{
98 QLLR log_apriori_prob_const_point = 0;
99 int b = 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));
104 b++;
105 }
106 }
107
108 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);
114 }
115 else {
116 p0(b) = llrcalc.jaclog(p0(b), scaled_norm
117 + log_apriori_prob_const_point);
118 }
119 b++;
120 }
121 }
122}
123
124 void Modulator_ND::marginalize_bits(itpp::QLLRvec& llr, Soft_Demod_Method method) const
125{
126 if(method == FULL_ENUM_LOGMAP) {
127 // -- Demodulate the last 3 bits. The demodulation is hardcoded
128 // -- to avoid initialization of many but tiny inner-loops
129 demodllrbit0(llr[0]);
130 if(nb > 1) demodllrbit1(llr[1]);
131 if(nb > 2) demodllrbit2(llr[2]);
132 // -- Demodulate the remaining bits except the first one
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);
136 const QLLR *Qptr;
137 for(int bi = 3; bi < nb - 1 ; bi++) { // Run the loops for bits 3,...,nb-1.
138 logsum0 = -QLLR_MAX;
139 logsum1 = -QLLR_MAX;
140 const int forhalfdiff = 1 << bi, fordiff = 2 * forhalfdiff, fordbldiff = 2 * fordiff;
141 Qptr = addrfirst;
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);
151 }
152 while(Qptr < addrlast) logsum0 = llrcalc.jaclog(*(Qptr++), logsum0);
153 llr[bi] = logsum0 - logsum1;
154 }
155 // -- Demodulate the first bit
156 logsum0 = -QLLR_MAX;
157 logsum1 = -QLLR_MAX;
158 Qptr = addrfirst;
159 while(Qptr < addrsemilast) logsum0 = llrcalc.jaclog(*(Qptr++), logsum0);
160 while(Qptr < addrlast) logsum1 = llrcalc.jaclog(*(Qptr++), logsum1);
161 llr[nb - 1] = logsum0 - logsum1;
162 }
163 else if(method == FULL_ENUM_MAXLOG) {
164 // -- Demodulate the last 3 bits. The demodulation is hardcoded to
165 // -- avoid initialization of many but tiny inner-loops.
166 demodmaxbit0(llr[0]);
167 if(nb > 1) demodmaxbit1(llr[1]);
168 if(nb > 2) demodmaxbit2(llr[2]);
169 // -- Demodulate the remaining bits except the first one
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);
173 const QLLR *Qptr;
174 for(int bi = 3; bi < nb - 1; bi++) { // Run the loops for bits nb-3,nb-4,...,2.
175 logmax0 = -QLLR_MAX;
176 logmax1 = -QLLR_MAX;
177 const int forhalfdiff = 1 << bi, fordiff = 2 * forhalfdiff, fordbldiff = 2 * fordiff;
178 Qptr = addrfirst;
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;
188 }
189 for(; Qptr < addrlast; Qptr++) logmax0 = *Qptr > logmax0 ? *Qptr : logmax0;
190 llr[bi] = logmax0 - logmax1;
191 }
192 // -- Demodulate the first bit
193 logmax0 = -QLLR_MAX;
194 logmax1 = -QLLR_MAX;
195 Qptr = addrfirst;
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;
199 }
200 else it_error("Improper soft demodulation method\n.");
201}
202
203 void Modulator_ND::demodllrbit0(itpp::QLLR& llr) const
204{
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);
217 }
218 logsum0 = llrcalc.jaclog(*(Qptr++), logsum0);
219 llr = logsum0 - logsum1;
220}
221
222void Modulator_ND::demodllrbit1(itpp::QLLR& llr) const
223{
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;
228
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);
244 }
245 logsum0 = llrcalc.jaclog(*(Qptr++), logsum0);
246 logsum0 = llrcalc.jaclog(*(Qptr++), logsum0);
247 llr = logsum0 - logsum1;
248}
249
250void Modulator_ND::demodllrbit2(itpp::QLLR& llr) const
251{
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;
256
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);
286 }
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;
292}
293
294void Modulator_ND::demodmaxbit0(itpp::QLLR& maxllr) const
295{
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;
301 Qptr++;
302 logmax1 = *Qptr > logmax1 ? *Qptr : logmax1;
303 Qptr++;
304 logmax1 = *Qptr > logmax1 ? *Qptr : logmax1;
305 Qptr++;
306 while(Qptr < addr3) {
307 logmax0 = *Qptr > logmax0 ? *Qptr : logmax0;
308 Qptr++;
309 logmax0 = *Qptr > logmax0 ? *Qptr : logmax0;
310 Qptr++;
311 logmax1 = *Qptr > logmax1 ? *Qptr : logmax1;
312 Qptr++;
313 logmax1 = *Qptr > logmax1 ? *Qptr : logmax1;
314 Qptr++;
315 }
316 logmax0 = *Qptr > logmax0 ? *Qptr : logmax0;
317 maxllr = logmax0 - logmax1;
318}
319
320void Modulator_ND::demodmaxbit1(itpp::QLLR& maxllr) const
321{
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;
327 Qptr++;
328 logmax0 = *Qptr > logmax0 ? *Qptr : logmax0;
329 Qptr++;
330 logmax1 = *Qptr > logmax1 ? *Qptr : logmax1;
331 Qptr++;
332 logmax1 = *Qptr > logmax1 ? *Qptr : logmax1;
333 Qptr++;
334 logmax1 = *Qptr > logmax1 ? *Qptr : logmax1;
335 Qptr++;
336 logmax1 = *Qptr > logmax1 ? *Qptr : logmax1;
337 Qptr++;
338 while(Qptr < addr3) {
339 logmax0 = *Qptr > logmax0 ? *Qptr : logmax0;
340 Qptr++;
341 logmax0 = *Qptr > logmax0 ? *Qptr : logmax0;
342 Qptr++;
343 logmax0 = *Qptr > logmax0 ? *Qptr : logmax0;
344 Qptr++;
345 logmax0 = *Qptr > logmax0 ? *Qptr : logmax0;
346 Qptr++;
347 logmax1 = *Qptr > logmax1 ? *Qptr : logmax1;
348 Qptr++;
349 logmax1 = *Qptr > logmax1 ? *Qptr : logmax1;
350 Qptr++;
351 logmax1 = *Qptr > logmax1 ? *Qptr : logmax1;
352 Qptr++;
353 logmax1 = *Qptr > logmax1 ? *Qptr : logmax1;
354 Qptr++;
355 }
356 logmax0 = *Qptr > logmax0 ? *Qptr : logmax0;
357 Qptr++;
358 logmax0 = *Qptr > logmax0 ? *Qptr : logmax0;
359 maxllr = logmax0 - logmax1;
360}
361
362void Modulator_ND::demodmaxbit2(itpp::QLLR& maxllr) const
363{
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;
369 Qptr++;
370 logmax0 = *Qptr > logmax0 ? *Qptr : logmax0;
371 Qptr++;
372 logmax0 = *Qptr > logmax0 ? *Qptr : logmax0;
373 Qptr++;
374 logmax0 = *Qptr > logmax0 ? *Qptr : logmax0;
375 Qptr++;
376 logmax1 = *Qptr > logmax1 ? *Qptr : logmax1;
377 Qptr++;
378 logmax1 = *Qptr > logmax1 ? *Qptr : logmax1;
379 Qptr++;
380 logmax1 = *Qptr > logmax1 ? *Qptr : logmax1;
381 Qptr++;
382 logmax1 = *Qptr > logmax1 ? *Qptr : logmax1;
383 Qptr++;
384 logmax1 = *Qptr > logmax1 ? *Qptr : logmax1;
385 Qptr++;
386 logmax1 = *Qptr > logmax1 ? *Qptr : logmax1;
387 Qptr++;
388 logmax1 = *Qptr > logmax1 ? *Qptr : logmax1;
389 Qptr++;
390 logmax1 = *Qptr > logmax1 ? *Qptr : logmax1;
391 Qptr++;
392 while(Qptr < addr3) {
393 logmax0 = *Qptr > logmax0 ? *Qptr : logmax0;
394 Qptr++;
395 logmax0 = *Qptr > logmax0 ? *Qptr : logmax0;
396 Qptr++;
397 logmax0 = *Qptr > logmax0 ? *Qptr : logmax0;
398 Qptr++;
399 logmax0 = *Qptr > logmax0 ? *Qptr : logmax0;
400 Qptr++;
401 logmax0 = *Qptr > logmax0 ? *Qptr : logmax0;
402 Qptr++;
403 logmax0 = *Qptr > logmax0 ? *Qptr : logmax0;
404 Qptr++;
405 logmax0 = *Qptr > logmax0 ? *Qptr : logmax0;
406 Qptr++;
407 logmax0 = *Qptr > logmax0 ? *Qptr : logmax0;
408 Qptr++;
409 logmax1 = *Qptr > logmax1 ? *Qptr : logmax1;
410 Qptr++;
411 logmax1 = *Qptr > logmax1 ? *Qptr : logmax1;
412 Qptr++;
413 logmax1 = *Qptr > logmax1 ? *Qptr : logmax1;
414 Qptr++;
415 logmax1 = *Qptr > logmax1 ? *Qptr : logmax1;
416 Qptr++;
417 logmax1 = *Qptr > logmax1 ? *Qptr : logmax1;
418 Qptr++;
419 logmax1 = *Qptr > logmax1 ? *Qptr : logmax1;
420 Qptr++;
421 logmax1 = *Qptr > logmax1 ? *Qptr : logmax1;
422 Qptr++;
423 logmax1 = *Qptr > logmax1 ? *Qptr : logmax1;
424 Qptr++;
425 }
426 logmax0 = *Qptr > logmax0 ? *Qptr : logmax0;
427 Qptr++;
428 logmax0 = *Qptr > logmax0 ? *Qptr : logmax0;
429 Qptr++;
430 logmax0 = *Qptr > logmax0 ? *Qptr : logmax0;
431 Qptr++;
432 logmax0 = *Qptr > logmax0 ? *Qptr : logmax0;
433 maxllr = logmax0 - logmax1;
434}
435
436
438{
439 Array<QLLRvec> result(length(l));
440 for(int i = 0; i < length(l); i++) {
441 result(i) = probabilities(l(i));
442 }
443 return result;
444}
445
446
447// ----------------------------------------------------------------------
448// Modulator_NRD
449// ----------------------------------------------------------------------
450
451
453{
454 Array<vec> retvec(nt);
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));
459 }
460 return retvec;
461}
462
463void Modulator_NRD::modulate_bits(const bvec &bits, vec &out_symbols) const
464{
465 it_assert(length(bits) == sum(k), "Modulator_NRD::modulate_bits(): "
466 "The number of input bits does not match.");
467
468 out_symbols.set_size(nt);
469
470 int b = 0;
471 for(int i = 0; i < nt; ++i) {
472 int symb = bin2dec(bits.mid(b, k(i)));
473 out_symbols(i) = symbols(i)(bits2symbols(i)(symb));
474 b += k(i);
475 }
476}
477
478vec Modulator_NRD::modulate_bits(const bvec &bits) const
479{
480 vec result(nt);
481 modulate_bits(bits, result);
482 return result;
483}
484
485
486void Modulator_NRD::init_soft_demodulator(const itpp::mat& H_in, const double& sigma2)
487{
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");
492 H = H_in;
493 bitcumsum = reverse(cumsum(reverse(k)) - reverse(k)); // Shifted cummulative sum
494 nb = sum(k);
495 hnorms.set_size(1 << nb);
496 Qnorms.set_size(1 << nb);
499 bpos2cpos.set_size(nb);
501 gaussnorm = 2 * sigma2;
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;
507 gray2dec(ci).set_size(M[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;
515 }
516 }
518 unsigned bitstring = 0, ind = 0;
519 hxnormupdate(Hx, bitstring, ind, nb - 1);
520 demod_initialized = true;
521}
522
523
524
525void Modulator_NRD::demodulate_soft_bits(const itpp::vec& y,
526 const itpp::QLLRvec& llr_apr,
527 itpp::QLLRvec& llr,
528 Soft_Demod_Method method)
529{
530 using namespace itpp;
531
532 it_assert_debug(demod_initialized, "You have to first run init_soft_demodulator().\n");
533 it_assert_debug(H.rows() == y.length(), "The dimensions are not correct.\n");
534 it_assert_debug(llr_apr.length() == nb, "The LLR_apr length is not correct.\n");
535
536 // -- Prepare all the norms with the newly received vectory y
537 llr.set_size(nb);
538 llrapr = reverse(llr_apr); /* The bits are reversed due to the
539 norm-updating functions having the rightmost bit
540 as the least significant*/
541
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);
546 QLLR lapr = 0;
547 for(int bi = 0; bi < nb; lapr -= llrcalc.jaclog(0, -llrapr[bi++]));
548
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);
553 }
554 unsigned bitstring = 0, ind = 0;
555 yxnormupdate(yx, lapr, bitstring, ind, nb - 1); // Recursive update of all the norms
556 marginalize_bits(llr,method); // Perform the appropriate bit marginalization
557 llr = reverse(llr);
558}
559
560void Modulator_NRD::demodulate_soft_bits(const vec &y, const mat &H,
561 double sigma2,
562 const QLLRvec &LLR_apriori,
563 QLLRvec &LLR_aposteriori,
564 Soft_Demod_Method method)
565{
566 switch(method) {
567 case ZF_LOGMAP: {
568 it_assert(H.rows() >= H.cols(), "Modulator_NRD::demodulate_soft_bits():"
569 " ZF demodulation impossible for undetermined systems");
570 // Set up the ZF detector
571 mat Ht = H.T();
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) {
576 // noise covariance of shat
577 double sigma_zf = std::sqrt(inv_HtH(i, i) * sigma2);
578 shat(i) /= sigma_zf;
579 h(i) /= sigma_zf;
580 }
581 demodulate_soft_bits(shat, h, 1.0, zeros_i(sum(k)), LLR_aposteriori);
582 }
583 break;
584 default: {
585 init_soft_demodulator(H, sigma2);
586 demodulate_soft_bits(y, LLR_apriori, LLR_aposteriori, method);
587 }
588 }
589}
590
591QLLRvec Modulator_NRD::demodulate_soft_bits(const vec &y, const mat &H,
592 double sigma2,
593 const QLLRvec &LLR_apriori,
594 Soft_Demod_Method method)
595{
596 QLLRvec result;
597 demodulate_soft_bits(y, H, sigma2, LLR_apriori, result, method);
598 return result;
599}
600
601void Modulator_NRD::demodulate_soft_bits(const vec &y, const vec &h,
602 double sigma2,
603 const QLLRvec &LLR_apriori,
604 QLLRvec &LLR_aposteriori)
605{
606 it_assert(length(LLR_apriori) == sum(k),
607 "Modulator_NRD::demodulate_soft_bits(): Wrong sizes");
608 it_assert((length(h) == length(y)) && (length(h) == nt),
609 "Modulator_NRD::demodulate_soft_bits(): Wrong sizes");
610
611 // set size of the output vector
612 LLR_aposteriori.set_size(LLR_apriori.size());
613
614 // normalisation constant "minus one over two sigma^2"
615 double moo2s2 = -1.0 / (2.0 * sigma2);
616
617 int b = 0;
618 for(int i = 0; i < nt; ++i) {
619 QLLRvec bnum = -QLLR_MAX * ones_i(k(i));
620 QLLRvec bdenom = bnum;
621 Array<QLLRvec> logP_apriori = probabilities(LLR_apriori(b, b + k(i) - 1));
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);
626 }
627 LLR_aposteriori.set_subvector(b, bnum - bdenom);
628 b += k(i);
629 }
630}
631
632void Modulator_NRD::hxnormupdate(itpp::vec& Hx, unsigned& bitstring, unsigned& ind, unsigned bit)
633{
634 using namespace itpp;
635 const unsigned col = bpos2cpos[bit];
636 if(bit < 1) {
637 hnorms[ind++] = Hx * Hx;
638 unsigned oldi = gray2dec(col)[bitstring & (M[col] - 1)];
639 bitstring ^= 1;
640 unsigned newi = gray2dec(col)[bitstring & (M[col] - 1)];
641 Hx += oldi > newi ? -hspacings(col)(newi) : hspacings(col)(oldi);
642 hnorms[ind++] = Hx * Hx;
643 return;
644 }
645 hxnormupdate(Hx, bitstring, ind, bit - 1);
646 unsigned oldi = gray2dec(col)[(bitstring >> bitcumsum[col]) & (M[col] - 1)];
647 bitstring ^= 1 << bit;
648 unsigned newi = gray2dec(col)[(bitstring >> bitcumsum[col]) & (M[col] - 1)];
649 Hx += oldi > newi ? -hspacings(col)(newi) : hspacings(col)(oldi);
650 hxnormupdate(Hx, bitstring, ind, bit - 1);
651}
652
653void Modulator_NRD::yxnormupdate(double& yx, itpp::QLLR& lapr, unsigned& bitstring, unsigned& ind, unsigned bit)
654{
655 using namespace itpp;
656 const unsigned col = bpos2cpos[bit];
657 if(bit < 1) {
658 Qnorms[ind] = llrcalc.to_qllr((yx - hnorms[ind]) / gaussnorm) + lapr;
659 ind++;
660 unsigned oldi = gray2dec(col)[bitstring & (M[col] - 1)];
661 bitstring ^= 1;
662 unsigned newi = gray2dec(col)[bitstring & (M[col] - 1)];
663 yx += oldi > newi ? -yspacings(col)[newi] : yspacings(col)[oldi];
664 lapr += (bitstring & 1) ? -llrapr[bit] : llrapr[bit];
665 Qnorms[ind] = llrcalc.to_qllr((yx - hnorms[ind]) / gaussnorm) + lapr;
666 ind++;
667 return;
668 }
669 yxnormupdate(yx, lapr, bitstring, ind, bit - 1);
670 unsigned oldi = gray2dec(col)[(bitstring >> bitcumsum[col]) & (M[col] - 1)];
671 bitstring ^= 1 << bit;
672 unsigned newi = gray2dec(col)[(bitstring >> bitcumsum[col]) & (M[col] - 1)];
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);
676}
677
678
679std::ostream &operator<<(std::ostream &os, const Modulator_NRD &mod)
680{
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;
687 // skip printing the trailing zero
688 os << "Symbol coordinates for dimension " << i << ": " << mod.symbols(i).left(mod.M(i)) << std::endl;
689 }
690 os << mod.get_llrcalc() << std::endl;
691 return os;
692}
693
694// ----------------------------------------------------------------------
695// Modulator_NCD
696// ----------------------------------------------------------------------
697
699{
700 Array<cvec> retvec(nt);
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));
705 }
706 return retvec;
707}
708
709void Modulator_NCD::modulate_bits(const bvec &bits, cvec &out_symbols) const
710{
711 it_assert(length(bits) == sum(k), "Modulator_NCD::modulate_bits(): "
712 "The number of input bits does not match.");
713
714 out_symbols.set_size(nt);
715
716 int b = 0;
717 for(int i = 0; i < nt; ++i) {
718 int symb = bin2dec(bits.mid(b, k(i)));
719 out_symbols(i) = symbols(i)(bits2symbols(i)(symb));
720 b += k(i);
721 }
722}
723
724cvec Modulator_NCD::modulate_bits(const bvec &bits) const
725{
726 cvec result(nt);
727 modulate_bits(bits, result);
728 return result;
729}
730
731void Modulator_NCD::init_soft_demodulator(const itpp::cmat& H_in, const double& sigma2)
732{
733 using namespace itpp;
734 it_assert_debug(H_in.cols() == nt, "The number of Tx antennas is wrong.\n");
735 it_assert_debug(sum(k) < 32, "Number of total bits per transmission can not be larger than 32.\n");
736 it_assert_debug(pow2i(sum(k)) == prod(M), "The modulater must use exhaustive constellations, i.e., #bits=log2(#symbs).\n");
737 H = H_in;
738 bitcumsum = reverse(cumsum(reverse(k)) - reverse(k)); // Shifted cummulative sum
739 nb = sum(k);
740 hnorms.set_size(1 << nb);
741 Qnorms.set_size(1 << nb);
744 bpos2cpos.set_size(nb);
746 gaussnorm = sigma2;
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;
752 gray2dec(ci).set_size(M[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;
760 }
761 }
763 unsigned bitstring = 0, ind = 0;
764 hxnormupdate(Hx, bitstring, ind, nb - 1);
765 demod_initialized = true;
766}
767
768void Modulator_NCD::demodulate_soft_bits(const itpp::cvec& y,
769 const itpp::QLLRvec& llr_apr,
770 itpp::QLLRvec& llr,
771 Soft_Demod_Method method)
772{
773 using namespace itpp;
774
775 it_assert_debug(demod_initialized, "You have to first run init_soft_demodulator().\n");
776 it_assert_debug(H.rows() == y.length(), "The dimensions are not correct.\n");
777 it_assert_debug(llr_apr.length() == nb, "The LLR_apr length is not correct.\n");
778
779 // -- Prepare all the norms with the newly received vectory y
780 llr.set_size(nb);
781 llrapr = reverse(llr_apr); /* The bits are reversed due to the
782 norm-updating functions having the rightmost bit
783 as the least significant*/
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();
788 QLLR lapr = 0;
789 for(int bi = 0; bi < nb; lapr -= llrcalc.jaclog(0, -llrapr[bi++]));
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();
794 }
795 unsigned bitstring = 0, ind = 0;
796 yxnormupdate(yx, lapr, bitstring, ind, nb - 1); // Recursive update of all the norms
797 marginalize_bits(llr,method);
798 llr=reverse(llr);
799}
800
801void Modulator_NCD::demodulate_soft_bits(const cvec &y, const cmat &H,
802 double sigma2,
803 const QLLRvec &LLR_apriori,
804 QLLRvec &LLR_aposteriori,
805 Soft_Demod_Method method)
806{
807 switch(method) {
808 case ZF_LOGMAP: {
809 it_assert(H.rows() >= H.cols(), "Modulator_NCD::demodulate_soft_bits():"
810 " ZF demodulation impossible for undetermined systems");
811 // Set up the ZF detector
812 cmat Hht = H.H();
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);
818 shat(i) /= sigma_zf;
819 h(i) /= sigma_zf;
820 }
821 demodulate_soft_bits(shat, h, 1.0, zeros_i(sum(k)), LLR_aposteriori);
822 }
823 break;
824 default: {
825 init_soft_demodulator(H, sigma2);
826 demodulate_soft_bits(y, LLR_apriori, LLR_aposteriori, method);
827 }
828 }
829}
830
831QLLRvec Modulator_NCD::demodulate_soft_bits(const cvec &y, const cmat &H,
832 double sigma2,
833 const QLLRvec &LLR_apriori,
834 Soft_Demod_Method method)
835{
836 QLLRvec result;
837 demodulate_soft_bits(y, H, sigma2, LLR_apriori, result, method);
838 return result;
839}
840
841
842void Modulator_NCD::demodulate_soft_bits(const cvec &y, const cvec &h,
843 double sigma2,
844 const QLLRvec &LLR_apriori,
845 QLLRvec &LLR_aposteriori)
846{
847 it_assert(length(LLR_apriori) == sum(k),
848 "Modulator_NCD::demodulate_soft_bits(): Wrong sizes");
849 it_assert((length(h) == length(y)) && (length(h) == nt),
850 "Modulator_NCD::demodulate_soft_bits(): Wrong sizes");
851
852 // set size of the output vector
853 LLR_aposteriori.set_size(LLR_apriori.size());
854
855 // normalisation constant "minus one over sigma^2"
856 double moos2 = -1.0 / sigma2;
857
858 int b = 0;
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));
862 Array<QLLRvec> logP_apriori = probabilities(LLR_apriori(b, b + k(i) - 1));
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);
867 }
868 LLR_aposteriori.set_subvector(b, bnum - bdenom);
869 b += k(i);
870 }
871}
872
873
874void Modulator_NCD::hxnormupdate(itpp::cvec& Hx, unsigned& bitstring, unsigned& ind, unsigned bit)
875{
876 using namespace itpp;
877 const unsigned col = bpos2cpos[bit];
878 if(bit < 1) {
879 hnorms[ind++] = sqr(norm(Hx));
880 unsigned oldi = gray2dec(col)[bitstring & (M[col] - 1)];
881 bitstring ^= 1;
882 unsigned newi = gray2dec(col)[bitstring & (M[col] - 1)];
883 Hx += oldi > newi ? -hspacings(col)(newi) : hspacings(col)(oldi);
884 hnorms[ind++] = sqr(norm(Hx));
885 return;
886 }
887 hxnormupdate(Hx, bitstring, ind, bit - 1);
888 unsigned oldi = gray2dec(col)[(bitstring >> bitcumsum[col]) & (M[col] - 1)];
889 bitstring ^= 1 << bit;
890 unsigned newi = gray2dec(col)[(bitstring >> bitcumsum[col]) & (M[col] - 1)];
891 Hx += oldi > newi ? -hspacings(col)(newi) : hspacings(col)(oldi);
892 hxnormupdate(Hx, bitstring, ind, bit - 1);
893}
894
895void Modulator_NCD::yxnormupdate(double& yx, itpp::QLLR& lapr, unsigned& bitstring, unsigned& ind, unsigned bit)
896{
897 using namespace itpp;
898 const unsigned col = bpos2cpos[bit];
899 if(bit < 1) {
900 Qnorms[ind] = llrcalc.to_qllr((yx - hnorms[ind]) / gaussnorm) + lapr;
901 //std::cerr << dec2bin(sum(k),(int)bitstring) << " " << Qnorms[ind] << " "
902 // << llrcalc.to_qllr((2*(rec.H()*H*modulate_bits(dec2bin(sum(k),(int)bitstring)))[0].real() - hnorms[ind]) / gaussnorm) + lapr << std::endl;
903 ind++;
904 unsigned oldi = gray2dec(col)[bitstring & (M[col] - 1)];
905 bitstring ^= 1;
906 unsigned newi = gray2dec(col)[bitstring & (M[col] - 1)];
907 yx += oldi > newi ? -yspacings(col)[newi] : yspacings(col)[oldi];
908 lapr += (bitstring & 1) ? -llrapr[bit] : llrapr[bit];
909 Qnorms[ind] = llrcalc.to_qllr((yx - hnorms[ind]) / gaussnorm) + lapr;
910 ind++;
911 return;
912 }
913 yxnormupdate(yx, lapr, bitstring, ind, bit - 1);
914 unsigned oldi = gray2dec(col)[(bitstring >> bitcumsum[col]) & (M[col] - 1)];
915 bitstring ^= 1 << bit;
916 unsigned newi = gray2dec(col)[(bitstring >> bitcumsum[col]) & (M[col] - 1)];
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);
920}
921
922
923std::ostream &operator<<(std::ostream &os, const Modulator_NCD &mod)
924{
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;
934 }
935 os << mod.get_llrcalc() << std::endl;
936 return os;
937}
938
939// ----------------------------------------------------------------------
940// ND_UPAM
941// ----------------------------------------------------------------------
942
943ND_UPAM::ND_UPAM(int nt, int Mary)
944{
945 set_M(nt, Mary);
946}
947
948void ND_UPAM::set_M(int nt_in, int Mary)
949{
950 nt = nt_in;
951 ivec Mary_temp(nt);
952 Mary_temp = Mary;
953 set_M(nt, Mary_temp);
954}
955
956void ND_UPAM::set_M(int nt_in, ivec Mary)
957{
958 nt = nt_in;
959 it_assert(length(Mary) == nt, "ND_UPAM::set_M(): Mary has wrong length");
960 k.set_size(nt);
961 M = Mary;
965 spacing.set_size(nt);
966
967 for(int i = 0; i < nt; i++) {
968 k(i) = round_i(::log2(static_cast<double>(M(i))));
969 it_assert((k(i) > 0) && ((1 << k(i)) == M(i)),
970 "ND_UPAM::set_M(): M is not a power of 2.");
971
972 symbols(i).set_size(M(i) + 1);
973 bits2symbols(i).set_size(M(i));
974 bitmap(i) = graycode(k(i));
975 double average_energy = (M(i) * M(i) - 1) / 3.0;
976 double scaling_factor = std::sqrt(average_energy);
977
978 for(int j = 0; j < M(i); ++j) {
979 symbols(i)(j) = ((M(i) - 1) - j * 2) / scaling_factor;
980 bits2symbols(i)(bin2dec(bitmap(i).get_row(j))) = j;
981 }
982
983 // the "symbols" vector must end with a zero; only for a trick
984 // exploited in update_norm()
985 symbols(i)(M(i)) = 0.0;
986
987 spacing(i) = 2.0 / scaling_factor;
988 }
989}
990
991int ND_UPAM::sphere_search_SE(const vec &y_in, const mat &H,
992 const imat &zrange, double r, ivec &zhat)
993{
994 // The implementation of this function basically follows the
995 // Schnorr-Eucner algorithm described in Agrell et al. (IEEE
996 // Trans. IT, 2002), but taking into account constellation
997 // boundaries, see the "accelerated sphere decoder" in Boutros et
998 // al. (IEEE Globecom, 2003). No lattice reduction is performed.
999 // Potentially the function can be speeded up by performing
1000 // lattice reduction, but it seems difficult to keep track of
1001 // constellation boundaries.
1002
1003 mat R = chol(H.transpose() * H);
1004 mat Ri = inv(R);
1005 mat Q = H * Ri;
1006 vec y = Q.transpose() * y_in;
1007 mat Vi = Ri.transpose();
1008
1009 int n = H.cols();
1010 vec dist(n);
1011 dist(n - 1) = 0;
1012 double bestdist = r * r;
1013 int status = -1; // search failed
1014
1015 mat E = zeros(n, n);
1016 for(int i = 0; i < n; i++) { // E(k,:) = y*Vi;
1017 for(int j = 0; j < n; j++) {
1018 E(i * n + n - 1) += y(j) * Vi(j + n * i);
1019 }
1020 }
1021
1022 ivec z(n);
1023 zhat.set_size(n);
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);
1028 ivec step(n);
1029 step(n - 1) = sign_nozero_i(p);
1030
1031 // Run search loop
1032 int k = n - 1; // k uses natural indexing, goes from 0 to n-1
1033
1034 while(true) {
1035 double newdist = dist(k) + p * p;
1036
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);
1040 }
1041
1042 k--;
1043 dist(k) = newdist;
1044 z(k) = floor_i(0.5 + E(k + k * 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);
1048
1049 step(k) = sign_nozero_i(p);
1050 }
1051 else {
1052 while(true) {
1053 if(newdist < bestdist) {
1054 zhat = z;
1055 bestdist = newdist;
1056 status = 0;
1057 }
1058 else if(k == n - 1) {
1059 goto exit_point;
1060 }
1061 else {
1062 k++;
1063 }
1064
1065 z(k) += step(k);
1066
1067 if((z(k) < zrange(k, 0)) || (z(k) > zrange(k, 1))) {
1068 step(k) = (-step(k) - sign_nozero_i(step(k)));
1069 z(k) += step(k);
1070 }
1071
1072 if((z(k) >= zrange(k, 0)) && (z(k) <= zrange(k, 1))) {
1073 break;
1074 }
1075 }
1076
1077 p = (E(k + k * n) - z(k)) / Vi(k + k * n);
1078 step(k) = (-step(k) - sign_nozero_i(step(k)));
1079 }
1080 }
1081
1082exit_point:
1083 return status;
1084}
1085
1086
1087int ND_UPAM::sphere_decoding(const vec &y, const mat &H, double rstart,
1088 double rmax, double stepup,
1089 QLLRvec &detected_bits)
1090{
1091 it_assert(H.rows() == length(y),
1092 "ND_UPAM::sphere_decoding(): dimension mismatch");
1093 it_assert(H.cols() == nt,
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");
1097
1098 // This function can be improved, e.g., by using an ordered search.
1099
1100 vec ytemp = y;
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);
1105 }
1106
1107 imat crange(nt, 2);
1108 for(int i = 0; i < nt; i++) {
1109 crange(i, 0) = 0;
1110 crange(i, 1) = M(i) - 1;
1111 }
1112
1113 int status = 0;
1114 double r = rstart;
1115 ivec s(sum(M));
1116 while(r <= rmax) {
1117 status = sphere_search_SE(ytemp, Htemp, crange, r, s);
1118
1119 if(status == 0) { // search successful
1120 detected_bits.set_size(sum(k));
1121 int b = 0;
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;
1126 }
1127 else {
1128 detected_bits(b) = -1000;
1129 }
1130 b++;
1131 }
1132 }
1133
1134 return status;
1135 }
1136 r = r * stepup;
1137 }
1138
1139 return status;
1140}
1141
1142// ----------------------------------------------------------------------
1143// ND_UQAM
1144// ----------------------------------------------------------------------
1145
1146// The ND_UQAM (MIMO with uniform QAM) class could alternatively
1147// have been implemented by using a ND_UPAM class of twice the
1148// dimension, but this does not fit as elegantly into the class
1149// structure
1150
1151ND_UQAM::ND_UQAM(int nt, int Mary)
1152{
1153 set_M(nt, Mary);
1154}
1155
1156void ND_UQAM::set_M(int nt_in, int Mary)
1157{
1158 nt = nt_in;
1159 ivec Mary_temp(nt);
1160 Mary_temp = Mary;
1161 set_M(nt, Mary_temp);
1162}
1163
1164void ND_UQAM::set_M(int nt_in, ivec Mary)
1165{
1166 nt = nt_in;
1167 it_assert(length(Mary) == nt, "ND_UQAM::set_M(): Mary has wrong length");
1168 k.set_size(nt);
1169 M = Mary;
1170 L.set_size(nt);
1174
1175 for(int i = 0; i < nt; ++i) {
1176 k(i) = round_i(::log2(static_cast<double>(M(i))));
1177 it_assert((k(i) > 0) && ((1 << k(i)) == M(i)),
1178 "ND_UQAM::set_M(): M is not a power of 2");
1179
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");
1182
1183 symbols(i).set_size(M(i) + 1);
1184 bitmap(i).set_size(M(i), k(i));
1185 bits2symbols(i).set_size(M(i));
1186 double average_energy = (M(i) - 1) * 2.0 / 3.0;
1187 double scaling_factor = std::sqrt(average_energy);
1189
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);
1195 bitmap(i).set_row(j1 * L(i) + j2, concat(gray_code.get_row(j1),
1196 gray_code.get_row(j2)));
1197 bits2symbols(i)(bin2dec(bitmap(i).get_row(j1 * L(i) + j2)))
1198 = j1 * L(i) + j2;
1199 }
1200 }
1201
1202 // must end with a zero; only for a trick exploited in
1203 // update_norm()
1204 symbols(i)(M(i)) = 0.0;
1205 }
1206}
1207
1208void ND_UQAM::set_constellation_points(const int nth, const cvec& inConstellation, const ivec& in_bit2symbols)
1209{
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");
1213 it_assert(is_even(inConstellation.size()) && (inConstellation.size() > 0),
1214 "ND_UQAM::set_constellation_points(): Number of symbols needs to be even and non-zero");
1215
1216 symbols(nth).replace_mid(0, inConstellation);
1217
1218 bits2symbols(nth) = in_bit2symbols;
1219
1220 for(int m = 0; m < M(nth); ++m) {
1221 bitmap(nth).set_row(bits2symbols(nth)(m), dec2bin(k(nth), m));
1222 }
1223
1224 // must end with a zero; only for a trick exploited in
1225 // update_norm()
1226 symbols(nth)(M(nth)) = 0.0;
1227};
1228
1229// ----------------------------------------------------------------------
1230// ND_UPSK
1231// ----------------------------------------------------------------------
1232
1233ND_UPSK::ND_UPSK(int nt, int Mary)
1234{
1235 set_M(nt, Mary);
1236}
1237
1238void ND_UPSK::set_M(int nt_in, int Mary)
1239{
1240 nt = nt_in;
1241 ivec Mary_temp(nt);
1242 Mary_temp = Mary;
1243 set_M(nt, Mary_temp);
1244}
1245
1246void ND_UPSK::set_M(int nt_in, ivec Mary)
1247{
1248 nt = nt_in;
1249 it_assert(length(Mary) == nt, "ND_UPSK::set_M() Mary has wrong length");
1250 k.set_size(nt);
1251 M = Mary;
1255
1256 for(int i = 0; i < nt; ++i) {
1257 k(i) = round_i(::log2(static_cast<double>(M(i))));
1258 it_assert((k(i) > 0) && ((1 << k(i)) == M(i)),
1259 "ND_UPSK::set_M(): M is not a power of 2");
1260
1261 symbols(i).set_size(M(i) + 1);
1262 bits2symbols(i).set_size(M(i));
1263 bitmap(i) = graycode(k(i));
1264
1265 double delta = 2.0 * pi / M(i);
1266 double epsilon = delta / 10000.0;
1267
1268 for(int j = 0; j < M(i); ++j) {
1269 std::complex<double> symb
1270 = std::complex<double>(std::polar(1.0, delta * j));
1271
1272 if(std::abs(std::real(symb)) < epsilon) {
1273 symbols(i)(j) = std::complex<double>(0.0, std::imag(symb));
1274 }
1275 else if(std::abs(std::imag(symb)) < epsilon) {
1276 symbols(i)(j) = std::complex<double>(std::real(symb), 0.0);
1277 }
1278 else {
1279 symbols(i)(j) = symb;
1280 }
1281
1282 bits2symbols(i)(bin2dec(bitmap(i).get_row(j))) = j;
1283 }
1284
1285 // must end with a zero; only for a trick exploited in
1286 // update_norm()
1287 symbols(i)(M(i)) = 0.0;
1288 }
1289}
1290
1291} // namespace itpp
Definitions of Cholesky factorisation functions.
General array class.
Definition factory.h:40
Array< T > left(int n) const
Get n left elements of the array.
Definition array.h:357
void set_size(int n, bool copy=false)
Resizing an Array<T>.
Definition array.h:257
int length() const
Returns the number of data elements in the array object.
Definition array.h:157
QLLR jaclog(QLLR a, QLLR b) const
Jacobian logarithm.
Definition llr.h:280
QLLR to_qllr(double l) const
Convert a "real" LLR value to an LLR type.
Definition llr.h:245
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...
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.
Base class for N-dimensional vector (MIMO) channel modulators/demodulators with real-valued component...
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.
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.
Definition itassert.h:126
#define it_assert_debug(t, s)
Abort if t is not true and NDEBUG is not defined.
Definition itassert.h:107
#define it_assert(t, s)
Abort if t is not true.
Definition itassert.h:94
bool inv(const mat &X, mat &Y)
Inverse of real square matrix.
Definition inv.cpp:87
int pow2i(int x)
Calculate two to the power of x (2^x); x is integer.
Definition log_exp.h:53
vec log2(const vec &x)
log-2 of the elements
Definition log_exp.cpp:36
int levels2bits(int n)
Calculate the number of bits needed to represent n different values (levels).
Definition log_exp.h:92
Vec< T > cumsum(const Vec< T > &v)
Cumulative sum of all elements in the vector.
Definition matfunc.h:157
T sum(const Vec< T > &v)
Sum of all elements in the vector.
Definition matfunc.h:59
int length(const Vec< T > &v)
Length of vector.
Definition matfunc.h:51
T prod(const Vec< T > &v)
The product of all elements in the vector.
Definition matfunc.h:195
bool chol(const mat &X, mat &F)
Cholesky factorisation of real symmetric and positive definite matrix.
Definition cholesky.cpp:86
bmat graycode(int m)
Generate Gray code of blocklength m.
Definition commfunc.cpp:39
int mod(int k, int n)
Calculates the modulus, i.e. the signed reminder after division.
Definition elem_math.h:166
bool is_even(int x)
Return true if x is an even integer.
Definition misc.h:122
vec real(const cvec &data)
Real part of complex values.
vec sqr(const cvec &data)
Absolute square of elements.
Definition elem_math.cpp:36
cvec conj(const cvec &x)
Conjugate of complex value.
Definition elem_math.h:226
Vec< T > reverse(const Vec< T > &in)
Reverse the input vector.
Definition matfunc.h:777
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))
Definition misc_stat.cpp:77
Definitions of matrix inversion routines.
IT++ compatibility types and functions.
Logarithmic and exponenential functions - header file.
Mat< bin > bmat
bin matrix
Definition mat.h:508
Various functions on vectors and matrices - header file.
Miscellaneous statistics functions and classes - header file.
Definition of vector (MIMO) modulator classes.
itpp namespace
Definition itmex.h:37
ITPP_EXPORT int round_i(double x)
Round to nearest integer.
std::ostream & operator<<(std::ostream &output, const bin &inbin)
Output stream of bin.
Definition binary.cpp:36
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.
Definition array.h:486
int gray_code(int x)
Convert to Gray Code.
Definition converters.h:421
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.
Definition misc.h:103
int floor_i(double x)
The nearest smaller integer.
Definition converters.h:350
int abs(const itpp::bin &inbin)
absolute value of bin
Definition binary.h:186
Sorting functions.

Generated on Tue Aug 17 2021 10:59:15 for IT++ by Doxygen 1.12.0