RMOL Logo Get Revenue Management Optimisation Library at SourceForge.net. Fast, secure and Free Open Source software downloads

mimoconv.cpp

Go to the documentation of this file.
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 }
SourceForge Logo

Generated on Tue Apr 14 17:57:51 2009 for RMOL by Doxygen 1.5.8