IT++ Logo Newcom Logo

modulator_nd.cpp

Go to the documentation of this file.
00001 
00033 #include <itpp/comm/modulator_nd.h>
00034 #include <itpp/base/converters.h>
00035 #include <itpp/comm/commfunc.h>
00036 #include <itpp/base/elmatfunc.h>
00037 #include <itpp/base/cholesky.h>
00038 #include <itpp/base/inv.h>
00039 
00040 namespace itpp {
00041 
00042   // ----------------------- MIMO GENERAL -----------------------------
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     } else {
00052       result(0) = -llrcalc.jaclog(0,l);
00053       result(1) = result(0) + l;
00054     }
00055     return result;
00056   }
00057   
00058   Vec<QLLRvec> Modulator_ND::probabilities(QLLRvec &l)
00059   {
00060     Vec<QLLRvec> result(length(l));   
00061     for (int i=0; i<length(l); i++) {
00062       result(i) = probabilities(l(i));
00063     }
00064     return result;
00065   }
00066 
00067   void Modulator_ND::update_LLR(Vec<QLLRvec> &logP_apriori, QLLRvec &p1, QLLRvec &p0, int s, QLLR scaled_norm, int j)
00068   {
00069     QLLR log_apriori_prob_const_point = 0;
00070     int b=0;
00071     for (int i=0; i<k(j); i++) {
00072       log_apriori_prob_const_point += ((bitmap(j)(s,i)==0) ? logP_apriori(b)(1) : logP_apriori(b)(0));
00073       b++;
00074     }
00075     
00076     b=0;
00077     for (int i=0; i<k(j); i++) {
00078       if (bitmap(j)(s,i)==0) {
00079         p1(b) =  llrcalc.jaclog(p1(b), scaled_norm + log_apriori_prob_const_point );
00080       }  else {
00081         p0(b) = llrcalc.jaclog(p0(b),  scaled_norm + log_apriori_prob_const_point );
00082       }
00083       b++;
00084     }
00085   }
00086 
00087   void Modulator_ND::update_LLR(Vec<QLLRvec> &logP_apriori, QLLRvec &p1, QLLRvec &p0, ivec &s, QLLR scaled_norm)
00088   {
00089     QLLR log_apriori_prob_const_point = 0;
00090     int b=0;
00091     for (int j=0; j<nt; j++) {
00092       for (int i=0; i<k(j); i++) {
00093         log_apriori_prob_const_point += ((bitmap(j)(s[j],i)==0) ? logP_apriori(b)(1) : logP_apriori(b)(0));
00094         b++;
00095       }
00096     }
00097     
00098     b=0;
00099     for (int j=0; j<nt; j++) {
00100       for (int i=0; i<k(j); i++) {
00101         if (bitmap(j)(s[j],i)==0) {
00102           p1(b) =  llrcalc.jaclog(p1(b), scaled_norm + log_apriori_prob_const_point );
00103         }  else {
00104           p0(b) = llrcalc.jaclog(p0(b),  scaled_norm + log_apriori_prob_const_point );
00105         }
00106         b++;
00107       }
00108     }       
00109   }
00110 
00111   void Modulator_NRD::update_norm(double &norm, int k, int sold, int snew, vec &ytH, mat &HtH, ivec &s)
00112   {
00113    
00114     int m = length(s);
00115     double cdiff = symbols(k)[snew]-symbols(k)[sold];
00116     
00117     norm += sqr(cdiff)*HtH(k,k);
00118     cdiff = 2.0*cdiff;
00119     norm -= cdiff*ytH[k];
00120      for (int i=0; i<m; i++) {
00121       norm += cdiff*HtH(i,k)*symbols(k)[s[i]];
00122      }
00123   }
00124 
00125   void Modulator_NCD::update_norm(double &norm, int k, int sold, int snew, cvec &ytH, cmat &HtH, ivec &s)
00126   {
00127     int m = length(s);
00128     std::complex<double> cdiff = symbols(k)[snew]-symbols(k)[sold];
00129     
00130     norm += sqr(cdiff)*(HtH(k,k).real());
00131     cdiff = 2.0*cdiff;
00132     norm -= (cdiff.real()*ytH[k].real() - cdiff.imag()*ytH[k].imag());
00133     for (int i=0; i<m; i++) {
00134       norm += (cdiff*HtH(i,k)*conj(symbols(k)[s[i]])).real();
00135     }
00136   }
00137 
00138   void Modulator_NRD::map_demod(QLLRvec &LLR_apriori,  QLLRvec &LLR_aposteriori,  
00139                            double sigma2,  vec &h, vec &y)
00140   {
00141     it_assert(length(LLR_apriori)==sum(k),"Modulator_NRD::map_demod()");
00142     it_assert(length(LLR_apriori)==length(LLR_aposteriori),"Modulator_NRD::map_demod()");
00143     it_assert(length(h)==length(y),"Modulator_NRD::map_demod()");
00144     it_assert(length(h)==nt,"Modulator_NRD:map_demod()");
00145 
00146     int b=0;
00147     double oo_2s2 = 1.0/(2.0*sigma2);
00148     double norm2;
00149     QLLRvec temp, bnum, bdenom;
00150     for (int i=0; i<nt; i++) {
00151       temp=LLR_apriori(b,b+k(i)-1);
00152       bnum = (-QLLR_MAX)*ones_i(k(i));
00153       bdenom = (-QLLR_MAX)*ones_i(k(i));    
00154       Vec<QLLRvec> logP_apriori = probabilities(temp);
00155       for (int j=0; j<M(i); j++) {
00156         norm2 = sqr(y(i)-h(i)*symbols(i)(j));
00157         QLLR scaled_norm = llrcalc.to_qllr(-norm2*oo_2s2);
00158         update_LLR(logP_apriori, bnum, bdenom, j, scaled_norm,i);
00159       }
00160       LLR_aposteriori.set_subvector(b,bnum-bdenom);
00161       b+=k(i);
00162     }
00163   };
00164 
00165   void Modulator_NCD::map_demod(QLLRvec &LLR_apriori,  QLLRvec &LLR_aposteriori,  
00166                                 double sigma2,  cvec &h, cvec &y)
00167   {
00168     it_assert(length(LLR_apriori)==sum(k),"Modulator_NCD::map_demod()");
00169     it_assert(length(LLR_apriori)==length(LLR_aposteriori),"Modulator_NCD::map_demod()");
00170     it_assert(length(h)==length(y),"Modulator_NCD::map_demod()");
00171     it_assert(length(h)==nt,"Modulator_NCD:map_demod()");
00172 
00173     int b=0;
00174     double oo_s2 = 1.0/sigma2;
00175     double norm2;
00176     QLLRvec temp, bnum, bdenom;
00177     for (int i=0; i<nt; i++) {
00178       temp=LLR_apriori(b,b+k(i)-1);
00179       bnum = (-QLLR_MAX)*ones_i(k(i));
00180       bdenom = (-QLLR_MAX)*ones_i(k(i));    
00181       Vec<QLLRvec> logP_apriori = probabilities(temp);
00182       for (int j=0; j<M(i); j++) {
00183         norm2 = sqr(y(i)-h(i)*symbols(i)(j)); 
00184         QLLR scaled_norm = llrcalc.to_qllr(-norm2*oo_s2);
00185         update_LLR(logP_apriori, bnum, bdenom, j, scaled_norm,i);
00186       }
00187       LLR_aposteriori.set_subvector(b,bnum-bdenom);
00188       b+=k(i);
00189     }
00190   };
00191 
00192   void Modulator_NRD::map_demod(QLLRvec &LLR_apriori,  QLLRvec &LLR_aposteriori,  
00193                            double sigma2,  mat &H, vec &y)
00194   {
00195     int nr = H.rows();
00196     int np=sum(k); // number of bits in total
00197     it_assert(length(LLR_apriori)==np,"Modulator_NRD::map_demod()");
00198     it_assert(length(LLR_apriori)==length(LLR_aposteriori),"Modulator_NRD::map_demod()");
00199     it_assert(H.rows()==length(y),"Modulator_NRD::map_demod()");
00200     it_assert(H.cols()==nt,"Modulator_NRD:map_demod()");
00201     
00202     int mode=0;
00203      for (int i=0; i<length(M); i++) {
00204        if (nt*M(i)>4) { mode = 1; }    // differential update only pays off for larger dimensions
00205      }
00206     
00207     Vec<QLLRvec> logP_apriori = probabilities(LLR_apriori);
00208     
00209     mat HtH = H.transpose()*H; 
00210     vec ytH = H.transpose()*y; 
00211     
00212     QLLRvec bnum = (-QLLR_MAX)*ones_i(np);
00213     QLLRvec bdenom = (-QLLR_MAX)*ones_i(np);    
00214     ivec s(nt);
00215     s.clear();
00216     double norm = 0.0;
00217     double oo_2s2 = 1.0/(2.0*sigma2);
00218     
00219     // Go over all constellation points  (r=dimension, s=vector of symbols)
00220     int r = nt-1;
00221     while (1==1) {
00222       
00223       if (mode==1) {
00224         update_norm(norm, r, s[r], 0, ytH, HtH, s);
00225       }
00226       s[r] = 0;      
00227       while (1==1) {
00228         if (s[r] > M(r)-1)  {
00229           if (r==nt-1) {
00230             goto exit_point;
00231           }
00232           r++;
00233         }   else {
00234           if (r==0) {
00235             if (mode==0) {
00236               norm = 0.0;
00237               for (int p=0; p<nr; p++) {
00238                 double d = y[p];
00239                 for (int i=0; i<nt; i++) {
00240                   d -= H(p,i)*symbols(i)[s[i]];
00241                 }
00242                 norm += sqr(d);
00243               }
00244             }
00245             QLLR scaled_norm = llrcalc.to_qllr(-norm*oo_2s2);
00246             update_LLR(logP_apriori, bnum, bdenom, s, scaled_norm);
00247           } else {
00248             r--;
00249             break;
00250           }
00251         }
00252         if (mode==1) {
00253           update_norm(norm, r, s[r], s[r]+1, ytH, HtH, s);
00254         }
00255         s[r]++;
00256       }
00257     }
00258     
00259   exit_point:
00260     LLR_aposteriori = bnum - bdenom;
00261 
00262   };
00263 
00264   void Modulator_NCD::map_demod(QLLRvec &LLR_apriori,  QLLRvec &LLR_aposteriori,  double sigma2,  
00265                                 cmat &H, cvec &y)
00266   {
00267     int nr = H.rows();
00268     int np=sum(k); // number of bits in total
00269     it_assert(length(LLR_apriori)==np,"Modulator_NCD::map_demod()");
00270     it_assert(length(LLR_apriori)==length(LLR_aposteriori),"Modulator_NCD::map_demod()");
00271     it_assert(H.rows()==length(y),"Modulator_NCD::map_demod()");
00272     it_assert(H.cols()==nt,"Modulator_NCD:map_demod()");
00273     
00274     int mode=0;  
00275     for (int i=0; i<length(M); i++) {
00276       if (nt*M(i)>4) { mode = 1; }    // differential update only pays off for larger dimensions
00277     }
00278     
00279     Vec<QLLRvec> logP_apriori = probabilities(LLR_apriori);
00280     
00281     cmat HtH = H.hermitian_transpose()*H; 
00282     cvec ytH = conj(H.hermitian_transpose()*y); 
00283     
00284     QLLRvec bnum = (-QLLR_MAX)*ones_i(np);
00285     QLLRvec bdenom = (-QLLR_MAX)*ones_i(np);   
00286     ivec s(nt);
00287     s.clear();
00288     double norm = 0.0;
00289     double oo_s2 = 1.0/sigma2;
00290     std::complex<double> d;
00291 
00292     // Go over all constellation points  (r=dimension, s=vector of symbols)
00293     int r = nt-1;
00294     while (1==1) {
00295       
00296       if (mode==1) {
00297         update_norm(norm, r, s[r], 0, ytH, HtH, s);
00298       }      
00299       s[r] = 0;
00300       while (1==1) {
00301         if (s[r] > M(r)-1)  {
00302           if (r==nt-1) {
00303             goto exit_point;
00304           }
00305           r++;
00306         }   else {
00307           if (r==0) {
00308             if (mode==0) {
00309               norm = 0.0;
00310               for (int p=0; p<nr; p++) {
00311                 d = y[p];
00312                 for (int i=0; i<nt; i++) {
00313                   d -= H(p,i)*symbols(i)[s[i]];
00314                 }
00315                 norm += sqr(d);
00316               }
00317             }
00318             QLLR scaled_norm = llrcalc.to_qllr(-norm*oo_s2);
00319             update_LLR(logP_apriori, bnum, bdenom, s, scaled_norm);
00320           } else {
00321             r--;
00322             break;
00323           }
00324         }
00325         if (mode==1) {
00326           update_norm(norm, r, s[r], s[r]+1, ytH, HtH, s);
00327         }
00328         s[r]++;
00329       }
00330     }
00331     
00332   exit_point:
00333     LLR_aposteriori = bnum - bdenom;
00334   };
00335 
00336 
00337 
00338   vec Modulator_NRD::modulate_bits(const bvec &bits) const
00339   {
00340     vec result(nt);
00341 
00342     it_assert(length(bits)==sum(k),"Modulator_NRD::modulate_bits(): The number of input bits does not match.");
00343 
00344     int b=0;
00345     for (int i=0; i<nt; i++) {
00346       int symb = bin2dec(bits.mid(b,k(i)));
00347       result(i) = symbols(i)(bits2symbols(i)(symb));
00348       b+=k(i);
00349     }
00350 
00351     return result;
00352   };
00353 
00354 
00355   cvec Modulator_NCD::modulate_bits(const bvec &bits) const
00356   {
00357     cvec result(nt);
00358 
00359     it_assert(length(bits)==sum(k),"Modulator_NCD::modulate_bits(): The number of input bits does not match.");
00360 
00361     int b=0;
00362     for (int i=0; i<nt; i++) {
00363       int symb = bin2dec(bits.mid(b,k(i)));
00364       result(i) = symbols(i)(bits2symbols(i)(symb));
00365       b+=k(i);
00366     }
00367 
00368     return result;
00369   };
00370   
00371 
00372 
00373 
00374   std::ostream &operator<<(std::ostream &os, const Modulator_NRD &mod)
00375   {
00376     os << "--- REAL MIMO (NRD) CHANNEL ---------" << std::endl;
00377     os << "Dimension (nt):           " << mod.nt << std::endl;
00378     os << "Bits per dimension (k):   " << mod.k << std::endl;
00379     os << "Symbols per dimension (M):" << mod.M << std::endl;
00380     for (int i=0; i<mod.nt; i++) { 
00381       os << "Bitmap for dimension " << i << ": " << mod.bitmap(i) << std::endl;
00382       // skip printing the trailing zero
00383       os << "Symbol coordinates for dimension " << i << ": " << mod.symbols(i).left(mod.M(i)) << std::endl;
00384     }
00385     os << mod.get_llrcalc() << std::endl;
00386     return os;
00387   };
00388   
00389   std::ostream &operator<<(std::ostream &os, const Modulator_NCD &mod)
00390   {
00391     os << "--- COMPLEX MIMO (NCD) CHANNEL --------" << std::endl;
00392     os << "Dimension (nt):           " << mod.nt << std::endl;
00393     os << "Bits per dimension (k):   " << mod.k << std::endl;
00394     os << "Symbols per dimension (M):" << mod.M << std::endl;
00395     for (int i=0; i<mod.nt; i++) { 
00396       os << "Bitmap for dimension " << i << ": " << mod.bitmap(i) << std::endl;
00397       os << "Symbol coordinates for dimension " << i << ": " << mod.symbols(i).left(mod.M(i)) << std::endl;
00398     }
00399     os << mod.get_llrcalc() << std::endl;
00400     return os;
00401   };
00402 
00403   // ------------------------- MIMO with uniform PAM ----------------------------
00404 
00405   ND_UPAM::ND_UPAM(int nt_in, int Mary)
00406   {
00407     nt=nt_in;
00408     set_Gray_PAM(nt,Mary);
00409   };
00410   
00411   void ND_UPAM::set_Gray_PAM(int nt_in, int Mary)
00412   {
00413     nt=nt_in;
00414     ivec Mary_temp(nt);
00415     for (int i=0; i<nt; i++) {
00416       Mary_temp(i)=Mary;
00417     }
00418     set_Gray_PAM(nt,Mary_temp);
00419   };
00420 
00421   void ND_UPAM::set_Gray_PAM(int nt_in, ivec Mary)
00422   {
00423     nt=nt_in;
00424     it_assert(length(Mary)==nt,"ND_UPAM::set_Gray_PAM() Mary has wrong length");
00425     k.set_size(nt);
00426     M=Mary;
00427     bitmap.set_size(nt);
00428     symbols.set_size(nt);
00429     bits2symbols.set_size(nt);    
00430     spacing.set_size(nt);
00431 
00432     for (int i=0; i<nt; i++) {
00433       k(i) = round_i(log2(double(M(i))));
00434       it_assert( ((k(i)>0) && ((1<<k(i))==M(i))),"ND_UPAM::set_Gray_PAM(): M is not a power of 2.");
00435        
00436       symbols(i).set_size(M(i)+1);
00437       bits2symbols(i).set_size(M(i));
00438       bitmap(i) = graycode(k(i));
00439       double average_energy = double(M(i)*M(i)-1.0)/3.0;
00440       double scaling_factor = std::sqrt(average_energy);
00441       
00442       for (int j=0; j<M(i); j++) {
00443         symbols(i)(j) = double((M(i)-1)-j*2) / scaling_factor;
00444         bits2symbols(i)(bin2dec(bitmap(i).get_row(j))) = j;
00445       }      
00446       
00447       symbols(i)(M(i))=0.0;  // must end with a zero; only for a trick exploited in update_norm()
00448 
00449       spacing(i)=2.0/scaling_factor;
00450     }
00451   };
00452   
00453 
00454 #define max(a,b) (((a)>(b))?(a):(b))
00455 #define min(a,b) (((a)<(b))?(a):(b))
00456 #define sign(a)  (((a)>0)?(1):(-1))
00457 
00458   int ND_UPAM::sphere_search_SE(vec &y_in, mat &H, imat &zrange, double r, ivec &zhat)
00459   {
00460     // The implementation of this function basically follows the
00461     // Schnorr-Eucner algorithm described in Agrell et al. (IEEE
00462     // Trans.  IT, 2002), but taking into account constellation
00463     // boundaries, see the "accelerated sphere decoder" in Boutros et
00464     // al. (IEEE Globecom, 2003).  No lattice reduction is performed.
00465     // Potentially the function can be speeded up by performing
00466     // lattice reduction, but it seems difficult to keep track of
00467     // constellation boundaries.
00468         
00469     mat R=chol(H.transpose()*H);
00470     mat Ri=inv(R);
00471     mat Q=H*Ri;
00472     vec y=Q.transpose()*y_in;
00473     mat Vi=Ri.transpose();
00474 
00475     int n=H.cols();
00476     vec dist(n);
00477     dist[n-1] = 0;
00478     double bestdist = r*r;
00479     int status = -1; // search failed 
00480     
00481     mat E=zeros(n,n);
00482     for (int i=0; i<n; i++) {   // E(k,:) = y*Vi; 
00483       for (int j=0; j<n; j++) {
00484         E(i*n+n-1) += y(j)*Vi(j+n*i);
00485       }
00486     }
00487     
00488     ivec z(n);
00489     zhat.set_size(n);
00490     z(n-1) = (int) std::floor(0.5 + E(n*n-1));
00491     z(n-1) = max(z(n-1),zrange(n-1,0));
00492     z(n-1) = min(z(n-1),zrange(n-1,1));
00493     double p = (E(n*n-1)-z(n-1))/Vi(n*n-1);
00494     ivec step(n);
00495     step(n-1) = sign(p);
00496     
00497     // Run search loop 
00498     int k=n-1;  // k uses natural indexing, goes from 0 to n-1 
00499     
00500     while (1==1) {
00501       double newdist = dist(k) + p*p;
00502       
00503       if ((newdist<bestdist) && (k!=0)) {
00504         
00505         for (int i=0; i<k; i++) {
00506           E(k-1+i*n) = E(k+i*n) - p*Vi(k+i*n);
00507         }
00508         
00509         k--;
00510         dist(k) = newdist;
00511         z(k) = (int) std::floor(0.5+E(k+k*n));
00512         z(k) = max(z(k),zrange(k,0));
00513         z(k) = min(z(k),zrange(k,1));
00514         p = (E(k+k*n)-z(k))/Vi(k+k*n);
00515         
00516         step[k] = sign(p);
00517       } else {
00518         while (1==1) {
00519           if (newdist<bestdist) {
00520             for (int j=0; j<n; j++) { zhat(j) = z(j); }
00521             bestdist = newdist; 
00522             status = 0;
00523           }
00524           else if (k==n-1) {
00525             goto exit_point; 
00526           } else {
00527             k++;
00528           }
00529           
00530           z[k] += step(k);
00531           
00532           if ((z(k)<zrange(k,0)) || (z(k)>zrange(k,1))) {
00533             step(k) = (-step(k) - sign(step(k)));
00534             z(k) += step(k);
00535           }
00536           
00537           if ((z(k)>=zrange(k,0)) && (z(k)<=zrange(k,1))) {
00538             break;
00539           }
00540         }
00541         
00542         p = (E(k+k*n)-z(k))/Vi(k+k*n);
00543         step(k) = (-step(k) - sign(step(k)));
00544       }
00545     }
00546     
00547   exit_point:
00548     return status;
00549   }
00550   
00551   
00552   int ND_UPAM::sphere_decoding(vec &y, mat &H, double rstart, double rmax, 
00553                                double stepup, QLLRvec &detected_bits)
00554   {
00555     it_assert(H.rows()==length(y),"ND_UPAM::sphere_decoding(): dimension mismatch");
00556     it_assert(H.cols()==nt,"ND_UPAM::sphere_decoding(): dimension mismatch");
00557     it_assert(rstart>0,"ND_UPAM::sphere_decoding(): radius error");
00558     it_assert(rmax>rstart,"ND_UPAM::sphere_decoding(): radius error");
00559     
00560     // This function can be improved, e.g., by using an ordered search.
00561 
00562     vec ytemp=y;
00563     mat Htemp(H.rows(),H.cols());
00564     for (int i=0; i<H.cols(); i++) {
00565       Htemp.set_col(i,H.get_col(i)*spacing(i));
00566       ytemp+=Htemp.get_col(i)*0.5*(M(i)-1.0);
00567     }
00568 
00569     imat crange(nt,2);
00570     for (int i=0; i<nt; i++) {
00571       crange(i,0) = 0;
00572       crange(i,1) = M(i)-1;
00573     }
00574 
00575     int status;  
00576     double r = rstart;
00577     ivec s(sum(M));
00578     while (r<=rmax) {
00579       status = sphere_search_SE(ytemp,Htemp,crange,r,s);
00580       
00581       if (status==0) { // search successful
00582         detected_bits.set_size(sum(k));
00583         int b=0;
00584         for (int j=0; j<nt; j++) {
00585           for (int i=0; i<k(j); i++) {
00586             if (bitmap(j)((M(j)-1-s[j]),i)==0) {
00587               detected_bits(b) = 1000;
00588             }  else {
00589               detected_bits(b) = -1000;
00590             }
00591             b++;
00592           }
00593         }           
00594         
00595         return status;
00596       }
00597       r = r*stepup;
00598     }
00599 
00600     return status;
00601     
00602   }
00603   
00604 
00605   // -------------------- MIMO with uniform QAM ---------------------------- 
00606 
00607   // The ND_UQAM class could alternatively have been implemented by
00608   // using a ND_UPAM class of twice the dimension, but this does not
00609   // fit as elegantly into the class structure
00610 
00611   ND_UQAM::ND_UQAM(int nt_in, int Mary)
00612   {
00613     nt=nt_in;
00614     set_Gray_QAM(nt,Mary);
00615   };
00616 
00617   void ND_UQAM::set_Gray_QAM(int nt_in, int Mary)
00618   {
00619     nt=nt_in;
00620     ivec Mary_temp(nt);
00621     for (int i=0; i<nt; i++) {
00622       Mary_temp(i)=Mary;
00623     }
00624     set_Gray_QAM(nt,Mary_temp);
00625   };
00626 
00627   void ND_UQAM::set_Gray_QAM(int nt_in, ivec Mary)
00628   {
00629     nt=nt_in;
00630     it_assert(length(Mary)==nt,"ND_UQAM::set_Gray_QAM() Mary has wrong length");
00631     k.set_size(nt);
00632     M=Mary;
00633     L.set_size(nt);
00634     bitmap.set_size(nt);
00635     symbols.set_size(nt);
00636     bits2symbols.set_size(nt);    
00637 
00638     for (int i=0; i<nt; i++) {
00639       k(i) = round_i(log2(double(M(i))));      
00640       it_assert( ((k(i)>0) && ((1<<k(i))==M(i))),"ND_UQAM::set_Gray_QAM(): M is not a power of 2.");
00641 
00642       L(i) = round_i(std::sqrt((double)M(i)));
00643       it_assert(L(i)*L(i)== M(i),"ND_UQAM: constellation M must be square");
00644       
00645       symbols(i).set_size(M(i)+1);
00646       bitmap(i).set_size(M(i),k(i));
00647       bits2symbols(i).set_size(M(i));
00648       double average_energy = double(M(i)-1)*2.0/3.0;
00649       double scaling_factor = std::sqrt(average_energy);
00650       bmat gray_code = graycode(levels2bits(L(i)));
00651       
00652       for (int j1=0; j1<L(i); j1++) {
00653         for (int j2=0; j2<L(i); j2++) {
00654           symbols(i)(j1*L(i)+j2) = std::complex<double>( ((L(i)-1)-j2*2.0)/scaling_factor,
00655                                                          ((L(i)-1)-j1*2.0)/scaling_factor );
00656           bitmap(i).set_row(j1*L(i)+j2, concat(gray_code.get_row(j1), 
00657                                                gray_code.get_row(j2)));
00658           bits2symbols(i)( bin2dec(bitmap(i).get_row(j1*L(i)+j2)) ) = j1*L(i)+j2;
00659         }
00660       }
00661       
00662       symbols(i)(M(i))=0.0;  // must end with a zero; only for a trick exploited in update_norm()
00663     }    
00664       
00665   };
00666 
00667   
00668 
00669 
00670   // ----------- MIMO with uniform PSK ---------------
00671 
00672   ND_UPSK::ND_UPSK(int nt_in, int Mary)
00673   {
00674     nt=nt_in;
00675     set_Gray_PSK(nt,Mary);
00676   };
00677 
00678   void ND_UPSK::set_Gray_PSK(int nt_in, int Mary) 
00679   { 
00680     nt=nt_in;
00681     ivec Mary_temp(nt);
00682     for (int i=0; i<nt; i++) {
00683       Mary_temp(i)=Mary;
00684     }
00685     set_Gray_PSK(nt,Mary_temp);
00686   };
00687   
00688   void ND_UPSK::set_Gray_PSK(int nt_in, ivec Mary) 
00689   { 
00690     nt=nt_in;
00691     it_assert(length(Mary)==nt,"ND_UPSK::set_Gray_PSK() Mary has wrong length");
00692     k.set_size(nt);
00693     M=Mary;
00694     bitmap.set_size(nt);
00695     symbols.set_size(nt);
00696     bits2symbols.set_size(nt);    
00697    
00698     for (int i=0; i<nt; i++) {
00699       k(i) = round_i(log2(double(M(i))));      
00700       it_assert( ((k(i)>0) && ((1<<k(i))==M(i))),"ND_UPSK::set_Gray_PSK(): M is not a power of 2.");
00701        
00702       symbols(i).set_size(M(i)+1);
00703       bits2symbols(i).set_size(M(i));
00704       bitmap(i) = graycode(k(i));
00705 
00706       double delta = 2.0*pi/M(i);
00707       double epsilon = delta/10000.0;
00708       
00709       for (int j=0; j<M(i); j++) {
00710         std::complex<double> symb = std::complex<double>(std::polar(1.0,delta*j));
00711 
00712         if (std::abs(std::real(symb)) < epsilon) { 
00713           symbols(i)(j) = std::complex<double>(0.0,std::imag(symb)); 
00714         } else if (std::abs(std::imag(symb)) < epsilon) { 
00715           symbols(i)(j) = std::complex<double>(std::real(symb),0.0); }
00716         else { 
00717           symbols(i)(j) = symb; 
00718         }
00719 
00720         bits2symbols(i)(bin2dec(bitmap(i).get_row(j))) = j;
00721       }      
00722       
00723       symbols(i)(M(i))=0.0;  // must end with a zero; only for a trick exploited in update_norm()
00724     }
00725   };
00726   
00727 
00728 } // namespace itpp
SourceForge Logo

Generated on Wed Apr 18 11:45:35 2007 for IT++ by Doxygen 1.5.2