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   // 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
SourceForge Logo

Generated on Sat Apr 19 10:42:05 2008 for IT++ by Doxygen 1.5.5