IT++ Logo Newcom Logo

modulator.cpp

Go to the documentation of this file.
00001 
00033 #include <itpp/comm/modulator.h>
00034 #include <itpp/base/converters.h>
00035 #include <itpp/base/stat.h>
00036 #include <itpp/comm/commfunc.h>
00037 
00038 
00039 namespace itpp {
00040 
00041 #ifndef M_SQRT1_2
00042 #define M_SQRT1_2 0.70710678118654752440
00043 #endif
00044 
00045   //------------- class: Modulator_1d ----------------
00046 
00047   Modulator_1d::Modulator_1d(const vec &symbols, const ivec &bitmap)
00048   {
00049     set(symbols, bitmap);
00050   }
00051 
00052   void Modulator_1d::set(const vec &in_symbols, const ivec &in_bitmap)
00053   {
00054     it_assert(in_symbols.size() == in_bitmap.size(), "Modulator_1d::set(): Number of symbols and bitmap does not match");
00055     M = in_bitmap.size();
00056     k = round_i(log2(double(M)));
00057     symbols = in_symbols;
00058     bitmap = in_bitmap; 
00059   }
00060 
00061   vec Modulator_1d::modulate(const ivec &symbolnumbers) const
00062   {
00063     vec temp(symbolnumbers.size());
00064     for (int i = 0; i < symbolnumbers.size(); i++)
00065       temp(i) = symbols(symbolnumbers(i));
00066     return temp;
00067   }
00068 
00069   vec Modulator_1d::modulate_bits(const bvec &bits) const
00070   {
00071     int no_symbols = bits.length() / k;
00072     int pos, symb;
00073     vec output = zeros(no_symbols);
00074 
00075     for (int i = 0; i < no_symbols; i++) {
00076       pos = 0;
00077       symb = bin2dec(bits.mid(i*k, k));
00078       while (bitmap(pos) != symb) { pos++; }
00079       output(i) = symbols(pos);
00080     }
00081     return output;
00082   }
00083 
00084   ivec Modulator_1d::demodulate(const vec &signal) const
00085   {
00086     double dist, mindist;
00087     int closest;
00088     ivec output(signal.size());
00089 
00090     for (int i = 0; i < signal.size(); i++) {
00091       mindist = std::fabs(symbols(0) - signal(i));
00092       closest = 0;
00093       for (int j = 1; j < M; j++) {
00094         dist = std::fabs(symbols(j) - signal(i));
00095         if (dist < mindist) {
00096           mindist = dist;
00097           closest = j;
00098         }
00099       }
00100       output(i) = closest;
00101     }
00102     return output;
00103   }
00104 
00105   bvec Modulator_1d::demodulate_bits(const vec &signal) const
00106   {
00107     double dist, mindist;
00108     int closest;
00109     bvec output(k*signal.size());
00110 
00111     for (int i = 0; i < signal.size(); i++) {
00112       mindist = std::fabs(symbols(0) - signal(i));
00113       closest = 0;
00114       for (int j = 1; j < M; j++) {
00115         dist = std::fabs(symbols(j) - signal(i));
00116         if (dist < mindist) { 
00117           mindist = dist;
00118           closest = j;
00119         }
00120       }
00121       output.replace_mid(i*k, dec2bin(k, bitmap(closest)));
00122     }
00123     return output;
00124   }
00125 
00126 
00127   //------------- class: Modulator_2d ----------------
00128   Modulator_2d::Modulator_2d(const cvec &symbols, const ivec &bitmap)
00129   {
00130     set(symbols, bitmap);
00131   }
00132 
00133   void Modulator_2d::set(const cvec &in_symbols, const ivec &in_bitmap)
00134   {
00135     it_assert(in_symbols.size() == in_bitmap.size(),
00136               "Modulator_2d::set() Number of symbols and bitmap does not match");
00137     it_assert(is_even(in_symbols.size()) && (in_symbols.size() > 0), 
00138               "Modulator_2d::set(): Number of symbols needs to be even and non-zero");
00139     it_assert((max(in_bitmap) == in_bitmap.size() - 1) && (min(in_bitmap) == 0),
00140               "Modulator_2d::set(): Improper bitmap vector");
00141     symbols = in_symbols;
00142     bitmap = in_bitmap; 
00143     M = bitmap.size();
00144     k = levels2bits(M);
00145     calculate_softbit_matricies(bitmap);
00146   }
00147 
00148   cvec Modulator_2d::modulate(const ivec &symbolnumbers) const
00149   {
00150     cvec temp(symbolnumbers.length());
00151     for(int i = 0; i < symbolnumbers.length(); i++)
00152       temp(i) = symbols(symbolnumbers(i));
00153     return temp;
00154   }
00155 
00156   void Modulator_2d::modulate_bits(const bvec &bits, cvec &output) const
00157   {
00158     int no_symbols = bits.length() / k;
00159     int pos, symb;
00160     output.set_size(no_symbols);
00161 
00162     for (int i = 0; i < no_symbols; i++) {
00163       pos = 0;
00164       symb = bin2dec(bits.mid(i*k, k));
00165       while (bitmap(pos) != symb) { pos++; }
00166       output(i) = symbols(pos);
00167     }
00168   }
00169 
00170   cvec Modulator_2d::modulate_bits(const bvec &bits) const
00171   {
00172     cvec output;
00173     modulate_bits(bits, output);
00174     return output;
00175   }
00176 
00177   ivec Modulator_2d::demodulate(const cvec &signal) const
00178   {
00179     double dist, mindist;
00180     int closest;
00181     ivec output(signal.size());
00182 
00183     for (int i = 0; i < signal.size(); i++) {
00184       mindist = std::abs(symbols(0) - signal(i));
00185       closest = 0;
00186       for (int j = 1; j < M; j++) {
00187         dist = std::abs(symbols(j) - signal(i));
00188         if (dist < mindist) {
00189           mindist = dist;
00190           closest = j;
00191         }
00192       }
00193       output(i) = closest;
00194     }
00195     return output;
00196   }
00197 
00198   void Modulator_2d::demodulate_bits(const cvec &signal, bvec &bits) const
00199   {
00200     double dist, mindist;
00201     int closest;
00202     bits.set_size(k*signal.size());
00203 
00204     for (int i = 0; i < signal.size(); i++) {
00205       mindist = std::abs(symbols(0) - signal(i));
00206       closest = 0;
00207       for (int j = 1; j < M; j++) {
00208         dist = std::abs(symbols(j) - signal(i));
00209         if (dist < mindist) {
00210           mindist = dist;
00211           closest = j;
00212         }
00213       }
00214       bits.replace_mid(i*k, dec2bin(k, bitmap(closest)));
00215     }
00216   }
00217 
00218   bvec Modulator_2d::demodulate_bits(const cvec &signal) const
00219   {
00220     bvec bits;
00221     demodulate_bits(signal, bits);
00222     return bits;
00223   }
00224 
00225   void Modulator_2d::demodulate_soft_bits(const cvec &rx_symbols, double N0,
00226                                           vec &soft_bits) const
00227   {
00228     // Definitions of local variables
00229     long no_symbols = rx_symbols.length();
00230     double P0, P1;
00231     vec soft_word(k);
00232 
00233     // Allocate storage space for the result vector:
00234     soft_bits.set_size(k*no_symbols, false);
00235 
00236     // For each symbol l do:
00237     for (int l = 0; l < no_symbols; l++) {
00238       // For each bit position i do:
00239       for (int i = 0; i < k; i++) {
00240         P0 = 0;
00241         P1 = 0;
00242         for (int j = 0; j < (M >> 1); j++) {
00243 //        s0 = symbols( S0(i,j) );
00244 //        s1 = symbols( S1(i,j) );
00245 //        p_z_s0 = (1.0/(pi*N0)) * std::exp(-sqr(std::abs(z - s0)) / N0);
00246 //        p_z_s1 = (1.0/(pi*N0)) * std::exp(-sqr(std::abs(z - s1)) / N0);
00247 //        P0(i) += p_z_s0;
00248 //        P1(i) += p_z_s1;
00249           // sqr(complex) = sqr(abs(complex))
00250           P0 += std::exp(-sqr(rx_symbols(l) - symbols(S0(i, j))) / N0);
00251           P1 += std::exp(-sqr(rx_symbols(l) - symbols(S1(i, j))) / N0);
00252         }
00253         // The soft bits for the l-th received symbol:
00254         soft_word(i) = trunc_log(P0) - trunc_log(P1);
00255       }
00256       // Put the results in the result vector:
00257       soft_bits.replace_mid(l*k, soft_word);
00258     }
00259 
00260   }
00261 
00262   void Modulator_2d::demodulate_soft_bits(const cvec &rx_symbols, 
00263                                           const cvec &chan, double N0, 
00264                                           vec &soft_bits) const
00265   {
00266     // Definitions of local variables
00267     long no_symbols = rx_symbols.length();
00268     double a2, P0, P1;
00269     vec soft_word(k);
00270 
00271     //Allocate storage space for the result vector:
00272     soft_bits.set_size(k*no_symbols,false);
00273 
00274     // For each symbol l do:
00275     for (int l = 0; l < no_symbols; l++) {
00276       a2 = sqr(chan(l));
00277       // For each bit position i do:
00278       for (int i = 0; i < k; i++) {
00279         P0 = 0;
00280         P1 = 0;
00281         for (int j = 0; j < (M >> 1); j++) {
00282 //        s0 = symbols( S0(i,j) );
00283 //        s1 = symbols( S1(i,j) );
00284 //        p_z_s0 = (1.0/(pi*N0*a2)) * std::exp(-sqr(std::abs(z-a2*s0)) / (N0*a2));
00285 //        p_z_s1 = (1.0/(pi*N0*a2)) * std::exp(-sqr(std::abs(z-a2*s1)) / (N0*a2));
00286 //        P0(i) += p_z_s0;
00287 //        P1(i) += p_z_s1;
00288           // sqr(complex) = sqr(abs(complex))
00289           P0 += std::exp(-sqr(rx_symbols(l) - a2*symbols(S0(i, j))) / (N0*a2));
00290           P1 += std::exp(-sqr(rx_symbols(l) - a2*symbols(S1(i, j))) / (N0*a2));
00291         }
00292         // The soft bits for the l-th received symbol:
00293         soft_word(i) = trunc_log(P0) - trunc_log(P1);
00294       }
00295       // Put the results in the result vector:
00296       soft_bits.replace_mid(l*k, soft_word);
00297     }
00298   }
00299 
00300   void Modulator_2d::demodulate_soft_bits_approx(const cvec &rx_symbols, 
00301                                                  double N0, 
00302                                                  vec &soft_bits) const
00303   {
00304     // Definitions of local variables
00305     long no_symbols = rx_symbols.length();
00306     double d_0, d_1;
00307     vec d(M);
00308 
00309     // Allocate storage space for the result vector:
00310     soft_bits.set_size(k*no_symbols, false);
00311 
00312     // For each symbol l do:
00313     for (int l = 0; l < no_symbols; l++) {
00314       // Calculate all distances:
00315       for (int i = 0; i < M; i++) { 
00316         d(i) = std::abs(rx_symbols(l) - symbols(i));
00317       }
00318       //For each of the k bits do:
00319       for (int i = 0; i < k; i++) {
00320         // Find the closest 0-point and the closest 1-point:
00321         d_0 = d(S0(i, 0));
00322         d_1 = d(S1(i, 0));
00323         for (int j = 1; j < (M >> 1); j++) {
00324           if (d(S0(i, j)) < d_0) { d_0 = d(S0(i, j)); }
00325           if (d(S1(i, j)) < d_1) { d_1 = d(S1(i, j)); }
00326         }
00327         //calculate the approximative metric:
00328         soft_bits(l*k+i) = (sqr(d_1) - sqr(d_0)) / N0;
00329       }
00330     }
00331   }
00332 
00333   void Modulator_2d::demodulate_soft_bits_approx(const cvec &rx_symbols,
00334                                                  const cvec &chan, double N0,
00335                                                  vec &soft_bits) const
00336   {
00337     // Definitions of local variables
00338     long no_symbols = rx_symbols.length();
00339     double d_0, d_1, Kf, c2;
00340     vec d(M);
00341 
00342     // Allocate storage space for the result vector:
00343     soft_bits.set_size(k*no_symbols, false);
00344 
00345     // For each symbol l do:
00346     for (int l = 0; l < no_symbols; l++) {
00347       c2 = sqr(chan(l)); 
00348       Kf = 1.0 / ( c2 * N0 );   
00349       // Calculate all distances:
00350       for (int i = 0; i < M; i++) {
00351         d(i) = std::abs(rx_symbols(l) - c2 * symbols(i));
00352       }
00353       // For each of the k bits do:
00354       for (int i = 0; i < k; i++) {
00355         //Find the closest 0-point and the closest 1-point:
00356         d_0 = d(S0(i, 0));
00357         d_1 = d(S1(i, 0));
00358         for (int j = 1; j < (M >> 1); j++) {
00359           if (d(S0(i, j)) < d_0) { d_0 = d(S0(i, j)); }
00360           if (d(S1(i, j)) < d_1) { d_1 = d(S1(i, j)); }
00361         }
00362         //calculate the approximative metric:
00363         soft_bits(l*k+i) = Kf * (sqr(d_1) - sqr(d_0));
00364       }
00365     }
00366   }
00367 
00368   void Modulator_2d::calculate_softbit_matricies(const ivec& inbitmap)
00369   {
00370     // Definitions of local variables
00371     int count0, count1;
00372     bvec bits(k);
00373 
00374     // Allocate storage space for the result matricies:
00375     S0.set_size(k, M/2, false);
00376     S1.set_size(k, M/2, false);
00377 
00378     for (int kk = 0; kk < k; kk++) {
00379       count0 = 0; 
00380       count1 = 0;
00381       for (int m = 0; m < M; m++) {
00382         bits = dec2bin(k, inbitmap(m));
00383         if (bits(kk) == bin(0)) {
00384           S0(kk, count0) = m;
00385           count0++; 
00386         } else {
00387           S1(kk, count1) = m;
00388           count1++;
00389         }
00390       }
00391     }
00392   }
00393 
00394 
00395   //------------- class: BPSK ----------------
00396 
00397   void BPSK::modulate_bits(const bvec &bits, vec &out) const
00398   {
00399     out.set_size(bits.size(), false);
00400     for (int i = 0; i < bits.size(); i++) { 
00401       out(i) = (bits(i) == 0 ? 1.0 : -1.0);
00402     }
00403   }
00404 
00405   // output is complex but symbols in real part only
00406   void BPSK::modulate_bits(const bvec &bits, cvec &out) const
00407   {
00408     out.set_size(bits.size(), false);
00409     for (int i = 0; i < bits.size(); i++) { 
00410       out(i) = (bits(i) == 0 ? 1.0 : -1.0);
00411     }
00412   }
00413 
00414   cvec BPSK::modulate_bits(const bvec &bits) const
00415   {
00416     cvec temp(bits.size());
00417     modulate_bits(bits, temp);
00418     return temp;
00419   }
00420 
00421   void BPSK::demodulate_bits(const vec &signal, bvec &out) const
00422   {
00423     out.set_size(signal.size(), false);
00424     for (int i = 0; i < signal.length(); i++) { 
00425       out(i) = (signal(i) > 0) ? bin(0) : bin(1); 
00426     }
00427   }
00428 
00429   bvec BPSK::demodulate_bits(const vec &signal) const
00430   {
00431     bvec temp(signal.size());
00432     demodulate_bits(signal, temp);
00433     return temp;
00434   }
00435 
00436   // Symbols are in real part. Channel estimation already applied to the
00437   // signal
00438   void BPSK::demodulate_bits(const cvec &signal, bvec &out) const
00439   {
00440     out.set_size(signal.size(), false);
00441     for (int i = 0; i < signal.length(); i++) { 
00442       out(i) = (std::real(signal(i)) > 0) ? bin(0) : bin(1); 
00443     }
00444   }
00445 
00446   bvec BPSK::demodulate_bits(const cvec &signal) const
00447   {
00448     bvec temp(signal.size());
00449     demodulate_bits(signal, temp);
00450     return temp;
00451   }
00452 
00453   // Outputs log-likelihood ratio of log(Pr(b=0|rx_symbols)/Pr(b=1|rx_symbols))
00454   void BPSK::demodulate_soft_bits(const vec &rx_symbols, double N0,
00455                                   vec &soft_bits) const
00456   {
00457     double factor = 4 / N0;
00458     soft_bits.set_size(rx_symbols.size(), false);
00459 
00460     for (int i = 0; i < rx_symbols.size(); i++) {
00461       soft_bits(i) = factor * rx_symbols(i);
00462     }
00463   }
00464 
00465   // Outputs log-likelihood ratio of log(Pr(b=0|rx_symbols)/Pr(b=1|rx_symbols))
00466   void BPSK::demodulate_soft_bits(const cvec &rx_symbols, double N0,
00467                                   vec &soft_bits) const
00468   {
00469     double factor = 4 / N0;
00470     soft_bits.set_size(rx_symbols.size(), false);
00471 
00472     for (int i = 0; i < rx_symbols.size(); i++) {
00473       soft_bits(i) = factor * std::real(rx_symbols(i));
00474     }
00475   }
00476 
00477   // Outputs log-likelihood ratio for fading channels
00478   void BPSK::demodulate_soft_bits(const cvec &rx_symbols, const cvec &channel,
00479                                   double N0, vec &soft_bits) const
00480   {
00481     double factor = 4 / N0;
00482     soft_bits.set_size(rx_symbols.size(), false);
00483 
00484     for (int i = 0; i < rx_symbols.size(); i++) {
00485       soft_bits(i) = factor * std::real(rx_symbols(i) * std::conj(channel(i)));
00486     }
00487   }
00488 
00489   void BPSK::demodulate_soft_bits_approx(const cvec &rx_symbols, double N0,
00490                                          vec &soft_bits) const
00491   {
00492     demodulate_soft_bits(rx_symbols, N0, soft_bits);
00493   }
00494 
00495   void BPSK::demodulate_soft_bits_approx(const cvec &rx_symbols, 
00496                                          const cvec &channel, double N0,
00497                                          vec &soft_bits) const
00498   {
00499     demodulate_soft_bits(rx_symbols, channel, N0, soft_bits);
00500   }
00501 
00502 
00503 
00504   //------------- class: PAM ----------------
00505 
00506   void PAM::modulate_bits(const bvec &bits, vec &out) const
00507   {
00508     // Check if some bits have to be cut and print warning message in
00509     // such case.
00510     if (bits.length() % k) {
00511       it_warning("PAM::modulate_bits(): The number of input bits is not a multiple of log2(M), where M is a constellation size. Remainder bits are not modulated.");
00512     }
00513     int no_symbols = bits.length() / k;
00514 
00515     out.set_size(no_symbols, false);
00516     for (int i = 0; i < no_symbols; i++) {
00517       out(i) = symbols(bits2symbols(bin2dec(bits.mid(i*k, k))));
00518     }
00519   }
00520 
00521   void PAM::modulate_bits(const bvec &bits, cvec &out) const
00522   {
00523     // Check if some bits have to be cut and print warning message in
00524     // such case.
00525     if (bits.length() % k) {
00526       it_warning("PAM::modulate_bits(): The number of input bits is not a multiple of log2(M), where M is a constellation size. Remainder bits are not modulated.");
00527     }
00528     int no_symbols = bits.length() / k;
00529 
00530     out.set_size(no_symbols, false);
00531     for (int i = 0; i < no_symbols; i++) {
00532       out(i) = symbols(bits2symbols(bin2dec(bits.mid(i*k, k))));
00533     }
00534   }
00535 
00536   cvec PAM::modulate_bits(const bvec &bits) const
00537   {
00538     cvec temp(bits.size());
00539     modulate_bits(bits, temp);
00540     return temp;
00541   }
00542 
00543   void PAM::demodulate_bits(const vec &signal, bvec &out) const
00544   {
00545     int est_symbol;
00546     out.set_size(k*signal.size(), false);
00547 
00548     for (int i = 0; i < signal.size(); i++) {
00549       est_symbol = round_i((M-1) - (signal(i) * scaling_factor + (M-1)) / 2);
00550       if (est_symbol < 0) 
00551         est_symbol = 0;
00552       else if (est_symbol > (M-1)) 
00553         est_symbol = M-1;
00554       out.replace_mid(i*k, bitmap.get_row(est_symbol));
00555     }
00556   }
00557 
00558   void PAM::demodulate_bits(const cvec &signal, bvec &out) const
00559   {
00560     int est_symbol;
00561     out.set_size(k*signal.size(), false);
00562 
00563     for (int i = 0; i < signal.size(); i++) {
00564       est_symbol = round_i((M-1) - (std::real(signal(i)) * scaling_factor 
00565                                     + (M-1)) / 2);
00566       if (est_symbol < 0)
00567         est_symbol = 0;
00568       else if (est_symbol > (M-1))
00569         est_symbol = M-1;
00570       out.replace_mid(i*k, bitmap.get_row(est_symbol));
00571     }
00572   }
00573 
00574   bvec PAM::demodulate_bits(const cvec &signal) const
00575   {
00576     bvec temp(signal.size());
00577     demodulate_bits(signal, temp);
00578     return temp;
00579   }
00580 
00581 
00582   void PAM::demodulate_soft_bits(const cvec &rx_symbols, double N0,
00583                                  vec &soft_bits) const
00584   {
00585     double P0, P1;
00586     vec expd(M);
00587 
00588     soft_bits.set_size(k*rx_symbols.size(), false);
00589 
00590     for (int l = 0; l < rx_symbols.size(); l++) {
00591       // calculate exponent of metric for each constellation point
00592       for (int j = 0; j < M; j++) {
00593         expd(j) = std::exp(-sqr(std::real(rx_symbols(l) - symbols(j))) / N0);
00594       }
00595 
00596       for (int i = 0; i < k; i++) {
00597         P0 = 0;
00598         P1 = 0;
00599         for (int j = 0; j < (M >> 1); j++) {
00600           P0 += expd(S0(i, j));
00601           P1 += expd(S1(i, j)); 
00602         }
00603         soft_bits(l*k+i) = trunc_log(P0) - trunc_log(P1);
00604       }
00605     }
00606   }
00607 
00608   void PAM::demodulate_soft_bits_approx(const cvec &rx_symbols, double N0,
00609                                         vec &soft_bits) const
00610   {
00611     double d0min, d1min;
00612     vec d(M);
00613 
00614     soft_bits.set_size(k*rx_symbols.size(), false);
00615 
00616     for (int l = 0; l < rx_symbols.size(); l++) {
00617       // calculate exponent of metric for each constellation point
00618       for (int j = 0; j < M; j++) {
00619         d(j) = sqr(std::real(rx_symbols(l) - symbols(j)));
00620       }
00621 
00622       for (int i = 0; i < k; i++) {
00623         d0min = std::numeric_limits<double>::max();
00624         d1min = std::numeric_limits<double>::max();
00625         for (int j = 0; j < (M >> 1); j++) {
00626           if (d(S0(i, j)) < d0min) { d0min = d(S0(i, j)); }
00627           if (d(S1(i, j)) < d1min) { d1min = d(S1(i, j)); }
00628         }
00629         soft_bits(l*k+i) = (-d0min + d1min) / N0;
00630       }
00631     }
00632   }
00633 
00634   void PAM::demodulate_soft_bits(const cvec &rx_symbols, const cvec &channel,
00635                                  double N0, vec &soft_bits) const
00636   {
00637     double P0, P1;
00638     vec expd(M);
00639 
00640     soft_bits.set_size(k*rx_symbols.size(), false);
00641 
00642     for (int l = 0; l < rx_symbols.size(); l++) {
00643       // calculate exponent of metric for each constellation point
00644       for (int j = 0; j < M; j++) {
00645         expd(j) = std::exp(-sqr(std::real(rx_symbols(l) 
00646                                           - channel(l) * symbols(j))) / N0);
00647       }
00648 
00649       for (int i = 0; i < k; i++) {
00650         P0 = 0;
00651         P1 = 0;
00652         for (int j = 0; j < (M >> 1); j++) {
00653           P0 += expd(S0(i, j));
00654           P1 += expd(S1(i, j)); 
00655         }
00656         soft_bits(l*k+i) = trunc_log(P0) - trunc_log(P1);
00657       }
00658     }
00659   }
00660 
00661   void PAM::demodulate_soft_bits_approx(const cvec &rx_symbols, 
00662                                         const cvec &channel, double N0,
00663                                         vec &soft_bits) const
00664   {
00665     double d0min, d1min;
00666     vec d(M);
00667 
00668     soft_bits.set_size(k*rx_symbols.size(), false);
00669 
00670     for (int l = 0; l < rx_symbols.size(); l++) {
00671       // calculate exponent of metric for each constellation point
00672       for (int j = 0; j < M; j++) {
00673         d(j) = sqr(std::real(rx_symbols(l) - channel(l) * symbols(j)));
00674       }
00675 
00676       for (int i = 0; i < k; i++) {
00677         d0min = std::numeric_limits<double>::max();
00678         d1min = std::numeric_limits<double>::max();
00679         for (int j = 0; j < (M >> 1); j++) {
00680           if (d(S0(i, j)) < d0min) { d0min = d(S0(i, j)); }
00681           if (d(S1(i, j)) < d1min) { d1min = d(S1(i, j)); }
00682         }
00683         soft_bits(l*k+i) = (-d0min + d1min) / N0;
00684       }
00685     }
00686   }
00687 
00688   void PAM::set_M(int Mary)
00689   {
00690     M = Mary;
00691     k = round_i(log2(double(M)));
00692 
00693     it_assert(pow2i(k) == Mary, "PAM::set_M(): M is not a power of 2");
00694 
00695     int count0, count1;
00696     bvec bits;
00697 
00698     symbols.set_size(M, false);
00699     bits2symbols.set_size(M, false);
00700     bitmap = graycode(k);
00701     average_energy = (sqr(M) - 1) / 3.0;
00702     scaling_factor = std::sqrt(average_energy);
00703 
00704     for (int i = 0; i < M; i++) {
00705       symbols(i) = ((M-1) - i*2) / scaling_factor;
00706       bits2symbols(bin2dec(bitmap.get_row(i))) = i;
00707     }
00708 
00709     // Calculate the soft bit mapping matrices S0 and S1
00710     S0.set_size(k, M/2, false);
00711     S1.set_size(k, M/2, false);
00712   
00713     for (int kk = 0; kk < k; kk++) {
00714       count0 = 0; 
00715       count1 = 0;
00716       for (int m = 0; m < M; m++) {
00717         bits = bitmap.get_row(m);
00718         if (bits(kk) == bin(0)) {
00719           S0(kk, count0) = m;
00720           count0++; 
00721         } else {
00722           S1(kk, count1) = m;
00723           count1++;
00724         }
00725       }
00726     }
00727   }
00728 
00729 
00730   //------------- class: QPSK ----------------
00731 
00732   void QPSK::modulate_bits(const bvec &bits, cvec &out) const
00733   {
00734     // Check if some bits have to be cut and print warning message in such
00735     // case.
00736     if (bits.length() & 0x1) {
00737       it_warning("QPSK::modulate_bits(): The number of input bits is not a multiple of 2. Remainder bits are not modulated.");
00738     }
00739     double real_part, imag_part;
00740     int no_symbols = bits.length() >> 1;
00741 
00742     out.set_size(no_symbols, false);
00743     for (int i = 0; i < no_symbols; i++) {
00744       real_part = (bits(2*i) == 0) ? M_SQRT1_2 : -M_SQRT1_2;
00745       imag_part = (bits(2*i+1) == 0) ? M_SQRT1_2 : -M_SQRT1_2;
00746       out(i) = std::complex<double>(real_part, imag_part);
00747     }
00748   }
00749 
00750   cvec QPSK::modulate_bits(const bvec &bits) const
00751   {
00752     cvec out;
00753     modulate_bits(bits, out);
00754     return out;
00755   }
00756 
00757   void QPSK::demodulate_bits(const cvec &signal, bvec &out) const
00758   {
00759     int no_symbols = signal.size();
00760     out.set_size(2*no_symbols, false);
00761     for (int i = 0; i < no_symbols; i++) {
00762       out(2*i) = ((std::real(signal(i)) > 0) ? 0 : 1);
00763       out(2*i+1) = ((std::imag(signal(i)) > 0) ? 0 : 1);
00764     }
00765   }
00766 
00767   bvec QPSK::demodulate_bits(const cvec &signal) const
00768   {
00769     bvec out;
00770     demodulate_bits(signal, out);
00771     return out;
00772   }
00773 
00774 
00775   // Outputs log-likelihood ratio of log (Pr(b=0|rx_symbols)/Pr(b=1|rx_symbols))
00776   void QPSK::demodulate_soft_bits(const cvec &rx_symbols, double N0,
00777                                   vec &soft_bits) const
00778   {
00779     soft_bits.set_size(2*rx_symbols.size(), false);
00780     double factor = 2 * std::sqrt(2.0) / N0;
00781     for (int i = 0; i < rx_symbols.size(); i++) {
00782       soft_bits(2*i) = std::real(rx_symbols(i)) * factor;
00783       soft_bits(2*i+1) = std::imag(rx_symbols(i)) * factor;
00784     }
00785   }
00786 
00787   // Outputs log-likelihood ratio for fading channels
00788   void QPSK::demodulate_soft_bits(const cvec &rx_symbols, const cvec &channel,
00789                                   double N0, vec &soft_bits) const
00790   {
00791     soft_bits.set_size(2*rx_symbols.size(), false);
00792     std::complex<double> temp;
00793     double factor = 2 * std::sqrt(2.0) / N0;
00794     
00795     for (int i = 0; i < rx_symbols.size(); i++) {
00796       temp = rx_symbols(i) * std::conj(channel(i));
00797       soft_bits(2*i) = std::real(temp) * factor;
00798       soft_bits(2*i+1) = std::imag(temp) * factor;
00799     }
00800   }
00801 
00802 
00803   void QPSK::demodulate_soft_bits_approx(const cvec &rx_symbols, double N0,
00804                                          vec &soft_bits) const
00805   {
00806     demodulate_soft_bits(rx_symbols, N0, soft_bits);
00807   }
00808 
00809   void QPSK::demodulate_soft_bits_approx(const cvec &rx_symbols,
00810                                          const cvec &channel, double N0,
00811                                          vec &soft_bits) const
00812   {
00813     demodulate_soft_bits(rx_symbols, channel, N0, soft_bits);
00814   }
00815 
00816 
00817   //------------- class: PSK ----------------
00818 
00819   void PSK::modulate_bits(const bvec &bits, cvec &out) const
00820   {
00821     // Check if some bits have to be cut and print warning message in such
00822     // case.
00823     if (bits.length() % k) {
00824       it_warning("PSK::modulate_bits(): The number of input bits is not a multiple of log2(M), where M is a constellation size. Remainder bits are not modulated.");
00825     }
00826     int no_symbols = bits.length() / k;
00827 
00828     out.set_size(no_symbols, false);
00829 
00830     for (int i = 0; i < no_symbols; i++) {
00831       out(i) = symbols(bits2symbols(bin2dec(bits.mid(i*k, k))));
00832     }
00833   }
00834 
00835   cvec PSK::modulate_bits(const bvec &bits) const
00836   {
00837     cvec temp(bits.size());
00838     modulate_bits(bits, temp);
00839     return temp;
00840   }
00841 
00842   void PSK::demodulate_bits(const cvec &signal, bvec &out) const
00843   {
00844     int est_symbol;
00845     double ang, temp;
00846 
00847     out.set_size(k*signal.size(), false);
00848 
00849     for (int i = 0; i < signal.size(); i++) {
00850       ang = std::arg(signal(i));
00851       temp = (ang < 0) ? (2*pi+ang) : ang;
00852       est_symbol = round_i(temp * (M/2) / pi) % M;
00853       out.replace_mid(i*k, bitmap.get_row(est_symbol));
00854     }
00855   }
00856 
00857   bvec PSK::demodulate_bits(const cvec &signal) const
00858   {
00859     bvec temp;
00860     demodulate_bits(signal, temp);
00861     return temp;
00862   }
00863 
00864   void PSK::demodulate_soft_bits(const cvec &rx_symbols, double N0,
00865                                  vec &soft_bits) const
00866   {
00867     double P0, P1;
00868     vec expd(M);
00869 
00870     soft_bits.set_size(k*rx_symbols.size(), false);
00871 
00872     for (int l = 0; l < rx_symbols.size(); l++) {
00873       for (int j = 0; j < M; j++) {
00874         expd(j) = std::exp(-sqr(rx_symbols(l) - symbols(j)) / N0);
00875       }
00876       for (int i = 0; i < k; i++) {
00877         P0 = 0;
00878         P1 = 0;
00879         for (int j = 0; j < (M >> 1); j++) {
00880           P0 += expd(S0(i, j));
00881           P1 += expd(S1(i, j));  
00882         }
00883         soft_bits(l*k+i) = trunc_log(P0) - trunc_log(P1);
00884       }
00885     }
00886   }
00887 
00888   void PSK::demodulate_soft_bits_approx(const cvec &rx_symbols, double N0,
00889                                         vec &soft_bits) const
00890   {
00891     double d0min, d1min;
00892     vec d(M);
00893 
00894     soft_bits.set_size(k*rx_symbols.size(), false);
00895 
00896     for (int l = 0; l < rx_symbols.size(); l++) {
00897       for (int j = 0; j < M; j++)
00898         d(j) = sqr(rx_symbols(l) - symbols(j));
00899       for (int i = 0; i < k; i++) {
00900         d0min = std::numeric_limits<double>::max();
00901         d1min = std::numeric_limits<double>::max();
00902         for (int j = 0; j < (M >> 1); j++) {
00903           if (d(S0(i, j)) < d0min) { d0min = d(S0(i, j)); }
00904           if (d(S1(i, j)) < d1min) { d1min = d(S1(i, j)); }
00905         }
00906         soft_bits(l*k+i) = (-d0min + d1min) / N0;
00907       }
00908     }
00909         }
00910 
00911   void PSK::demodulate_soft_bits(const cvec &rx_symbols, const cvec &channel,
00912                                  double N0, vec &soft_bits) const
00913   {
00914     double P0, P1;
00915     vec expd(M);
00916 
00917     soft_bits.set_size(k*rx_symbols.size(), false);
00918 
00919     for (int l = 0; l < rx_symbols.size(); l++) {
00920       for (int j = 0; j < M; j++) {
00921         expd(j) = std::exp(-sqr(rx_symbols(l) - channel(l) * symbols(j)) / N0);
00922       }
00923       for (int i = 0; i < k; i++) {
00924         P0 = 0;
00925         P1 = 0;
00926         for (int j = 0; j < (M >> 1); j++) {
00927           P0 += expd(S0(i, j));
00928           P1 += expd(S1(i, j));  
00929         }
00930         soft_bits(l*k+i) = trunc_log(P0) - trunc_log(P1);
00931       }
00932     }
00933   }
00934 
00935   void PSK::demodulate_soft_bits_approx(const cvec &rx_symbols,
00936                                         const cvec &channel, double N0,
00937                                         vec &soft_bits) const
00938   {
00939     double d0min, d1min;
00940     vec d(M);
00941 
00942     soft_bits.set_size(k*rx_symbols.size(), false);
00943 
00944     for (int l = 0; l < rx_symbols.size(); l++) {
00945       for (int j = 0; j < M; j++)
00946         d(j) = sqr(rx_symbols(l) - channel(l) * symbols(j));
00947       for (int i = 0; i < k; i++) {
00948         d0min = std::numeric_limits<double>::max();
00949         d1min = std::numeric_limits<double>::max();
00950         for (int j = 0; j < (M >> 1); j++) {
00951           if (d(S0(i, j)) < d0min) { d0min = d(S0(i, j)); }
00952           if (d(S1(i, j)) < d1min) { d1min = d(S1(i, j)); }
00953         }
00954         soft_bits(l*k+i) = (-d0min + d1min) / N0;
00955       }
00956     }
00957   }
00958 
00959   void PSK::set_M(int Mary)
00960   {
00961     k = round_i(log2(double(Mary)));
00962     M = Mary;
00963     it_assert(pow2i(k) == Mary, "PSK::set_M(): M is not a power of 2");
00964     symbols.set_size(M, false);
00965     bits2symbols.set_size(M, false);
00966     bitmap = graycode(k);
00967 
00968     double delta = 2.0 * pi / M;
00969     double epsilon = delta / 10000.0;
00970     std::complex<double> symb;
00971     for (int i = 0; i < M; i++) {
00972       symb = std::complex<double>(std::polar(1.0, delta*i));
00973       if (std::fabs(std::real(symb)) < epsilon) { 
00974         symbols(i) = std::complex<double>(0.0, std::imag(symb)); 
00975       }
00976       else if (std::fabs(std::imag(symb)) < epsilon) { 
00977         symbols(i) = std::complex<double>(std::real(symb), 0.0);
00978       }
00979       else { symbols(i) = symb; }
00980 
00981       bits2symbols(bin2dec(bitmap.get_row(i))) = i;
00982     }
00983 
00984     // Calculate the soft bit mapping matrices S0 and S1
00985     S0.set_size(k, M/2,false);
00986     S1.set_size(k, M/2,false);
00987     int count0, count1;
00988     bvec bits;
00989   
00990     for (int kk = 0; kk < k; kk++) {
00991       count0 = 0; 
00992       count1 = 0;
00993       for (int m = 0; m < M; m++) {
00994         bits = bitmap.get_row(m);
00995         if (bits(kk) == bin(0)) {
00996           S0(kk, count0) = m;
00997           count0++; 
00998         } else {
00999           S1(kk, count1) = m;
01000           count1++;
01001         }
01002       }
01003     }
01004   }
01005 
01006 
01007   //------------- class: QAM ----------------
01008 
01009   void QAM::modulate_bits(const bvec &bits, cvec &out) const
01010   {
01011     // Check if some bits have to be cut and print warning message in such
01012     // case.
01013     if (bits.length() % k) {
01014       it_warning("QAM::modulate_bits(): The number of input bits is not a multiple of log2(M), where M is a constellation size. Remainder bits are not modulated.");
01015     }
01016     int no_symbols = bits.length() / k;
01017     out.set_size(no_symbols, false);
01018     for (int i = 0; i < no_symbols; i++) {
01019       out(i) = symbols(bits2symbols(bin2dec(bits.mid(i*k, k))));
01020     }
01021   }
01022 
01023   cvec QAM::modulate_bits(const bvec &bits) const
01024   {
01025     cvec temp(bits.size());
01026     modulate_bits(bits, temp);
01027     return temp;
01028   }
01029 
01030   void QAM::demodulate_bits(const cvec &signal, bvec &out) const
01031   {
01032     out.set_size(k*signal.size(), false);
01033 
01034     int temp_real, temp_imag;
01035 
01036     for (int i = 0; i < signal.size(); i++) {
01037       temp_real = round_i((L-1) - (std::real(signal(i) * scaling_factor)
01038                                    + (L-1)) / 2.0);
01039       temp_imag = round_i((L-1) - (std::imag(signal(i) * scaling_factor)
01040                                    + (L-1)) / 2.0);
01041       if (temp_real < 0) 
01042         temp_real = 0; 
01043       else if (temp_real > (L-1)) 
01044         temp_real = (L-1);
01045       if (temp_imag < 0) 
01046         temp_imag = 0; 
01047       else if (temp_imag > (L-1)) 
01048         temp_imag = (L-1);
01049       out.replace_mid(k*i, bitmap.get_row(temp_imag * L + temp_real));
01050     }
01051   }
01052 
01053   bvec QAM::demodulate_bits(const cvec &signal) const
01054   {
01055     bvec temp;
01056     demodulate_bits(signal, temp);
01057     return temp;
01058   }
01059 
01060   void QAM::demodulate_soft_bits(const cvec &rx_symbols, double N0,
01061                                  vec &soft_bits) const
01062   {
01063     double P0, P1;
01064     vec expd(M);
01065 
01066     soft_bits.set_size(k*rx_symbols.size(), false);
01067 
01068     for (int l = 0; l < rx_symbols.size(); l++) {
01069       for (int j = 0; j < M; j++) {
01070         expd(j) = std::exp(-sqr(rx_symbols(l) - symbols(j)) / N0);
01071       }
01072       for (int i = 0; i < k; i++) {
01073         P0 = 0;
01074         P1 = 0;
01075         for (int j = 0; j < (M >> 1); j++) {
01076           P0 += expd(S0(i, j));
01077           P1 += expd(S1(i, j));  
01078         }
01079         soft_bits(l*k+i) = trunc_log(P0) - trunc_log(P1);
01080       }
01081     }
01082   }
01083 
01084 
01085   void QAM::demodulate_soft_bits_approx(const cvec &rx_symbols, double N0,
01086                                         vec &soft_bits) const
01087   {
01088     double d0min, d1min;
01089     vec d(M);
01090 
01091     soft_bits.set_size(k*rx_symbols.size(), false);
01092 
01093     for (int l = 0; l < rx_symbols.size(); l++) {
01094       for (int j = 0; j < M; j++)
01095         d(j) = sqr(rx_symbols(l) - symbols(j));
01096 
01097       for (int i = 0; i < k; i++) {
01098         d0min = std::numeric_limits<double>::max();
01099         d1min = std::numeric_limits<double>::max();
01100         for (int j = 0; j < (M >> 1); j++) {
01101           if (d(S0(i, j)) < d0min) { d0min = d(S0(i, j)); }
01102           if (d(S1(i, j)) < d1min) { d1min = d(S1(i, j)); }
01103         }
01104         soft_bits(l*k+i) = (-d0min + d1min) / N0;
01105       }
01106     }
01107   }
01108 
01109   void QAM::demodulate_soft_bits(const cvec &rx_symbols, const cvec &channel,
01110                                  double N0, vec &soft_bits) const
01111   {
01112     double P0, P1;
01113     vec expd(M);
01114 
01115     soft_bits.set_size(k*rx_symbols.size(), false);
01116 
01117     for (int l = 0; l < rx_symbols.size(); l++) {
01118       for (int j = 0; j < M; j++) {
01119         expd(j) = std::exp(-sqr(rx_symbols(l) - channel(l) * symbols(j)) / N0);
01120       }
01121       for (int i = 0; i < k; i++) {
01122         P0 = 0;
01123         P1 = 0;
01124         for (int j = 0; j < (M >> 1); j++) {
01125           P0 += expd(S0(i, j));
01126           P1 += expd(S1(i, j));  
01127         }
01128         soft_bits(l*k+i) = trunc_log(P0) - trunc_log(P1);
01129       }
01130     }
01131   }
01132 
01133 
01134   void QAM::demodulate_soft_bits_approx(const cvec &rx_symbols, 
01135                                         const cvec &channel, double N0,
01136                                         vec &soft_bits) const
01137   {
01138     double d0min, d1min;
01139     vec d(M);
01140 
01141     soft_bits.set_size(k*rx_symbols.size(), false);
01142 
01143     for (int l = 0; l < rx_symbols.size(); l++) {
01144       for (int j = 0; j < M; j++)
01145         d(j) = sqr(rx_symbols(l) - channel(l) * symbols(j));
01146 
01147       for (int i = 0; i < k; i++) {
01148         d0min = std::numeric_limits<double>::max();
01149         d1min = std::numeric_limits<double>::max();
01150         for (int j = 0; j < (M >> 1); j++) {
01151           if (d(S0(i, j)) < d0min) { d0min = d(S0(i, j)); }
01152           if (d(S1(i, j)) < d1min) { d1min = d(S1(i, j)); }
01153         }
01154         soft_bits(l*k+i) = (-d0min + d1min) / N0;
01155       }
01156     }
01157   }
01158 
01159 
01160   void QAM::set_M(int Mary)
01161   {
01162     k = round_i(log2(double(Mary)));
01163     M = Mary;
01164     L = round_i(std::sqrt(static_cast<double>(M)));
01165     it_assert((pow2i(k) == Mary) && (is_even(k)), 
01166               "QAM::set_M(): M is not an even power of 2");
01167 
01168     int count0, count1;
01169     bvec bits;
01170 
01171     symbols.set_size(M, false);
01172     bitmap.set_size(M, k, false);
01173     bits2symbols.set_size(M, false);
01174     bmat gray_code=graycode(round_i(log2(double(L))));
01175     average_energy = (M-1) * 2.0 / 3.0;
01176     scaling_factor = std::sqrt(average_energy);
01177 
01178     for (int i = 0; i < L; i++) {
01179       for (int j = 0; j < L; j++) {
01180         symbols(i*L+j) = std::complex<double>(((L-1) - j*2) / scaling_factor,
01181                                               ((L-1) - i*2) / scaling_factor);
01182         bitmap.set_row(i*L+j, concat(gray_code.get_row(i), 
01183                                      gray_code.get_row(j)));
01184         bits2symbols(bin2dec(bitmap.get_row(i*L+j))) = i*L+j;
01185       }
01186     }
01187 
01188     // Calculate the soft bit mapping matrices S0 and S1
01189     S0.set_size(k, M/2, false);
01190     S1.set_size(k, M/2, false);
01191   
01192     for (int kk = 0; kk < k; kk++) {
01193       count0 = 0; 
01194       count1 = 0;
01195       for (int m = 0; m < M; m++) {
01196         bits = bitmap.get_row(m);
01197         if (bits(kk) == bin(0)) {
01198           S0(kk, count0) = m;
01199           count0++; 
01200         } else {
01201           S1(kk, count1) = m;
01202           count1++;
01203         }
01204       }
01205     }
01206   }
01207 
01208 } // namespace itpp
SourceForge Logo

Generated on Wed Apr 18 11:19:59 2007 for IT++ by Doxygen 1.5.2