00001 00030 #include <itpp/comm/modulator_nd.h> 00031 #include <itpp/comm/commfunc.h> 00032 #include <itpp/base/algebra/cholesky.h> 00033 #include <itpp/base/algebra/inv.h> 00034 #include <itpp/base/math/elem_math.h> 00035 #include <itpp/base/converters.h> 00036 00037 namespace itpp 00038 { 00039 00040 // ---------------------------------------------------------------------- 00041 // Modulator_ND 00042 // ---------------------------------------------------------------------- 00043 00044 QLLRvec Modulator_ND::probabilities(QLLR l) 00045 { 00046 QLLRvec result(2); 00047 00048 if (l < 0) { // this can be done more efficiently 00049 result(1) = -llrcalc.jaclog(0, -l); 00050 result(0) = result(1) - l; 00051 } 00052 else { 00053 result(0) = -llrcalc.jaclog(0, l); 00054 result(1) = result(0) + l; 00055 } 00056 return result; 00057 } 00058 00059 Array<QLLRvec> Modulator_ND::probabilities(const QLLRvec &l) 00060 { 00061 Array<QLLRvec> result(length(l)); 00062 for (int i = 0; i < length(l); i++) { 00063 result(i) = probabilities(l(i)); 00064 } 00065 return result; 00066 } 00067 00068 void Modulator_ND::update_LLR(const Array<QLLRvec> &logP_apriori, int s, 00069 QLLR scaled_norm, int j, QLLRvec &p1, 00070 QLLRvec &p0) 00071 { 00072 QLLR log_apriori_prob_const_point = 0; 00073 int b = 0; 00074 for (int i = 0; i < k(j); i++) { 00075 log_apriori_prob_const_point += 00076 ((bitmap(j)(s, i) == 0) ? logP_apriori(b)(1) : logP_apriori(b)(0)); 00077 b++; 00078 } 00079 00080 b = 0; 00081 for (int i = 0; i < k(j); i++) { 00082 if (bitmap(j)(s, i) == 0) { 00083 p1(b) = llrcalc.jaclog(p1(b), scaled_norm 00084 + log_apriori_prob_const_point); 00085 } 00086 else { 00087 p0(b) = llrcalc.jaclog(p0(b), scaled_norm 00088 + log_apriori_prob_const_point); 00089 } 00090 b++; 00091 } 00092 } 00093 00094 void Modulator_ND::update_LLR(const Array<QLLRvec> &logP_apriori, 00095 const ivec &s, QLLR scaled_norm, 00096 QLLRvec &p1, QLLRvec &p0) 00097 { 00098 QLLR log_apriori_prob_const_point = 0; 00099 int b = 0; 00100 for (int j = 0; j < nt; j++) { 00101 for (int i = 0; i < k(j); i++) { 00102 log_apriori_prob_const_point += 00103 ((bitmap(j)(s[j], i) == 0) ? logP_apriori(b)(1) : logP_apriori(b)(0)); 00104 b++; 00105 } 00106 } 00107 00108 b = 0; 00109 for (int j = 0; j < nt; j++) { 00110 for (int i = 0; i < k(j); i++) { 00111 if (bitmap(j)(s[j], i) == 0) { 00112 p1(b) = llrcalc.jaclog(p1(b), scaled_norm 00113 + log_apriori_prob_const_point); 00114 } 00115 else { 00116 p0(b) = llrcalc.jaclog(p0(b), scaled_norm 00117 + log_apriori_prob_const_point); 00118 } 00119 b++; 00120 } 00121 } 00122 } 00123 00124 // ---------------------------------------------------------------------- 00125 // Modulator_NRD 00126 // ---------------------------------------------------------------------- 00127 00128 void Modulator_NRD::modulate_bits(const bvec &bits, vec &out_symbols) const 00129 { 00130 it_assert(length(bits) == sum(k), "Modulator_NRD::modulate_bits(): " 00131 "The number of input bits does not match."); 00132 00133 out_symbols.set_size(nt); 00134 00135 int b = 0; 00136 for (int i = 0; i < nt; ++i) { 00137 int symb = bin2dec(bits.mid(b, k(i))); 00138 out_symbols(i) = symbols(i)(bits2symbols(i)(symb)); 00139 b += k(i); 00140 } 00141 } 00142 00143 vec Modulator_NRD::modulate_bits(const bvec &bits) const 00144 { 00145 vec result(nt); 00146 modulate_bits(bits, result); 00147 return result; 00148 } 00149 00150 00151 void Modulator_NRD::demodulate_soft_bits(const vec &y, const mat &H, 00152 double sigma2, 00153 const QLLRvec &LLR_apriori, 00154 QLLRvec &LLR_aposteriori, 00155 Soft_Demod_Method method) 00156 { 00157 switch (method) { 00158 case FULL_ENUM_LOGMAP: 00159 demodulate_soft_bits(y, H, sigma2, LLR_apriori, LLR_aposteriori); 00160 break; 00161 case ZF_LOGMAP: { 00162 it_assert(H.rows() >= H.cols(), "Modulator_NRD::demodulate_soft_bits():" 00163 " ZF demodulation impossible for undetermined systems"); 00164 // Set up the ZF detector 00165 mat Ht = H.T(); 00166 mat inv_HtH = inv(Ht * H); 00167 vec shat = inv_HtH * Ht * y; 00168 vec h = ones(shat.size()); 00169 for (int i = 0; i < shat.size(); ++i) { 00170 // noise covariance of shat 00171 double sigma_zf = std::sqrt(inv_HtH(i, i) * sigma2); 00172 shat(i) /= sigma_zf; 00173 h(i) /= sigma_zf; 00174 } 00175 demodulate_soft_bits(shat, h, 1.0, zeros_i(sum(k)), LLR_aposteriori); 00176 } 00177 break; 00178 default: 00179 it_error("Modulator_NRD::demodulate_soft_bits(): Improper soft " 00180 "demodulation method"); 00181 } 00182 } 00183 00184 QLLRvec Modulator_NRD::demodulate_soft_bits(const vec &y, const mat &H, 00185 double sigma2, 00186 const QLLRvec &LLR_apriori, 00187 Soft_Demod_Method method) 00188 { 00189 QLLRvec result; 00190 demodulate_soft_bits(y, H, sigma2, LLR_apriori, result, method); 00191 return result; 00192 } 00193 00194 void Modulator_NRD::demodulate_soft_bits(const vec &y, const vec &h, 00195 double sigma2, 00196 const QLLRvec &LLR_apriori, 00197 QLLRvec &LLR_aposteriori) 00198 { 00199 it_assert(length(LLR_apriori) == sum(k), 00200 "Modulator_NRD::demodulate_soft_bits(): Wrong sizes"); 00201 it_assert((length(h) == length(y)) && (length(h) == nt), 00202 "Modulator_NRD::demodulate_soft_bits(): Wrong sizes"); 00203 00204 // set size of the output vector 00205 LLR_aposteriori.set_size(LLR_apriori.size()); 00206 00207 // normalisation constant "minus one over two sigma^2" 00208 double moo2s2 = -1.0 / (2.0 * sigma2); 00209 00210 int b = 0; 00211 for (int i = 0; i < nt; ++i) { 00212 QLLRvec bnum = -QLLR_MAX * ones_i(k(i)); 00213 QLLRvec bdenom = bnum; 00214 Array<QLLRvec> logP_apriori = probabilities(LLR_apriori(b, b + k(i) - 1)); 00215 for (int j = 0; j < M(i); ++j) { 00216 double norm2 = moo2s2 * sqr(y(i) - h(i) * symbols(i)(j)); 00217 QLLR scaled_norm = llrcalc.to_qllr(norm2); 00218 update_LLR(logP_apriori, j, scaled_norm, i, bnum, bdenom); 00219 } 00220 LLR_aposteriori.set_subvector(b, bnum - bdenom); 00221 b += k(i); 00222 } 00223 } 00224 00225 void Modulator_NRD::demodulate_soft_bits(const vec &y, const mat &H, 00226 double sigma2, 00227 const QLLRvec &LLR_apriori, 00228 QLLRvec &LLR_aposteriori) 00229 { 00230 int np = sum(k); // number of bits in total 00231 int nr = H.rows(); 00232 it_assert(length(LLR_apriori) == np, 00233 "Modulator_NRD::demodulate_soft_bits(): Wrong sizes"); 00234 it_assert((H.rows() == length(y)) && (H.cols() == nt), 00235 "Modulator_NRD::demodulate_soft_bits(): Wrong sizes"); 00236 00237 // set size of the output vector 00238 LLR_aposteriori.set_size(LLR_apriori.size()); 00239 00240 // normalisation constant "minus one over two sigma^2" 00241 double moo2s2 = -1.0 / (2.0 * sigma2); 00242 00243 bool diff_update = false; 00244 for (int i = 0; i < length(M); ++i) { 00245 // differential update only pays off for larger dimensions 00246 if (nt * M(i) > 4) { 00247 diff_update = true; 00248 break; 00249 } 00250 } 00251 00252 Array<QLLRvec> logP_apriori = probabilities(LLR_apriori); 00253 00254 mat Ht = H.T(); 00255 mat HtH = Ht * H; 00256 vec ytH = Ht * y; 00257 00258 QLLRvec bnum = -QLLR_MAX * ones_i(np); 00259 QLLRvec bdenom = bnum; 00260 ivec s = zeros_i(nt); 00261 double norm = 0.0; 00262 00263 // Go over all constellation points (r=dimension, s=vector of symbols) 00264 int r = nt - 1; 00265 while (true) { 00266 00267 if (diff_update) { 00268 update_norm(norm, r, s(r), 0, ytH, HtH, s); 00269 } 00270 s(r) = 0; 00271 00272 while (true) { 00273 if (s(r) > M(r) - 1) { 00274 if (r == nt - 1) { 00275 goto exit_point; 00276 } 00277 r++; 00278 } 00279 else { 00280 if (r == 0) { 00281 if (!diff_update) { 00282 norm = 0.0; 00283 for (int p = 0; p < nr; ++p) { 00284 double d = y(p); 00285 for (int i = 0; i < nt; ++i) { 00286 d -= H(p, i) * symbols(i)[s[i]]; 00287 } 00288 norm += sqr(d); 00289 } 00290 } 00291 QLLR scaled_norm = llrcalc.to_qllr(norm * moo2s2); 00292 update_LLR(logP_apriori, s, scaled_norm, bnum, bdenom); 00293 } 00294 else { 00295 r--; 00296 break; 00297 } 00298 } 00299 if (diff_update) { 00300 update_norm(norm, r, s(r), s(r) + 1, ytH, HtH, s); 00301 } 00302 s(r)++; 00303 } 00304 } 00305 00306 exit_point: 00307 LLR_aposteriori = bnum - bdenom; 00308 00309 } 00310 00311 00312 void Modulator_NRD::update_norm(double &norm, int k, int sold, int snew, 00313 const vec &ytH, const mat &HtH, 00314 const ivec &s) 00315 { 00316 00317 int m = length(s); 00318 double cdiff = symbols(k)[snew] - symbols(k)[sold]; 00319 00320 norm += sqr(cdiff) * HtH(k, k); 00321 cdiff *= 2.0; 00322 norm -= cdiff * ytH[k]; 00323 for (int i = 0; i < m; ++i) { 00324 norm += cdiff * HtH(i, k) * symbols(k)[s[i]]; 00325 } 00326 } 00327 00328 00329 std::ostream &operator<<(std::ostream &os, const Modulator_NRD &mod) 00330 { 00331 os << "--- REAL MIMO (NRD) CHANNEL ---------" << std::endl; 00332 os << "Dimension (nt): " << mod.nt << std::endl; 00333 os << "Bits per dimension (k): " << mod.k << std::endl; 00334 os << "Symbols per dimension (M):" << mod.M << std::endl; 00335 for (int i = 0; i < mod.nt; i++) { 00336 os << "Bitmap for dimension " << i << ": " << mod.bitmap(i) << std::endl; 00337 // skip printing the trailing zero 00338 os << "Symbol coordinates for dimension " << i << ": " << mod.symbols(i).left(mod.M(i)) << std::endl; 00339 } 00340 os << mod.get_llrcalc() << std::endl; 00341 return os; 00342 } 00343 00344 // ---------------------------------------------------------------------- 00345 // Modulator_NCD 00346 // ---------------------------------------------------------------------- 00347 00348 void Modulator_NCD::modulate_bits(const bvec &bits, cvec &out_symbols) const 00349 { 00350 it_assert(length(bits) == sum(k), "Modulator_NCD::modulate_bits(): " 00351 "The number of input bits does not match."); 00352 00353 out_symbols.set_size(nt); 00354 00355 int b = 0; 00356 for (int i = 0; i < nt; ++i) { 00357 int symb = bin2dec(bits.mid(b, k(i))); 00358 out_symbols(i) = symbols(i)(bits2symbols(i)(symb)); 00359 b += k(i); 00360 } 00361 } 00362 00363 cvec Modulator_NCD::modulate_bits(const bvec &bits) const 00364 { 00365 cvec result(nt); 00366 modulate_bits(bits, result); 00367 return result; 00368 } 00369 00370 00371 void Modulator_NCD::demodulate_soft_bits(const cvec &y, const cmat &H, 00372 double sigma2, 00373 const QLLRvec &LLR_apriori, 00374 QLLRvec &LLR_aposteriori, 00375 Soft_Demod_Method method) 00376 { 00377 switch (method) { 00378 case FULL_ENUM_LOGMAP: 00379 demodulate_soft_bits(y, H, sigma2, LLR_apriori, LLR_aposteriori); 00380 break; 00381 case ZF_LOGMAP: { 00382 it_assert(H.rows() >= H.cols(), "Modulator_NCD::demodulate_soft_bits():" 00383 " ZF demodulation impossible for undetermined systems"); 00384 // Set up the ZF detector 00385 cmat Hht = H.H(); 00386 cmat inv_HhtH = inv(Hht * H); 00387 cvec shat = inv_HhtH * Hht * y; 00388 cvec h = ones_c(shat.size()); 00389 for (int i = 0; i < shat.size(); ++i) { 00390 double sigma_zf = std::sqrt(real(inv_HhtH(i, i)) * sigma2); 00391 shat(i) /= sigma_zf; 00392 h(i) /= sigma_zf; 00393 } 00394 demodulate_soft_bits(shat, h, 1.0, zeros_i(sum(k)), LLR_aposteriori); 00395 } 00396 break; 00397 default: 00398 it_error("Modulator_NCD::demodulate_soft_bits(): Improper soft " 00399 "demodulation method"); 00400 } 00401 } 00402 00403 QLLRvec Modulator_NCD::demodulate_soft_bits(const cvec &y, const cmat &H, 00404 double sigma2, 00405 const QLLRvec &LLR_apriori, 00406 Soft_Demod_Method method) 00407 { 00408 QLLRvec result; 00409 demodulate_soft_bits(y, H, sigma2, LLR_apriori, result, method); 00410 return result; 00411 } 00412 00413 00414 void Modulator_NCD::demodulate_soft_bits(const cvec &y, const cvec &h, 00415 double sigma2, 00416 const QLLRvec &LLR_apriori, 00417 QLLRvec &LLR_aposteriori) 00418 { 00419 it_assert(length(LLR_apriori) == sum(k), 00420 "Modulator_NCD::demodulate_soft_bits(): Wrong sizes"); 00421 it_assert((length(h) == length(y)) && (length(h) == nt), 00422 "Modulator_NCD::demodulate_soft_bits(): Wrong sizes"); 00423 00424 // set size of the output vector 00425 LLR_aposteriori.set_size(LLR_apriori.size()); 00426 00427 // normalisation constant "minus one over sigma^2" 00428 double moos2 = -1.0 / sigma2; 00429 00430 int b = 0; 00431 for (int i = 0; i < nt; ++i) { 00432 QLLRvec bnum = -QLLR_MAX * ones_i(k(i)); 00433 QLLRvec bdenom = -QLLR_MAX * ones_i(k(i)); 00434 Array<QLLRvec> logP_apriori = probabilities(LLR_apriori(b, b + k(i) - 1)); 00435 for (int j = 0; j < M(i); ++j) { 00436 double norm2 = moos2 * sqr(y(i) - h(i) * symbols(i)(j)); 00437 QLLR scaled_norm = llrcalc.to_qllr(norm2); 00438 update_LLR(logP_apriori, j, scaled_norm, i, bnum, bdenom); 00439 } 00440 LLR_aposteriori.set_subvector(b, bnum - bdenom); 00441 b += k(i); 00442 } 00443 } 00444 00445 void Modulator_NCD::demodulate_soft_bits(const cvec &y, const cmat &H, 00446 double sigma2, 00447 const QLLRvec &LLR_apriori, 00448 QLLRvec &LLR_aposteriori) 00449 { 00450 int np = sum(k); // number of bits in total 00451 int nr = H.rows(); 00452 it_assert(length(LLR_apriori) == np, 00453 "Modulator_NRD::demodulate_soft_bits(): Wrong sizes"); 00454 it_assert((H.rows() == length(y)) && (H.cols() == nt), 00455 "Modulator_NRD::demodulate_soft_bits(): Wrong sizes"); 00456 00457 // set size of the output vector 00458 LLR_aposteriori.set_size(LLR_apriori.size()); 00459 00460 // normalisation constant "minus one over sigma^2" 00461 double moos2 = -1.0 / sigma2; 00462 00463 bool diff_update = false; 00464 for (int i = 0; i < length(M); ++i) { 00465 // differential update only pays off for larger dimensions 00466 if (nt * M(i) > 4) { 00467 diff_update = true; 00468 } 00469 } 00470 00471 Array<QLLRvec> logP_apriori = probabilities(LLR_apriori); 00472 00473 cmat HtH = H.hermitian_transpose() * H; 00474 cvec ytH = conj(H.hermitian_transpose() * y); 00475 00476 QLLRvec bnum = -QLLR_MAX * ones_i(np); 00477 QLLRvec bdenom = -QLLR_MAX * ones_i(np); 00478 ivec s(nt); 00479 s.clear(); 00480 double norm = 0.0; 00481 00482 // Go over all constellation points (r=dimension, s=vector of symbols) 00483 int r = nt - 1; 00484 while (true) { 00485 00486 if (diff_update) { 00487 update_norm(norm, r, s(r), 0, ytH, HtH, s); 00488 } 00489 s(r) = 0; 00490 00491 while (true) { 00492 if (s(r) > M(r) - 1) { 00493 if (r == nt - 1) { 00494 goto exit_point; 00495 } 00496 r++; 00497 } 00498 else { 00499 if (r == 0) { 00500 if (!diff_update) { 00501 norm = 0.0; 00502 for (int p = 0; p < nr; ++p) { 00503 std::complex<double> d = y(p); 00504 for (int i = 0; i < nt; ++i) { 00505 d -= H(p, i) * symbols(i)[s[i]]; 00506 } 00507 norm += sqr(d); 00508 } 00509 } 00510 QLLR scaled_norm = llrcalc.to_qllr(norm * moos2); 00511 update_LLR(logP_apriori, s, scaled_norm, bnum, bdenom); 00512 } 00513 else { 00514 r--; 00515 break; 00516 } 00517 } 00518 if (diff_update) { 00519 update_norm(norm, r, s(r), s(r) + 1, ytH, HtH, s); 00520 } 00521 s(r)++; 00522 } 00523 } 00524 00525 exit_point: 00526 LLR_aposteriori = bnum - bdenom; 00527 } 00528 00529 00530 void Modulator_NCD::update_norm(double &norm, int k, int sold, int snew, 00531 const cvec &ytH, const cmat &HtH, 00532 const ivec &s) 00533 { 00534 int m = length(s); 00535 std::complex<double> cdiff = symbols(k)[snew] - symbols(k)[sold]; 00536 00537 norm += sqr(cdiff) * (HtH(k, k).real()); 00538 cdiff *= 2.0; 00539 norm -= (cdiff.real() * ytH[k].real() - cdiff.imag() * ytH[k].imag()); 00540 for (int i = 0; i < m; i++) { 00541 norm += (cdiff * HtH(i, k) * conj(symbols(k)[s[i]])).real(); 00542 } 00543 } 00544 00545 00546 std::ostream &operator<<(std::ostream &os, const Modulator_NCD &mod) 00547 { 00548 os << "--- COMPLEX MIMO (NCD) CHANNEL --------" << std::endl; 00549 os << "Dimension (nt): " << mod.nt << std::endl; 00550 os << "Bits per dimension (k): " << mod.k << std::endl; 00551 os << "Symbols per dimension (M):" << mod.M << std::endl; 00552 for (int i = 0; i < mod.nt; i++) { 00553 os << "Bitmap for dimension " << i << ": " 00554 << mod.bitmap(i) << std::endl; 00555 os << "Symbol coordinates for dimension " << i << ": " 00556 << mod.symbols(i).left(mod.M(i)) << std::endl; 00557 } 00558 os << mod.get_llrcalc() << std::endl; 00559 return os; 00560 } 00561 00562 // ---------------------------------------------------------------------- 00563 // ND_UPAM 00564 // ---------------------------------------------------------------------- 00565 00566 ND_UPAM::ND_UPAM(int nt, int Mary) 00567 { 00568 set_M(nt, Mary); 00569 } 00570 00571 void ND_UPAM::set_M(int nt_in, int Mary) 00572 { 00573 nt = nt_in; 00574 ivec Mary_temp(nt); 00575 Mary_temp = Mary; 00576 set_M(nt, Mary_temp); 00577 } 00578 00579 void ND_UPAM::set_M(int nt_in, ivec Mary) 00580 { 00581 nt = nt_in; 00582 it_assert(length(Mary) == nt, "ND_UPAM::set_M(): Mary has wrong length"); 00583 k.set_size(nt); 00584 M = Mary; 00585 bitmap.set_size(nt); 00586 symbols.set_size(nt); 00587 bits2symbols.set_size(nt); 00588 spacing.set_size(nt); 00589 00590 for (int i = 0; i < nt; i++) { 00591 k(i) = round_i(::log2(static_cast<double>(M(i)))); 00592 it_assert((k(i) > 0) && ((1 << k(i)) == M(i)), 00593 "ND_UPAM::set_M(): M is not a power of 2."); 00594 00595 symbols(i).set_size(M(i) + 1); 00596 bits2symbols(i).set_size(M(i)); 00597 bitmap(i) = graycode(k(i)); 00598 double average_energy = (M(i) * M(i) - 1) / 3.0; 00599 double scaling_factor = std::sqrt(average_energy); 00600 00601 for (int j = 0; j < M(i); ++j) { 00602 symbols(i)(j) = ((M(i) - 1) - j * 2) / scaling_factor; 00603 bits2symbols(i)(bin2dec(bitmap(i).get_row(j))) = j; 00604 } 00605 00606 // the "symbols" vector must end with a zero; only for a trick 00607 // exploited in update_norm() 00608 symbols(i)(M(i)) = 0.0; 00609 00610 spacing(i) = 2.0 / scaling_factor; 00611 } 00612 } 00613 00614 int ND_UPAM::sphere_search_SE(const vec &y_in, const mat &H, 00615 const imat &zrange, double r, ivec &zhat) 00616 { 00617 // The implementation of this function basically follows the 00618 // Schnorr-Eucner algorithm described in Agrell et al. (IEEE 00619 // Trans. IT, 2002), but taking into account constellation 00620 // boundaries, see the "accelerated sphere decoder" in Boutros et 00621 // al. (IEEE Globecom, 2003). No lattice reduction is performed. 00622 // Potentially the function can be speeded up by performing 00623 // lattice reduction, but it seems difficult to keep track of 00624 // constellation boundaries. 00625 00626 mat R = chol(H.transpose() * H); 00627 mat Ri = inv(R); 00628 mat Q = H * Ri; 00629 vec y = Q.transpose() * y_in; 00630 mat Vi = Ri.transpose(); 00631 00632 int n = H.cols(); 00633 vec dist(n); 00634 dist(n - 1) = 0; 00635 double bestdist = r * r; 00636 int status = -1; // search failed 00637 00638 mat E = zeros(n, n); 00639 for (int i = 0; i < n; i++) { // E(k,:) = y*Vi; 00640 for (int j = 0; j < n; j++) { 00641 E(i*n + n - 1) += y(j) * Vi(j + n * i); 00642 } 00643 } 00644 00645 ivec z(n); 00646 zhat.set_size(n); 00647 z(n - 1) = floor_i(0.5 + E(n * n - 1)); 00648 z(n - 1) = std::max(z(n - 1), zrange(n - 1, 0)); 00649 z(n - 1) = std::min(z(n - 1), zrange(n - 1, 1)); 00650 double p = (E(n * n - 1) - z(n - 1)) / Vi(n * n - 1); 00651 ivec step(n); 00652 step(n - 1) = sign_nozero_i(p); 00653 00654 // Run search loop 00655 int k = n - 1; // k uses natural indexing, goes from 0 to n-1 00656 00657 while (true) { 00658 double newdist = dist(k) + p * p; 00659 00660 if ((newdist < bestdist) && (k != 0)) { 00661 for (int i = 0; i < k; i++) { 00662 E(k - 1 + i*n) = E(k + i * n) - p * Vi(k + i * n); 00663 } 00664 00665 k--; 00666 dist(k) = newdist; 00667 z(k) = floor_i(0.5 + E(k + k * n)); 00668 z(k) = std::max(z(k), zrange(k, 0)); 00669 z(k) = std::min(z(k), zrange(k, 1)); 00670 p = (E(k + k * n) - z(k)) / Vi(k + k * n); 00671 00672 step(k) = sign_nozero_i(p); 00673 } 00674 else { 00675 while (true) { 00676 if (newdist < bestdist) { 00677 zhat = z; 00678 bestdist = newdist; 00679 status = 0; 00680 } 00681 else if (k == n - 1) { 00682 goto exit_point; 00683 } 00684 else { 00685 k++; 00686 } 00687 00688 z(k) += step(k); 00689 00690 if ((z(k) < zrange(k, 0)) || (z(k) > zrange(k, 1))) { 00691 step(k) = (-step(k) - sign_nozero_i(step(k))); 00692 z(k) += step(k); 00693 } 00694 00695 if ((z(k) >= zrange(k, 0)) && (z(k) <= zrange(k, 1))) { 00696 break; 00697 } 00698 } 00699 00700 p = (E(k + k * n) - z(k)) / Vi(k + k * n); 00701 step(k) = (-step(k) - sign_nozero_i(step(k))); 00702 } 00703 } 00704 00705 exit_point: 00706 return status; 00707 } 00708 00709 00710 int ND_UPAM::sphere_decoding(const vec &y, const mat &H, double rstart, 00711 double rmax, double stepup, 00712 QLLRvec &detected_bits) 00713 { 00714 it_assert(H.rows() == length(y), 00715 "ND_UPAM::sphere_decoding(): dimension mismatch"); 00716 it_assert(H.cols() == nt, 00717 "ND_UPAM::sphere_decoding(): dimension mismatch"); 00718 it_assert(rstart > 0, "ND_UPAM::sphere_decoding(): radius error"); 00719 it_assert(rmax > rstart, "ND_UPAM::sphere_decoding(): radius error"); 00720 00721 // This function can be improved, e.g., by using an ordered search. 00722 00723 vec ytemp = y; 00724 mat Htemp(H.rows(), H.cols()); 00725 for (int i = 0; i < H.cols(); i++) { 00726 Htemp.set_col(i, H.get_col(i)*spacing(i)); 00727 ytemp += Htemp.get_col(i) * 0.5 * (M(i) - 1.0); 00728 } 00729 00730 imat crange(nt, 2); 00731 for (int i = 0; i < nt; i++) { 00732 crange(i, 0) = 0; 00733 crange(i, 1) = M(i) - 1; 00734 } 00735 00736 int status = 0; 00737 double r = rstart; 00738 ivec s(sum(M)); 00739 while (r <= rmax) { 00740 status = sphere_search_SE(ytemp, Htemp, crange, r, s); 00741 00742 if (status == 0) { // search successful 00743 detected_bits.set_size(sum(k)); 00744 int b = 0; 00745 for (int j = 0; j < nt; j++) { 00746 for (int i = 0; i < k(j); i++) { 00747 if (bitmap(j)((M(j) - 1 - s[j]), i) == 0) { 00748 detected_bits(b) = 1000; 00749 } 00750 else { 00751 detected_bits(b) = -1000; 00752 } 00753 b++; 00754 } 00755 } 00756 00757 return status; 00758 } 00759 r = r * stepup; 00760 } 00761 00762 return status; 00763 } 00764 00765 // ---------------------------------------------------------------------- 00766 // ND_UQAM 00767 // ---------------------------------------------------------------------- 00768 00769 // The ND_UQAM (MIMO with uniform QAM) class could alternatively 00770 // have been implemented by using a ND_UPAM class of twice the 00771 // dimension, but this does not fit as elegantly into the class 00772 // structure 00773 00774 ND_UQAM::ND_UQAM(int nt, int Mary) 00775 { 00776 set_M(nt, Mary); 00777 } 00778 00779 void ND_UQAM::set_M(int nt_in, int Mary) 00780 { 00781 nt = nt_in; 00782 ivec Mary_temp(nt); 00783 Mary_temp = Mary; 00784 set_M(nt, Mary_temp); 00785 } 00786 00787 void ND_UQAM::set_M(int nt_in, ivec Mary) 00788 { 00789 nt = nt_in; 00790 it_assert(length(Mary) == nt, "ND_UQAM::set_M(): Mary has wrong length"); 00791 k.set_size(nt); 00792 M = Mary; 00793 L.set_size(nt); 00794 bitmap.set_size(nt); 00795 symbols.set_size(nt); 00796 bits2symbols.set_size(nt); 00797 00798 for (int i = 0; i < nt; ++i) { 00799 k(i) = round_i(::log2(static_cast<double>(M(i)))); 00800 it_assert((k(i) > 0) && ((1 << k(i)) == M(i)), 00801 "ND_UQAM::set_M(): M is not a power of 2"); 00802 00803 L(i) = round_i(std::sqrt(static_cast<double>(M(i)))); 00804 it_assert(L(i)*L(i) == M(i), "ND_UQAM: constellation M must be square"); 00805 00806 symbols(i).set_size(M(i) + 1); 00807 bitmap(i).set_size(M(i), k(i)); 00808 bits2symbols(i).set_size(M(i)); 00809 double average_energy = (M(i) - 1) * 2.0 / 3.0; 00810 double scaling_factor = std::sqrt(average_energy); 00811 bmat gray_code = graycode(levels2bits(L(i))); 00812 00813 for (int j1 = 0; j1 < L(i); ++j1) { 00814 for (int j2 = 0; j2 < L(i); ++j2) { 00815 symbols(i)(j1*L(i) + j2) = 00816 std::complex<double>(((L(i) - 1) - j2 * 2.0) / scaling_factor, 00817 ((L(i) - 1) - j1 * 2.0) / scaling_factor); 00818 bitmap(i).set_row(j1*L(i) + j2, concat(gray_code.get_row(j1), 00819 gray_code.get_row(j2))); 00820 bits2symbols(i)(bin2dec(bitmap(i).get_row(j1*L(i) + j2))) 00821 = j1 * L(i) + j2; 00822 } 00823 } 00824 00825 // must end with a zero; only for a trick exploited in 00826 // update_norm() 00827 symbols(i)(M(i)) = 0.0; 00828 } 00829 } 00830 00831 // ---------------------------------------------------------------------- 00832 // ND_UPSK 00833 // ---------------------------------------------------------------------- 00834 00835 ND_UPSK::ND_UPSK(int nt, int Mary) 00836 { 00837 set_M(nt, Mary); 00838 } 00839 00840 void ND_UPSK::set_M(int nt_in, int Mary) 00841 { 00842 nt = nt_in; 00843 ivec Mary_temp(nt); 00844 Mary_temp = Mary; 00845 set_M(nt, Mary_temp); 00846 } 00847 00848 void ND_UPSK::set_M(int nt_in, ivec Mary) 00849 { 00850 nt = nt_in; 00851 it_assert(length(Mary) == nt, "ND_UPSK::set_M() Mary has wrong length"); 00852 k.set_size(nt); 00853 M = Mary; 00854 bitmap.set_size(nt); 00855 symbols.set_size(nt); 00856 bits2symbols.set_size(nt); 00857 00858 for (int i = 0; i < nt; ++i) { 00859 k(i) = round_i(::log2(static_cast<double>(M(i)))); 00860 it_assert((k(i) > 0) && ((1 << k(i)) == M(i)), 00861 "ND_UPSK::set_M(): M is not a power of 2"); 00862 00863 symbols(i).set_size(M(i) + 1); 00864 bits2symbols(i).set_size(M(i)); 00865 bitmap(i) = graycode(k(i)); 00866 00867 double delta = 2.0 * pi / M(i); 00868 double epsilon = delta / 10000.0; 00869 00870 for (int j = 0; j < M(i); ++j) { 00871 std::complex<double> symb 00872 = std::complex<double>(std::polar(1.0, delta * j)); 00873 00874 if (std::abs(std::real(symb)) < epsilon) { 00875 symbols(i)(j) = std::complex<double>(0.0, std::imag(symb)); 00876 } 00877 else if (std::abs(std::imag(symb)) < epsilon) { 00878 symbols(i)(j) = std::complex<double>(std::real(symb), 0.0); 00879 } 00880 else { 00881 symbols(i)(j) = symb; 00882 } 00883 00884 bits2symbols(i)(bin2dec(bitmap(i).get_row(j))) = j; 00885 } 00886 00887 // must end with a zero; only for a trick exploited in 00888 // update_norm() 00889 symbols(i)(M(i)) = 0.0; 00890 } 00891 } 00892 00893 } // namespace itpp
Generated on Tue Feb 2 09:33:31 2010 for IT++ by Doxygen 1.6.2