IT++ Logo
modulator_nd.cpp
Go to the documentation of this file.
1 
29 #include <itpp/comm/modulator_nd.h>
30 #include <itpp/comm/commfunc.h>
32 #include <itpp/base/algebra/inv.h>
33 #include <itpp/base/matfunc.h>
35 #include <itpp/base/math/log_exp.h>
36 #include <itpp/base/converters.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 
45 namespace 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 
94 void Modulator_ND::update_LLR(const Array<QLLRvec> &logP_apriori,
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 
222 void 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 
250 void 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 
294 void 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 
320 void 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 
362 void 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 
463 void 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 
478 vec Modulator_NRD::modulate_bits(const bvec &bits) const
479 {
480  vec result(nt);
481  modulate_bits(bits, result);
482  return result;
483 }
484 
485 
486 void 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);
497  hspacings.set_size(nt);
498  yspacings.set_size(nt);
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 
525 void 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 
560 void 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 
591 QLLRvec 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 
601 void 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 
632 void 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 
653 void 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 
679 std::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 
709 void 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 
724 cvec Modulator_NCD::modulate_bits(const bvec &bits) const
725 {
726  cvec result(nt);
727  modulate_bits(bits, result);
728  return result;
729 }
730 
731 void 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);
742  hspacings.set_size(nt);
743  yspacings.set_size(nt);
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 
768 void 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 
801 void 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 
831 QLLRvec 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 
842 void 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 
874 void 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 
895 void 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 
923 std::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 
943 ND_UPAM::ND_UPAM(int nt, int Mary)
944 {
945  set_M(nt, Mary);
946 }
947 
948 void 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 
956 void 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;
962  bitmap.set_size(nt);
963  symbols.set_size(nt);
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 
991 int 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 
1082 exit_point:
1083  return status;
1084 }
1085 
1086 
1087 int 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 
1151 ND_UQAM::ND_UQAM(int nt, int Mary)
1152 {
1153  set_M(nt, Mary);
1154 }
1155 
1156 void 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 
1164 void 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);
1171  bitmap.set_size(nt);
1172  symbols.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 
1208 void 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 
1233 ND_UPSK::ND_UPSK(int nt, int Mary)
1234 {
1235  set_M(nt, Mary);
1236 }
1237 
1238 void 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 
1246 void 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;
1252  bitmap.set_size(nt);
1253  symbols.set_size(nt);
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
SourceForge Logo

Generated on Sat Jul 6 2013 10:54:23 for IT++ by Doxygen 1.8.2