#include "stdafx.h"
#include "TestFunction.h"


void Test_STBC(void)
{

  MIMO  mimo;
  QAM qpsk(16);
  CSTBC STCoding;     //Space-Time Coding
  MIMOChannel channel;
  const vec EbN0dB = "-5:1:20";        // SNR range
  int Number_of_bits;
  double Ec, Eb;
  vec  EbN0, N0, noise_variance, bit_error_rate; //vec is a vector containing double


  Ec = 1.0;                      //The transmitted energy per QPSK symbol is 1.
  Eb = Ec / 2.0;                 //The transmitted energy per bit is 0.5.
  //EbN0dB = linspace(0.0,9.0,10); //Simulate for 10 Eb/N0 values from 0 to 9 dB.
  EbN0 = inv_dB(EbN0dB);         //Calculate Eb/N0 in a linear scale instead of dB.
  N0 = Eb * pow(EbN0,-1.0);      //N0 is the variance of the (complex valued) noise.
  Number_of_bits = 1000000;       //One hundred thousand bits is transmitted for each Eb/N0 value

  Array<BERC> beruc(length(EbN0dB));  // counter for uncoded BER
  bit_error_rate.set_size(EbN0dB.length(),false);

  for (int nsnr=0; nsnr<length(EbN0dB); nsnr++)
  {
    bvec transmitted_bits,received_bits;
    cvec transmitted_signal,received_signal;
    cmat tr_STSignal, rec_STSignal;
    //general the transmitted bits 
    transmitted_bits = randb(Number_of_bits);
    
    //modulate the transmitted bits
    transmitted_signal = qpsk.modulate_bits(transmitted_bits);
    
    //ST encoder
    STCoding.Encoder(transmitted_signal,tr_STSignal);
    
    //through the channel
    channel.IniChannel(2,2);
    channel.set_noise(N0(nsnr));
    rec_STSignal = channel(tr_STSignal);
    //cmat rec_STSignal = tr_STSignal;
    
    //ST decoder
    STCoding.Decoder(rec_STSignal,channel.GetChannel(),received_signal);
    //STCoding.Decoder(tr_STSignal,channel.GetChannel(),received_signal);  
    
    //Demodulate the bits
    received_bits = qpsk.demodulate_bits(received_signal/sqrt(2.0)); 
	//Count the number of errors
    beruc(nsnr).count(transmitted_bits,received_bits);
    //bit_error_rate(nsnr) = beruc(nsnr).get_errorrate();

	cout<<endl<<"the SNR:"<<EbN0dB(nsnr)<<"has been simulated"<<endl;
	cout<<"bit_error_rate"<<beruc(nsnr).get_errorrate();
  }

  for (int nsnr=0; nsnr<length(EbN0dB); nsnr++)
  {
    bit_error_rate(nsnr) = beruc(nsnr).get_errorrate();

  }

  //////////////////////////////////////////////////
  //Save simulation result to File.
  //////////////////////////////////////////////////
  it_file ff("Test_STBC1.it");

  // Save channel coefficients to the output file
  //ff << Name("transmitted_bits") << transmitted_bits;
  //ff << Name("transmitted_signal") << transmitted_signal;

  //ff << Name("channel") << channel.GetChannel();
  //ff << Name("tr_STSignal") << tr_STSignal;
  //ff << Name("rec_STSignal") << rec_STSignal;


  //ff << Name("received_signal") << received_signal;
  //ff << Name("received_bits") << received_bits;
  ff << Name("bit_error_rate")<<bit_error_rate;
  // Close the output file
  ff.close();
}

void Test_AWGNChannel()
{
  //Initiate the AWGN_Channel class
  double noisevar = 0.001;
  AWGN_Channel awgn_channel(noisevar);
  
  //Initiate a QPSK-modulator, and generate the transmitted signal
  QAM qpsk(16);
  bvec transmitted_bits = randb(500);
  cvec transmitted_signal = qpsk.modulate_bits(transmitted_bits);
  
  //Usage of the member operator ()
  cvec received_signal = awgn_channel(transmitted_signal);
  
  //Demodulate the bits
  bvec received_bits = qpsk.demodulate_bits(received_signal);

  //////////////////////////////////////////////////
  //Save simulation result to File.
  //////////////////////////////////////////////////
  it_file ff("Test_AWGNChannel.it");

  // Save channel coefficients to the output file
  ff << Name("transmitted_bits") << transmitted_bits;
  ff << Name("transmitted_signal") << transmitted_signal;
  ff << Name("received_signal") << received_signal;
  ff << Name("received_bits") << received_bits;
  // Close the output file
  ff.close();

}


int Test_MIMO(int nC,int nRx,int nTx,int Tc)
{
  // -- 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)

  ////add by ssw
  //nC = 2;
  //nRx = 2;
  //nTx = 2;
  //Tc = 1;
 ////add by ssw


 /* 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 int Nbitsmax = 50000000;  // 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_M(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 = inv_dB(-EbN0db(nsnr));
    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)) 
		{
          chan.demodulate_soft_bits(Y(k), H(k/Tc), sigma2, llr_apr, llr_apost,
                                    ND_UQAM::ZF_LOGMAP);
          LLRin(0).set_subvector(k*Nbitspvec,(k+1)*Nbitspvec-1,llr_apost);
        }

        // ML demodulation
        if (Contflag(1)) 
		{
          chan.demodulate_soft_bits(Y(k), H(k/Tc), sigma2, llr_apr, llr_apost);
          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;
}