mimoconv.cpp

这是fast fix-point algorithm 的C＋＋版本

#include <itpp/itcomm.h>

using std::cout;
using std::endl;
using namespace itpp;
using namespace std;

/* Zero-forcing detector with ad hoc soft information. This function
   applies the ZF (pseudoinverse) linear filter to the received
   data. This results in effective noise with covariance matrix
   inv(H'H)*sigma2.  The diagonal elements of this noise covariance
   matrix are taken as noise variances per component in the processed
   received data but the noise correlation is ignored.

 */

void ZF_demod(ND_UQAM &channel, ivec &LLR_apr,  ivec &LLR_apost, double sigma2, cmat &H, cvec &y)
{
  it_assert(H.rows()>=H.cols(),"ZF_demod() - underdetermined systems not tolerated");
  cvec shat=ls_solve_od(H,y);                     // the ZF solution
  vec Sigma2=real(diag(inv(H.hermitian_transpose()*H)))*sigma2;  // noise covariance of shat
  cvec h(length(shat));
  for (int i=0; i<length(shat); i++) {
    shat(i) = shat(i)/sqrt(Sigma2(i));
    h(i) = 1.0/sqrt(Sigma2(i));
  }
  channel.map_demod(LLR_apr,LLR_apost,1.0,h,shat);
}


