00001 00002 #include <itpp/itcomm.h> 00003 00004 using std::cout; 00005 using std::endl; 00006 using namespace itpp; 00007 using namespace std; 00008 00009 /* Zero-forcing detector with ad hoc soft information. This function 00010 applies the ZF (pseudoinverse) linear filter to the received 00011 data. This results in effective noise with covariance matrix 00012 inv(H'H)*sigma2. The diagonal elements of this noise covariance 00013 matrix are taken as noise variances per component in the processed 00014 received data but the noise correlation is ignored. 00015 00016 */ 00017 00018 void ZF_demod(ND_UQAM &channel, ivec &LLR_apr, ivec &LLR_apost, double sigma2, cmat &H, cvec &y) 00019 { 00020 it_assert(H.rows()>=H.cols(),"ZF_demod() - underdetermined systems not tolerated"); 00021 cvec shat=ls_solve_od(H,y); // the ZF solution 00022 vec Sigma2=real(diag(inv(H.hermitian_transpose()*H)))*sigma2; // noise covariance of shat 00023 cvec h(length(shat)); 00024 for (int i=0; i<length(shat); i++) { 00025 shat(i) = shat(i)/sqrt(Sigma2(i)); 00026 h(i) = 1.0/sqrt(Sigma2(i)); 00027 } 00028 channel.map_demod(LLR_apr,LLR_apost,1.0,h,shat); 00029 } 00030 00031 00032 extern int main(int argc, char **argv) 00033 { 00034 // -- modulation and channel parameters (taken from command line input) -- 00035 int nC; // type of constellation (1=QPSK, 2=16-QAM, 3=64-QAM) 00036 int nRx; // number of receive antennas 00037 int nTx; // number of transmit antennas 00038 int Tc; // coherence time (number of channel vectors with same H) 00039 00040 if (argc!=5) { 00041 cout << "Usage: cm nTx nRx nC Tc" << endl << "Example: cm 2 2 1 100000 (2x2 QPSK MIMO on slow fading channel)" << endl; 00042 exit(1); 00043 } else { 00044 sscanf(argv[1],"%i",&nTx); 00045 sscanf(argv[2],"%i",&nRx); 00046 sscanf(argv[3],"%i",&nC); 00047 sscanf(argv[4],"%i",&Tc); 00048 } 00049 00050 cout << "Initializing.. " << nTx << " TX antennas, " << nRx << " RX antennas, " 00051 << (1<<nC) << "-PAM per dimension, coherence time " << Tc << endl; 00052 00053 // -- simulation control parameters -- 00054 const vec EbN0db = "-5:0.5:50"; // SNR range 00055 const int Nmethods =2; // number of demodulators to try 00056 const long int Nbitsmax=50*1000*1000; // maximum number of bits to ever simulate per SNR point 00057 const int Nu = 1000; // length of data packet (before applying channel coding) 00058 00059 int Nbers, Nfers; // target number of bit/frame errors per SNR point 00060 double BERmin, FERmin; // BER/FER at which to terminate simulation 00061 if (Tc==1) { // Fast fading channel, BER is of primary interest 00062 BERmin = 0.001; // stop simulating a given method if BER<this value 00063 FERmin = 1.0e-10; // stop simulating a given method if FER<this value 00064 Nbers = 1000; // move to next SNR point after counting 1000 bit errors 00065 Nfers = 200; // do not stop on this condition 00066 } else { // Slow fading channel, FER is of primary interest here 00067 BERmin = 1.0e-15; // stop simulating a given method if BER<this value 00068 FERmin = 0.01; // stop simulating a given method if FER<this value 00069 Nbers = -1; // do not stop on this condition 00070 Nfers = 200; // move to next SNR point after counting 200 frame errors 00071 } 00072 00073 // -- Channel code parameters -- 00074 Convolutional_Code code; 00075 ivec generator(3); 00076 generator(0)=0133; // use rate 1/3 code 00077 generator(1)=0165; 00078 generator(2)=0171; 00079 double rate=1.0/3.0; 00080 code.set_generator_polynomials(generator, 7); 00081 bvec dummy; 00082 code.encode_tail(randb(Nu),dummy); 00083 const int Nc = length(dummy); // find out how long the coded blocks are 00084 00085 // ============= Initialize ==================================== 00086 00087 const int Nctx = (int) (2*nC*nTx*ceil(double(Nc)/double(2*nC*nTx))); // Total number of bits to transmit 00088 const int Nvec = Nctx/(2*nC*nTx); // Number of channel vectors to transmit 00089 const int Nbitspvec = 2*nC*nTx; // Number of bits per channel vector 00090 00091 // initialize MIMO channel with uniform QAM per complex dimension and Gray coding 00092 ND_UQAM chan; 00093 chan.set_Gray_QAM(nTx,1<<(2*nC)); 00094 cout << chan << endl; 00095 00096 // initialize interleaver 00097 Sequence_Interleaver<bin> sequence_interleaver_b(Nctx); 00098 Sequence_Interleaver<int> sequence_interleaver_i(Nctx); 00099 sequence_interleaver_b.randomize_interleaver_sequence(); 00100 sequence_interleaver_i.set_interleaver_sequence(sequence_interleaver_b.get_interleaver_sequence()); 00101 00102 // RNG_randomize(); 00103 00104 Array<cvec> Y(Nvec); // received data 00105 Array<cmat> H(Nvec/Tc+1); // channel matrix (new matrix for each coherence interval) 00106 00107 ivec Contflag = ones_i(Nmethods); // flag to determine whether to run a given demodulator 00108 if (pow(2.0,nC*2.0*nTx)>256) { // ML decoder too complex.. 00109 Contflag(1)=0; 00110 } 00111 if (nTx>nRx) { 00112 Contflag(0)=0; // ZF not for underdetermined systems 00113 } 00114 cout << "Running methods: " << Contflag << endl; 00115 00116 cout.setf(ios::fixed, ios::floatfield); 00117 cout.setf(ios::showpoint); 00118 cout.precision(5); 00119 00120 // ================== Run simulation ======================= 00121 for (int nsnr=0; nsnr<length(EbN0db); nsnr++) { 00122 const double Eb=1.0; // transmitted energy per information bit 00123 const double N0 = pow(10.0,-EbN0db(nsnr)/10.0); 00124 const double sigma2=N0; // Variance of each scalar complex noise sample 00125 const double Es=rate*2*nC*Eb; // Energy per complex scalar symbol 00126 // (Each transmitted scalar complex symbol contains rate*2*nC 00127 // information bits.) 00128 const double Ess=sqrt(Es); 00129 00130 Array<BERC> berc(Nmethods); // counter for coded BER 00131 Array<BERC> bercu(Nmethods); // counter for uncoded BER 00132 Array<BLERC> ferc(Nmethods); // counter for coded FER 00133 00134 for (int i=0; i<Nmethods; i++) { 00135 ferc(i).set_blocksize(Nu); 00136 } 00137 00138 long int nbits=0; 00139 while (nbits<Nbitsmax) { 00140 nbits += Nu; 00141 00142 // generate and encode random data 00143 bvec inputbits = randb(Nu); 00144 bvec txbits; 00145 code.encode_tail(inputbits, txbits); 00146 // coded block length is not always a multiple of the number of 00147 // bits per channel vector 00148 txbits=concat(txbits,randb(Nctx-Nc)); 00149 txbits = sequence_interleaver_b.interleave(txbits); 00150 00151 // -- generate channel and data ---- 00152 for (int k=0; k<Nvec; k++) { 00153 /* A complex valued channel matrix is used here. An 00154 alternative (with equivalent result) would be to use a 00155 real-valued (structured) channel matrix of twice the 00156 dimension. 00157 */ 00158 if (k%Tc==0) { // generate a new channel realization every Tc intervals 00159 H(k/Tc) = Ess*randn_c(nRx,nTx); 00160 } 00161 00162 // modulate and transmit bits 00163 bvec bitstmp = txbits(k*2*nTx*nC,(k+1)*2*nTx*nC-1); 00164 cvec x=chan.modulate_bits(bitstmp); 00165 cvec e=sqrt(sigma2)*randn_c(nRx); 00166 Y(k) = H(k/Tc)*x+e; 00167 } 00168 00169 // -- demodulate -- 00170 Array<QLLRvec> LLRin(Nmethods); 00171 for (int i=0; i<Nmethods; i++) { 00172 LLRin(i) = zeros_i(Nctx); 00173 } 00174 00175 QLLRvec llr_apr =zeros_i(nC*2*nTx); // no a priori input to demodulator 00176 QLLRvec llr_apost=zeros_i(nC*2*nTx); 00177 for (int k=0; k<Nvec; k++) { 00178 // zero forcing demodulation 00179 if (Contflag(0)) { 00180 ZF_demod(chan,llr_apr,llr_apost,sigma2,H(k/Tc),Y(k)); 00181 LLRin(0).set_subvector(k*Nbitspvec,(k+1)*Nbitspvec-1,llr_apost); 00182 } 00183 00184 // ML demodulation 00185 if (Contflag(1)) { 00186 chan.map_demod(llr_apr, llr_apost, sigma2, H(k/Tc), Y(k)); 00187 LLRin(1).set_subvector(k*Nbitspvec,(k+1)*Nbitspvec-1,llr_apost); 00188 } 00189 } 00190 00191 // -- decode and count errors -- 00192 for (int i=0; i<Nmethods; i++) { 00193 bvec decoded_bits; 00194 if (Contflag(i)) { 00195 bercu(i).count(txbits(0,Nc-1),LLRin(i)(0,Nc-1)<0); // uncoded BER 00196 LLRin(i) = sequence_interleaver_i.deinterleave(LLRin(i),0); 00197 // QLLR values must be converted to real numbers since the convolutional decoder wants this 00198 vec llr=chan.get_llrcalc().to_double(LLRin(i).left(Nc)); 00199 // llr=-llr; // UNCOMMENT THIS LINE IF COMPILING WITH 3.10.5 OR EARLIER (BEFORE HARMONIZING LLR CONVENTIONS) 00200 code.decode_tail(llr,decoded_bits); 00201 berc(i).count(inputbits(0,Nu-1),decoded_bits(0,Nu-1)); // coded BER 00202 ferc(i).count(inputbits(0,Nu-1),decoded_bits(0,Nu-1)); // coded FER 00203 } 00204 } 00205 00206 /* Check whether it is time to terminate the simulation. 00207 Terminate when all demodulators that are still running have 00208 counted at least Nbers or Nfers bit/frame errors. */ 00209 int minber=1000000; 00210 int minfer=1000000; 00211 for (int i=0; i<Nmethods; i++) { 00212 if (Contflag(i)) { 00213 minber=min(minber,round_i(berc(i).get_errors())); 00214 minfer=min(minfer,round_i(ferc(i).get_errors())); 00215 } 00216 } 00217 if (Nbers>0 && minber>Nbers) { break;} 00218 if (Nfers>0 && minfer>Nfers) { break;} 00219 } 00220 00221 cout << "-----------------------------------------------------" << endl; 00222 cout << "Eb/N0: " << EbN0db(nsnr) << " dB. Simulated " << nbits << " bits." << endl; 00223 cout << " Uncoded BER: " << bercu(0).get_errorrate() << " (ZF); " << bercu(1).get_errorrate() << " (ML)" << endl; 00224 cout << " Coded BER: " << berc(0).get_errorrate() << " (ZF); " << berc(1).get_errorrate() << " (ML)" << endl; 00225 cout << " Coded FER: " << ferc(0).get_errorrate() << " (ZF); " << ferc(1).get_errorrate() << " (ML)" << endl; 00226 cout.flush(); 00227 00228 /* Check wheter it is time to terminate simulation. Stop when all 00229 methods have reached the min BER/FER of interest. */ 00230 int contflag=0; 00231 for (int i=0; i<Nmethods; i++) { 00232 if (Contflag(i)) { 00233 if (berc(i).get_errorrate()>BERmin) { contflag=1; } else { Contflag(i)= 0; } 00234 if (ferc(i).get_errorrate()>FERmin) { contflag=1; } else { Contflag(i)= 0; } 00235 } 00236 } 00237 if (contflag) { continue; } else {break; } 00238 } 00239 00240 return 0; 00241 }
Generated on Fri Sep 11 06:31:59 2009 for RMOL by Doxygen 1.5.8