IT++ Logo

modulator_nd.cpp

Go to the documentation of this file.
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
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
SourceForge Logo

Generated on Tue Nov 23 08:47:56 2010 for IT++ by Doxygen 1.6.1