extern int main(int argc, char **argv)
{ 
  // -- modulation and channel parameters (taken from command line input) --
  int nC;                    // type of constellation  (1=QPSK, 2=16-QAM, 3=64-QAM)
  int nRx;                   // number of receive antennas
  int nTx;                   // number of transmit antennas 
  int Tc;                    // coherence time (number of channel vectors with same H)

  if (argc!=5) {
    cout << "Usage: cm nTx nRx nC Tc" << endl << "Example: cm 2 2 1 100000 (2x2 QPSK MIMO on slow fading channel)" << endl;
    exit(1); 
  } else {
    sscanf(argv[1],"%i",&nTx);
    sscanf(argv[2],"%i",&nRx);
    sscanf(argv[3],"%i",&nC);
    sscanf(argv[4],"%i",&Tc);
  }

  cout << "Initializing.. " << nTx << " TX antennas, " << nRx << " RX antennas, " 
       << (1<<nC) << "-PAM per dimension, coherence time " << Tc << endl;

  // -- simulation control parameters --
  const vec EbN0db = "-5:0.5:50";        // SNR range
  const int Nmethods =2;                 // number of demodulators to try
  const long int Nbitsmax=50*1000*1000;  // maximum number of bits to ever simulate per SNR point
  const int Nu = 1000;                   // length of data packet (before applying channel coding)

  int Nbers, Nfers;              // target number of bit/frame errors per SNR point
  double BERmin, FERmin;         // BER/FER at which to terminate simulation
  if (Tc==1) {           // Fast fading channel, BER is of primary interest 
    BERmin = 0.001;      // stop simulating a given method if BER<this value
    FERmin = 1.0e-10;    // stop simulating a given method if FER<this value
    Nbers = 1000;        // move to next SNR point after counting 1000 bit errors
    Nfers = 200;         // do not stop on this condition
  } else {               // Slow fading channel, FER is of primary interest here
    BERmin = 1.0e-15;    // stop simulating a given method if BER<this value
    FERmin = 0.01;       // stop simulating a given method if FER<this value
    Nbers = -1;          // do not stop on this condition
    Nfers = 200;         // move to next SNR point after counting 200 frame errors
  }

  // -- Channel code parameters --
  Convolutional_Code code;
  ivec generator(3);
  generator(0)=0133;  // use rate 1/3 code
  generator(1)=0165;
  generator(2)=0171;
  double rate=1.0/3.0;
  code.set_generator_polynomials(generator, 7);
  bvec dummy;
  code.encode_tail(randb(Nu),dummy);
  const int Nc = length(dummy);      // find out how long the coded blocks are

  // ============= Initialize ====================================

  const int Nctx = (int) (2*nC*nTx*ceil(double(Nc)/double(2*nC*nTx)));   // Total number of bits to transmit   
  const int Nvec = Nctx/(2*nC*nTx);          // Number of channel vectors to transmit
  const int Nbitspvec = 2*nC*nTx;            // Number of bits per channel vector

  // initialize MIMO channel with uniform QAM per complex dimension and Gray coding
  ND_UQAM chan;
  chan.set_Gray_QAM(nTx,1<<(2*nC));  
  cout << chan << endl;

  // initialize interleaver
  Sequence_Interleaver<bin> sequence_interleaver_b(Nctx);
  Sequence_Interleaver<int> sequence_interleaver_i(Nctx);
  sequence_interleaver_b.randomize_interleaver_sequence();
  sequence_interleaver_i.set_interleaver_sequence(sequence_interleaver_b.get_interleaver_sequence());

  //  RNG_randomize();

  Array<cvec> Y(Nvec);        // received data
  Array<cmat> H(Nvec/Tc+1);   // channel matrix (new matrix for each coherence interval)

  ivec Contflag = ones_i(Nmethods);   // flag to determine whether to run a given demodulator
  if (pow(2.0,nC*2.0*nTx)>256) {      // ML decoder too complex..
    Contflag(1)=0;  
  }
  if (nTx>nRx) {
    Contflag(0)=0;                    // ZF not for underdetermined systems
  }
  cout << "Running methods: " << Contflag << endl;

  cout.setf(ios::fixed, ios::floatfield); 
  cout.setf(ios::showpoint); 
  cout.precision(5);    
  
  // ================== Run simulation =======================
  for (int nsnr=0; nsnr<length(EbN0db); nsnr++) {
    const double Eb=1.0; // transmitted energy per information bit
    const double N0 = pow(10.0,-EbN0db(nsnr)/10.0);
    const double sigma2=N0;   // Variance of each scalar complex noise sample
    const double Es=rate*2*nC*Eb; // Energy per complex scalar symbol
    // (Each transmitted scalar complex symbol contains rate*2*nC
    // information bits.)
    const double Ess=sqrt(Es);

    Array<BERC> berc(Nmethods);  // counter for coded BER
    Array<BERC> bercu(Nmethods); // counter for uncoded BER
    Array<BLERC> ferc(Nmethods); // counter for coded FER

    for (int i=0; i<Nmethods; i++) {
      ferc(i).set_blocksize(Nu);
    }

    long int nbits=0;
    while (nbits<Nbitsmax) {
      nbits += Nu;

      // generate and encode random data
      bvec inputbits = randb(Nu);
      bvec txbits;
      code.encode_tail(inputbits, txbits);
      // coded block length is not always a multiple of the number of
      // bits per channel vector
      txbits=concat(txbits,randb(Nctx-Nc));   
      txbits = sequence_interleaver_b.interleave(txbits);

      // -- generate channel and data ----
      for (int k=0; k<Nvec; k++) {
	/* A complex valued channel matrix is used here. An
	   alternative (with equivalent result) would be to use a
	   real-valued (structured) channel matrix of twice the
	   dimension.
	*/
	if (k%Tc==0) {       // generate a new channel realization every Tc intervals
	  H(k/Tc) = Ess*randn_c(nRx,nTx);
	}
	
	// modulate and transmit bits
	bvec bitstmp = txbits(k*2*nTx*nC,(k+1)*2*nTx*nC-1);
	cvec x=chan.modulate_bits(bitstmp);
	cvec e=sqrt(sigma2)*randn_c(nRx);  
	Y(k) = H(k/Tc)*x+e;
      }

      // -- demodulate --
      Array<QLLRvec> LLRin(Nmethods);     
      for (int i=0; i<Nmethods; i++) {
	LLRin(i) = zeros_i(Nctx);
      }

      QLLRvec llr_apr =zeros_i(nC*2*nTx);  // no a priori input to demodulator
      QLLRvec llr_apost=zeros_i(nC*2*nTx);
      for (int k=0; k<Nvec; k++) {                
	// zero forcing demodulation
	if (Contflag(0)) {
	  ZF_demod(chan,llr_apr,llr_apost,sigma2,H(k/Tc),Y(k));
	  LLRin(0).set_subvector(k*Nbitspvec,(k+1)*Nbitspvec-1,llr_apost);
	}
	  
	// ML demodulation
	if (Contflag(1)) { 
	  chan.map_demod(llr_apr, llr_apost, sigma2, H(k/Tc), Y(k));
	  LLRin(1).set_subvector(k*Nbitspvec,(k+1)*Nbitspvec-1,llr_apost);
	}	  
      }
            
      // -- decode and count errors --
      for (int i=0; i<Nmethods; i++) {
	bvec decoded_bits;
	if (Contflag(i)) {
	  bercu(i).count(txbits(0,Nc-1),LLRin(i)(0,Nc-1)<0);  // uncoded BER
	  LLRin(i) = sequence_interleaver_i.deinterleave(LLRin(i),0);
	  // QLLR values must be converted to real numbers since the convolutional decoder wants this
	  vec llr=chan.get_llrcalc().to_double(LLRin(i).left(Nc)); 
	  //	  llr=-llr; // UNCOMMENT THIS LINE IF COMPILING WITH 3.10.5 OR EARLIER (BEFORE HARMONIZING LLR CONVENTIONS)
	  code.decode_tail(llr,decoded_bits);
	  berc(i).count(inputbits(0,Nu-1),decoded_bits(0,Nu-1));  // coded BER
	  ferc(i).count(inputbits(0,Nu-1),decoded_bits(0,Nu-1));  // coded FER
	}
      }	
      
      /* Check whether it is time to terminate the simulation.
       Terminate when all demodulators that are still running have
       counted at least Nbers or Nfers bit/frame errors. */
      int minber=1000000;
      int minfer=1000000;
      for (int i=0; i<Nmethods; i++) {
	if (Contflag(i)) {
	  minber=min(minber,round_i(berc(i).get_errors()));   
	  minfer=min(minfer,round_i(ferc(i).get_errors()));  
	}
      }
      if (Nbers>0 && minber>Nbers) { break;}
      if (Nfers>0 && minfer>Nfers) { break;}
    }
    
    cout << "-----------------------------------------------------" << endl;
    cout << "Eb/N0: " << EbN0db(nsnr) << " dB. Simulated " << nbits << " bits." << endl;
    cout << " Uncoded BER: " << bercu(0).get_errorrate() << " (ZF);     " << bercu(1).get_errorrate() << " (ML)" << endl;
    cout << " Coded BER:   " << berc(0).get_errorrate()  << " (ZF);     " << berc(1).get_errorrate()  << " (ML)" << endl;
    cout << " Coded FER:   " << ferc(0).get_errorrate()  << " (ZF);     " << ferc(1).get_errorrate()  << " (ML)" << endl;
    cout.flush();

    /* Check wheter it is time to terminate simulation. Stop when all
    methods have reached the min BER/FER of interest. */
    int contflag=0;
    for (int i=0; i<Nmethods; i++) {
      if (Contflag(i)) {
	if (berc(i).get_errorrate()>BERmin)  {  contflag=1;  } else { Contflag(i)= 0; } 
	if (ferc(i).get_errorrate()>FERmin)  {  contflag=1;  } else { Contflag(i)= 0; } 
      }
    }
    if (contflag) { continue; } else {break; }
  }
  
  return 0;
}