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
Generated on Wed Apr 18 11:19:59 2007 for IT++ by Doxygen 1.5.2