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