IT++ Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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>
34 #include <itpp/base/converters.h>
35 #include <itpp/base/itcompat.h>
36 
37 namespace itpp
38 {
39 
40 // ----------------------------------------------------------------------
41 // Modulator_ND
42 // ----------------------------------------------------------------------
43 
45 {
46  QLLRvec result(2);
47 
48  if (l < 0) { // this can be done more efficiently
49  result(1) = -llrcalc.jaclog(0, -l);
50  result(0) = result(1) - l;
51  }
52  else {
53  result(0) = -llrcalc.jaclog(0, l);
54  result(1) = result(0) + l;
55  }
56  return result;
57 }
58 
60 {
61  Array<QLLRvec> result(length(l));
62  for (int i = 0; i < length(l); i++) {
63  result(i) = probabilities(l(i));
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 // ----------------------------------------------------------------------
125 // Modulator_NRD
126 // ----------------------------------------------------------------------
127 
128 void Modulator_NRD::modulate_bits(const bvec &bits, vec &out_symbols) const
129 {
130  it_assert(length(bits) == sum(k), "Modulator_NRD::modulate_bits(): "
131  "The number of input bits does not match.");
132 
133  out_symbols.set_size(nt);
134 
135  int b = 0;
136  for (int i = 0; i < nt; ++i) {
137  int symb = bin2dec(bits.mid(b, k(i)));
138  out_symbols(i) = symbols(i)(bits2symbols(i)(symb));
139  b += k(i);
140  }
141 }
142 
143 vec Modulator_NRD::modulate_bits(const bvec &bits) const
144 {
145  vec result(nt);
146  modulate_bits(bits, result);
147  return result;
148 }
149 
150 
151 void Modulator_NRD::demodulate_soft_bits(const vec &y, const mat &H,
152  double sigma2,
153  const QLLRvec &LLR_apriori,
154  QLLRvec &LLR_aposteriori,
155  Soft_Demod_Method method)
156 {
157  switch (method) {
158  case FULL_ENUM_LOGMAP:
159  demodulate_soft_bits(y, H, sigma2, LLR_apriori, LLR_aposteriori);
160  break;
161  case ZF_LOGMAP: {
162  it_assert(H.rows() >= H.cols(), "Modulator_NRD::demodulate_soft_bits():"
163  " ZF demodulation impossible for undetermined systems");
164  // Set up the ZF detector
165  mat Ht = H.T();
166  mat inv_HtH = inv(Ht * H);
167  vec shat = inv_HtH * Ht * y;
168  vec h = ones(shat.size());
169  for (int i = 0; i < shat.size(); ++i) {
170  // noise covariance of shat
171  double sigma_zf = std::sqrt(inv_HtH(i, i) * sigma2);
172  shat(i) /= sigma_zf;
173  h(i) /= sigma_zf;
174  }
175  demodulate_soft_bits(shat, h, 1.0, zeros_i(sum(k)), LLR_aposteriori);
176  }
177  break;
178  default:
179  it_error("Modulator_NRD::demodulate_soft_bits(): Improper soft "
180  "demodulation method");
181  }
182 }
183 
184 QLLRvec Modulator_NRD::demodulate_soft_bits(const vec &y, const mat &H,
185  double sigma2,
186  const QLLRvec &LLR_apriori,
187  Soft_Demod_Method method)
188 {
189  QLLRvec result;
190  demodulate_soft_bits(y, H, sigma2, LLR_apriori, result, method);
191  return result;
192 }
193 
194 void Modulator_NRD::demodulate_soft_bits(const vec &y, const vec &h,
195  double sigma2,
196  const QLLRvec &LLR_apriori,
197  QLLRvec &LLR_aposteriori)
198 {
199  it_assert(length(LLR_apriori) == sum(k),
200  "Modulator_NRD::demodulate_soft_bits(): Wrong sizes");
201  it_assert((length(h) == length(y)) && (length(h) == nt),
202  "Modulator_NRD::demodulate_soft_bits(): Wrong sizes");
203 
204  // set size of the output vector
205  LLR_aposteriori.set_size(LLR_apriori.size());
206 
207  // normalisation constant "minus one over two sigma^2"
208  double moo2s2 = -1.0 / (2.0 * sigma2);
209 
210  int b = 0;
211  for (int i = 0; i < nt; ++i) {
212  QLLRvec bnum = -QLLR_MAX * ones_i(k(i));
213  QLLRvec bdenom = bnum;
214  Array<QLLRvec> logP_apriori = probabilities(LLR_apriori(b, b + k(i) - 1));
215  for (int j = 0; j < M(i); ++j) {
216  double norm2 = moo2s2 * sqr(y(i) - h(i) * symbols(i)(j));
217  QLLR scaled_norm = llrcalc.to_qllr(norm2);
218  update_LLR(logP_apriori, j, scaled_norm, i, bnum, bdenom);
219  }
220  LLR_aposteriori.set_subvector(b, bnum - bdenom);
221  b += k(i);
222  }
223 }
224 
225 void Modulator_NRD::demodulate_soft_bits(const vec &y, const mat &H,
226  double sigma2,
227  const QLLRvec &LLR_apriori,
228  QLLRvec &LLR_aposteriori)
229 {
230  int np = sum(k); // number of bits in total
231  int nr = H.rows();
232  it_assert(length(LLR_apriori) == np,
233  "Modulator_NRD::demodulate_soft_bits(): Wrong sizes");
234  it_assert((H.rows() == length(y)) && (H.cols() == nt),
235  "Modulator_NRD::demodulate_soft_bits(): Wrong sizes");
236 
237  // set size of the output vector
238  LLR_aposteriori.set_size(LLR_apriori.size());
239 
240  // normalisation constant "minus one over two sigma^2"
241  double moo2s2 = -1.0 / (2.0 * sigma2);
242 
243  bool diff_update = false;
244  for (int i = 0; i < length(M); ++i) {
245  // differential update only pays off for larger dimensions
246  if (nt * M(i) > 4) {
247  diff_update = true;
248  break;
249  }
250  }
251 
252  Array<QLLRvec> logP_apriori = probabilities(LLR_apriori);
253 
254  mat Ht = H.T();
255  mat HtH = Ht * H;
256  vec ytH = Ht * y;
257 
258  QLLRvec bnum = -QLLR_MAX * ones_i(np);
259  QLLRvec bdenom = bnum;
260  ivec s = zeros_i(nt);
261  double norm = 0.0;
262 
263  // Go over all constellation points (r=dimension, s=vector of symbols)
264  int r = nt - 1;
265  while (true) {
266 
267  if (diff_update) {
268  update_norm(norm, r, s(r), 0, ytH, HtH, s);
269  }
270  s(r) = 0;
271 
272  while (true) {
273  if (s(r) > M(r) - 1) {
274  if (r == nt - 1) {
275  goto exit_point;
276  }
277  r++;
278  }
279  else {
280  if (r == 0) {
281  if (!diff_update) {
282  norm = 0.0;
283  for (int p = 0; p < nr; ++p) {
284  double d = y(p);
285  for (int i = 0; i < nt; ++i) {
286  d -= H(p, i) * symbols(i)[s[i]];
287  }
288  norm += sqr(d);
289  }
290  }
291  QLLR scaled_norm = llrcalc.to_qllr(norm * moo2s2);
292  update_LLR(logP_apriori, s, scaled_norm, bnum, bdenom);
293  }
294  else {
295  r--;
296  break;
297  }
298  }
299  if (diff_update) {
300  update_norm(norm, r, s(r), s(r) + 1, ytH, HtH, s);
301  }
302  s(r)++;
303  }
304  }
305 
306 exit_point:
307  LLR_aposteriori = bnum - bdenom;
308 
309 }
310 
311 
312 void Modulator_NRD::update_norm(double &norm, int k, int sold, int snew,
313  const vec &ytH, const mat &HtH,
314  const ivec &s)
315 {
316 
317  int m = length(s);
318  double cdiff = symbols(k)[snew] - symbols(k)[sold];
319 
320  norm += sqr(cdiff) * HtH(k, k);
321  cdiff *= 2.0;
322  norm -= cdiff * ytH[k];
323  for (int i = 0; i < m; ++i) {
324  norm += cdiff * HtH(i, k) * symbols(k)[s[i]];
325  }
326 }
327 
328 
329 std::ostream &operator<<(std::ostream &os, const Modulator_NRD &mod)
330 {
331  os << "--- REAL MIMO (NRD) CHANNEL ---------" << std::endl;
332  os << "Dimension (nt): " << mod.nt << std::endl;
333  os << "Bits per dimension (k): " << mod.k << std::endl;
334  os << "Symbols per dimension (M):" << mod.M << std::endl;
335  for (int i = 0; i < mod.nt; i++) {
336  os << "Bitmap for dimension " << i << ": " << mod.bitmap(i) << std::endl;
337  // skip printing the trailing zero
338  os << "Symbol coordinates for dimension " << i << ": " << mod.symbols(i).left(mod.M(i)) << std::endl;
339  }
340  os << mod.get_llrcalc() << std::endl;
341  return os;
342 }
343 
344 // ----------------------------------------------------------------------
345 // Modulator_NCD
346 // ----------------------------------------------------------------------
347 
348 void Modulator_NCD::modulate_bits(const bvec &bits, cvec &out_symbols) const
349 {
350  it_assert(length(bits) == sum(k), "Modulator_NCD::modulate_bits(): "
351  "The number of input bits does not match.");
352 
353  out_symbols.set_size(nt);
354 
355  int b = 0;
356  for (int i = 0; i < nt; ++i) {
357  int symb = bin2dec(bits.mid(b, k(i)));
358  out_symbols(i) = symbols(i)(bits2symbols(i)(symb));
359  b += k(i);
360  }
361 }
362 
363 cvec Modulator_NCD::modulate_bits(const bvec &bits) const
364 {
365  cvec result(nt);
366  modulate_bits(bits, result);
367  return result;
368 }
369 
370 
371 void Modulator_NCD::demodulate_soft_bits(const cvec &y, const cmat &H,
372  double sigma2,
373  const QLLRvec &LLR_apriori,
374  QLLRvec &LLR_aposteriori,
375  Soft_Demod_Method method)
376 {
377  switch (method) {
378  case FULL_ENUM_LOGMAP:
379  demodulate_soft_bits(y, H, sigma2, LLR_apriori, LLR_aposteriori);
380  break;
381  case ZF_LOGMAP: {
382  it_assert(H.rows() >= H.cols(), "Modulator_NCD::demodulate_soft_bits():"
383  " ZF demodulation impossible for undetermined systems");
384  // Set up the ZF detector
385  cmat Hht = H.H();
386  cmat inv_HhtH = inv(Hht * H);
387  cvec shat = inv_HhtH * Hht * y;
388  cvec h = ones_c(shat.size());
389  for (int i = 0; i < shat.size(); ++i) {
390  double sigma_zf = std::sqrt(real(inv_HhtH(i, i)) * sigma2);
391  shat(i) /= sigma_zf;
392  h(i) /= sigma_zf;
393  }
394  demodulate_soft_bits(shat, h, 1.0, zeros_i(sum(k)), LLR_aposteriori);
395  }
396  break;
397  default:
398  it_error("Modulator_NCD::demodulate_soft_bits(): Improper soft "
399  "demodulation method");
400  }
401 }
402 
403 QLLRvec Modulator_NCD::demodulate_soft_bits(const cvec &y, const cmat &H,
404  double sigma2,
405  const QLLRvec &LLR_apriori,
406  Soft_Demod_Method method)
407 {
408  QLLRvec result;
409  demodulate_soft_bits(y, H, sigma2, LLR_apriori, result, method);
410  return result;
411 }
412 
413 
414 void Modulator_NCD::demodulate_soft_bits(const cvec &y, const cvec &h,
415  double sigma2,
416  const QLLRvec &LLR_apriori,
417  QLLRvec &LLR_aposteriori)
418 {
419  it_assert(length(LLR_apriori) == sum(k),
420  "Modulator_NCD::demodulate_soft_bits(): Wrong sizes");
421  it_assert((length(h) == length(y)) && (length(h) == nt),
422  "Modulator_NCD::demodulate_soft_bits(): Wrong sizes");
423 
424  // set size of the output vector
425  LLR_aposteriori.set_size(LLR_apriori.size());
426 
427  // normalisation constant "minus one over sigma^2"
428  double moos2 = -1.0 / sigma2;
429 
430  int b = 0;
431  for (int i = 0; i < nt; ++i) {
432  QLLRvec bnum = -QLLR_MAX * ones_i(k(i));
433  QLLRvec bdenom = -QLLR_MAX * ones_i(k(i));
434  Array<QLLRvec> logP_apriori = probabilities(LLR_apriori(b, b + k(i) - 1));
435  for (int j = 0; j < M(i); ++j) {
436  double norm2 = moos2 * sqr(y(i) - h(i) * symbols(i)(j));
437  QLLR scaled_norm = llrcalc.to_qllr(norm2);
438  update_LLR(logP_apriori, j, scaled_norm, i, bnum, bdenom);
439  }
440  LLR_aposteriori.set_subvector(b, bnum - bdenom);
441  b += k(i);
442  }
443 }
444 
445 void Modulator_NCD::demodulate_soft_bits(const cvec &y, const cmat &H,
446  double sigma2,
447  const QLLRvec &LLR_apriori,
448  QLLRvec &LLR_aposteriori)
449 {
450  int np = sum(k); // number of bits in total
451  int nr = H.rows();
452  it_assert(length(LLR_apriori) == np,
453  "Modulator_NRD::demodulate_soft_bits(): Wrong sizes");
454  it_assert((H.rows() == length(y)) && (H.cols() == nt),
455  "Modulator_NRD::demodulate_soft_bits(): Wrong sizes");
456 
457  // set size of the output vector
458  LLR_aposteriori.set_size(LLR_apriori.size());
459 
460  // normalisation constant "minus one over sigma^2"
461  double moos2 = -1.0 / sigma2;
462 
463  bool diff_update = false;
464  for (int i = 0; i < length(M); ++i) {
465  // differential update only pays off for larger dimensions
466  if (nt * M(i) > 4) {
467  diff_update = true;
468  }
469  }
470 
471  Array<QLLRvec> logP_apriori = probabilities(LLR_apriori);
472 
473  cmat HtH = H.hermitian_transpose() * H;
474  cvec ytH = conj(H.hermitian_transpose() * y);
475 
476  QLLRvec bnum = -QLLR_MAX * ones_i(np);
477  QLLRvec bdenom = -QLLR_MAX * ones_i(np);
478  ivec s(nt);
479  s.clear();
480  double norm = 0.0;
481 
482  // Go over all constellation points (r=dimension, s=vector of symbols)
483  int r = nt - 1;
484  while (true) {
485 
486  if (diff_update) {
487  update_norm(norm, r, s(r), 0, ytH, HtH, s);
488  }
489  s(r) = 0;
490 
491  while (true) {
492  if (s(r) > M(r) - 1) {
493  if (r == nt - 1) {
494  goto exit_point;
495  }
496  r++;
497  }
498  else {
499  if (r == 0) {
500  if (!diff_update) {
501  norm = 0.0;
502  for (int p = 0; p < nr; ++p) {
503  std::complex<double> d = y(p);
504  for (int i = 0; i < nt; ++i) {
505  d -= H(p, i) * symbols(i)[s[i]];
506  }
507  norm += sqr(d);
508  }
509  }
510  QLLR scaled_norm = llrcalc.to_qllr(norm * moos2);
511  update_LLR(logP_apriori, s, scaled_norm, bnum, bdenom);
512  }
513  else {
514  r--;
515  break;
516  }
517  }
518  if (diff_update) {
519  update_norm(norm, r, s(r), s(r) + 1, ytH, HtH, s);
520  }
521  s(r)++;
522  }
523  }
524 
525 exit_point:
526  LLR_aposteriori = bnum - bdenom;
527 }
528 
529 
530 void Modulator_NCD::update_norm(double &norm, int k, int sold, int snew,
531  const cvec &ytH, const cmat &HtH,
532  const ivec &s)
533 {
534  int m = length(s);
535  std::complex<double> cdiff = symbols(k)[snew] - symbols(k)[sold];
536 
537  norm += sqr(cdiff) * (HtH(k, k).real());
538  cdiff *= 2.0;
539  norm -= (cdiff.real() * ytH[k].real() - cdiff.imag() * ytH[k].imag());
540  for (int i = 0; i < m; i++) {
541  norm += (cdiff * HtH(i, k) * conj(symbols(k)[s[i]])).real();
542  }
543 }
544 
545 
546 std::ostream &operator<<(std::ostream &os, const Modulator_NCD &mod)
547 {
548  os << "--- COMPLEX MIMO (NCD) CHANNEL --------" << std::endl;
549  os << "Dimension (nt): " << mod.nt << std::endl;
550  os << "Bits per dimension (k): " << mod.k << std::endl;
551  os << "Symbols per dimension (M):" << mod.M << std::endl;
552  for (int i = 0; i < mod.nt; i++) {
553  os << "Bitmap for dimension " << i << ": "
554  << mod.bitmap(i) << std::endl;
555  os << "Symbol coordinates for dimension " << i << ": "
556  << mod.symbols(i).left(mod.M(i)) << std::endl;
557  }
558  os << mod.get_llrcalc() << std::endl;
559  return os;
560 }
561 
562 // ----------------------------------------------------------------------
563 // ND_UPAM
564 // ----------------------------------------------------------------------
565 
566 ND_UPAM::ND_UPAM(int nt, int Mary)
567 {
568  set_M(nt, Mary);
569 }
570 
571 void ND_UPAM::set_M(int nt_in, int Mary)
572 {
573  nt = nt_in;
574  ivec Mary_temp(nt);
575  Mary_temp = Mary;
576  set_M(nt, Mary_temp);
577 }
578 
579 void ND_UPAM::set_M(int nt_in, ivec Mary)
580 {
581  nt = nt_in;
582  it_assert(length(Mary) == nt, "ND_UPAM::set_M(): Mary has wrong length");
583  k.set_size(nt);
584  M = Mary;
585  bitmap.set_size(nt);
588  spacing.set_size(nt);
589 
590  for (int i = 0; i < nt; i++) {
591  k(i) = round_i(::log2(static_cast<double>(M(i))));
592  it_assert((k(i) > 0) && ((1 << k(i)) == M(i)),
593  "ND_UPAM::set_M(): M is not a power of 2.");
594 
595  symbols(i).set_size(M(i) + 1);
596  bits2symbols(i).set_size(M(i));
597  bitmap(i) = graycode(k(i));
598  double average_energy = (M(i) * M(i) - 1) / 3.0;
599  double scaling_factor = std::sqrt(average_energy);
600 
601  for (int j = 0; j < M(i); ++j) {
602  symbols(i)(j) = ((M(i) - 1) - j * 2) / scaling_factor;
603  bits2symbols(i)(bin2dec(bitmap(i).get_row(j))) = j;
604  }
605 
606  // the "symbols" vector must end with a zero; only for a trick
607  // exploited in update_norm()
608  symbols(i)(M(i)) = 0.0;
609 
610  spacing(i) = 2.0 / scaling_factor;
611  }
612 }
613 
614 int ND_UPAM::sphere_search_SE(const vec &y_in, const mat &H,
615  const imat &zrange, double r, ivec &zhat)
616 {
617  // The implementation of this function basically follows the
618  // Schnorr-Eucner algorithm described in Agrell et al. (IEEE
619  // Trans. IT, 2002), but taking into account constellation
620  // boundaries, see the "accelerated sphere decoder" in Boutros et
621  // al. (IEEE Globecom, 2003). No lattice reduction is performed.
622  // Potentially the function can be speeded up by performing
623  // lattice reduction, but it seems difficult to keep track of
624  // constellation boundaries.
625 
626  mat R = chol(H.transpose() * H);
627  mat Ri = inv(R);
628  mat Q = H * Ri;
629  vec y = Q.transpose() * y_in;
630  mat Vi = Ri.transpose();
631 
632  int n = H.cols();
633  vec dist(n);
634  dist(n - 1) = 0;
635  double bestdist = r * r;
636  int status = -1; // search failed
637 
638  mat E = zeros(n, n);
639  for (int i = 0; i < n; i++) { // E(k,:) = y*Vi;
640  for (int j = 0; j < n; j++) {
641  E(i*n + n - 1) += y(j) * Vi(j + n * i);
642  }
643  }
644 
645  ivec z(n);
646  zhat.set_size(n);
647  z(n - 1) = floor_i(0.5 + E(n * n - 1));
648  z(n - 1) = std::max(z(n - 1), zrange(n - 1, 0));
649  z(n - 1) = std::min(z(n - 1), zrange(n - 1, 1));
650  double p = (E(n * n - 1) - z(n - 1)) / Vi(n * n - 1);
651  ivec step(n);
652  step(n - 1) = sign_nozero_i(p);
653 
654  // Run search loop
655  int k = n - 1; // k uses natural indexing, goes from 0 to n-1
656 
657  while (true) {
658  double newdist = dist(k) + p * p;
659 
660  if ((newdist < bestdist) && (k != 0)) {
661  for (int i = 0; i < k; i++) {
662  E(k - 1 + i*n) = E(k + i * n) - p * Vi(k + i * n);
663  }
664 
665  k--;
666  dist(k) = newdist;
667  z(k) = floor_i(0.5 + E(k + k * n));
668  z(k) = std::max(z(k), zrange(k, 0));
669  z(k) = std::min(z(k), zrange(k, 1));
670  p = (E(k + k * n) - z(k)) / Vi(k + k * n);
671 
672  step(k) = sign_nozero_i(p);
673  }
674  else {
675  while (true) {
676  if (newdist < bestdist) {
677  zhat = z;
678  bestdist = newdist;
679  status = 0;
680  }
681  else if (k == n - 1) {
682  goto exit_point;
683  }
684  else {
685  k++;
686  }
687 
688  z(k) += step(k);
689 
690  if ((z(k) < zrange(k, 0)) || (z(k) > zrange(k, 1))) {
691  step(k) = (-step(k) - sign_nozero_i(step(k)));
692  z(k) += step(k);
693  }
694 
695  if ((z(k) >= zrange(k, 0)) && (z(k) <= zrange(k, 1))) {
696  break;
697  }
698  }
699 
700  p = (E(k + k * n) - z(k)) / Vi(k + k * n);
701  step(k) = (-step(k) - sign_nozero_i(step(k)));
702  }
703  }
704 
705 exit_point:
706  return status;
707 }
708 
709 
710 int ND_UPAM::sphere_decoding(const vec &y, const mat &H, double rstart,
711  double rmax, double stepup,
712  QLLRvec &detected_bits)
713 {
714  it_assert(H.rows() == length(y),
715  "ND_UPAM::sphere_decoding(): dimension mismatch");
716  it_assert(H.cols() == nt,
717  "ND_UPAM::sphere_decoding(): dimension mismatch");
718  it_assert(rstart > 0, "ND_UPAM::sphere_decoding(): radius error");
719  it_assert(rmax > rstart, "ND_UPAM::sphere_decoding(): radius error");
720 
721  // This function can be improved, e.g., by using an ordered search.
722 
723  vec ytemp = y;
724  mat Htemp(H.rows(), H.cols());
725  for (int i = 0; i < H.cols(); i++) {
726  Htemp.set_col(i, H.get_col(i)*spacing(i));
727  ytemp += Htemp.get_col(i) * 0.5 * (M(i) - 1.0);
728  }
729 
730  imat crange(nt, 2);
731  for (int i = 0; i < nt; i++) {
732  crange(i, 0) = 0;
733  crange(i, 1) = M(i) - 1;
734  }
735 
736  int status = 0;
737  double r = rstart;
738  ivec s(sum(M));
739  while (r <= rmax) {
740  status = sphere_search_SE(ytemp, Htemp, crange, r, s);
741 
742  if (status == 0) { // search successful
743  detected_bits.set_size(sum(k));
744  int b = 0;
745  for (int j = 0; j < nt; j++) {
746  for (int i = 0; i < k(j); i++) {
747  if (bitmap(j)((M(j) - 1 - s[j]), i) == 0) {
748  detected_bits(b) = 1000;
749  }
750  else {
751  detected_bits(b) = -1000;
752  }
753  b++;
754  }
755  }
756 
757  return status;
758  }
759  r = r * stepup;
760  }
761 
762  return status;
763 }
764 
765 // ----------------------------------------------------------------------
766 // ND_UQAM
767 // ----------------------------------------------------------------------
768 
769 // The ND_UQAM (MIMO with uniform QAM) class could alternatively
770 // have been implemented by using a ND_UPAM class of twice the
771 // dimension, but this does not fit as elegantly into the class
772 // structure
773 
774 ND_UQAM::ND_UQAM(int nt, int Mary)
775 {
776  set_M(nt, Mary);
777 }
778 
779 void ND_UQAM::set_M(int nt_in, int Mary)
780 {
781  nt = nt_in;
782  ivec Mary_temp(nt);
783  Mary_temp = Mary;
784  set_M(nt, Mary_temp);
785 }
786 
787 void ND_UQAM::set_M(int nt_in, ivec Mary)
788 {
789  nt = nt_in;
790  it_assert(length(Mary) == nt, "ND_UQAM::set_M(): Mary has wrong length");
791  k.set_size(nt);
792  M = Mary;
793  L.set_size(nt);
794  bitmap.set_size(nt);
797 
798  for (int i = 0; i < nt; ++i) {
799  k(i) = round_i(::log2(static_cast<double>(M(i))));
800  it_assert((k(i) > 0) && ((1 << k(i)) == M(i)),
801  "ND_UQAM::set_M(): M is not a power of 2");
802 
803  L(i) = round_i(std::sqrt(static_cast<double>(M(i))));
804  it_assert(L(i)*L(i) == M(i), "ND_UQAM: constellation M must be square");
805 
806  symbols(i).set_size(M(i) + 1);
807  bitmap(i).set_size(M(i), k(i));
808  bits2symbols(i).set_size(M(i));
809  double average_energy = (M(i) - 1) * 2.0 / 3.0;
810  double scaling_factor = std::sqrt(average_energy);
812 
813  for (int j1 = 0; j1 < L(i); ++j1) {
814  for (int j2 = 0; j2 < L(i); ++j2) {
815  symbols(i)(j1*L(i) + j2) =
816  std::complex<double>(((L(i) - 1) - j2 * 2.0) / scaling_factor,
817  ((L(i) - 1) - j1 * 2.0) / scaling_factor);
818  bitmap(i).set_row(j1*L(i) + j2, concat(gray_code.get_row(j1),
819  gray_code.get_row(j2)));
820  bits2symbols(i)(bin2dec(bitmap(i).get_row(j1*L(i) + j2)))
821  = j1 * L(i) + j2;
822  }
823  }
824 
825  // must end with a zero; only for a trick exploited in
826  // update_norm()
827  symbols(i)(M(i)) = 0.0;
828  }
829 }
830 
831 // ----------------------------------------------------------------------
832 // ND_UPSK
833 // ----------------------------------------------------------------------
834 
835 ND_UPSK::ND_UPSK(int nt, int Mary)
836 {
837  set_M(nt, Mary);
838 }
839 
840 void ND_UPSK::set_M(int nt_in, int Mary)
841 {
842  nt = nt_in;
843  ivec Mary_temp(nt);
844  Mary_temp = Mary;
845  set_M(nt, Mary_temp);
846 }
847 
848 void ND_UPSK::set_M(int nt_in, ivec Mary)
849 {
850  nt = nt_in;
851  it_assert(length(Mary) == nt, "ND_UPSK::set_M() Mary has wrong length");
852  k.set_size(nt);
853  M = Mary;
854  bitmap.set_size(nt);
857 
858  for (int i = 0; i < nt; ++i) {
859  k(i) = round_i(::log2(static_cast<double>(M(i))));
860  it_assert((k(i) > 0) && ((1 << k(i)) == M(i)),
861  "ND_UPSK::set_M(): M is not a power of 2");
862 
863  symbols(i).set_size(M(i) + 1);
864  bits2symbols(i).set_size(M(i));
865  bitmap(i) = graycode(k(i));
866 
867  double delta = 2.0 * pi / M(i);
868  double epsilon = delta / 10000.0;
869 
870  for (int j = 0; j < M(i); ++j) {
871  std::complex<double> symb
872  = std::complex<double>(std::polar(1.0, delta * j));
873 
874  if (std::abs(std::real(symb)) < epsilon) {
875  symbols(i)(j) = std::complex<double>(0.0, std::imag(symb));
876  }
877  else if (std::abs(std::imag(symb)) < epsilon) {
878  symbols(i)(j) = std::complex<double>(std::real(symb), 0.0);
879  }
880  else {
881  symbols(i)(j) = symb;
882  }
883 
884  bits2symbols(i)(bin2dec(bitmap(i).get_row(j))) = j;
885  }
886 
887  // must end with a zero; only for a trick exploited in
888  // update_norm()
889  symbols(i)(M(i)) = 0.0;
890  }
891 }
892 
893 } // namespace itpp
SourceForge Logo

Generated on Fri Mar 21 2014 17:14:13 for IT++ by Doxygen 1.8.1.2