modem/AMC.cpp

// Copyright 2011-2013, Per Zetterberg, KTH Royal Institute of Technology
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program.  If not, see <http://www.gnu.org/licenses/>.
//


#include "LDPC1.hpp"
#include "AMC.hpp"
#include <itpp/itcomm.h>
#include <iostream>
#include <fstream>
#include <math.h>
#include <cmath>
#include <string>
#include <unistd.h>

using namespace itpp;
using namespace std;

AMC::AMC(void) {};

void AMC::init(uint32_t code_word_size,bool load_from_file) {

  b0.init(0.25,code_word_size,10,load_from_file,"code_data1.it");
  b1.init(0.5,code_word_size,10,load_from_file,"code_data2.it");
  b2.init(0.625,code_word_size,10,load_from_file,"code_data3.it");
  b3.init(0.75,code_word_size,10,load_from_file,"code_data4.it");

  m_code_word_size=code_word_size;

  soft_bit.set_length(code_word_size);
  only_raw_bits=false; 
};

int AMC::codec_message_size(uint32_t codec_ix)  {
  return codec(codec_ix)->message_length();
}

void AMC::encode(uint32_t codec_ix, bvec message, bvec &transmitted) {
  codec(codec_ix)->chan_encode(message,transmitted);
};

void AMC::modulate(uint32_t modulation_order, uint32_t codec_ix,  bvec message, bvec &transmitted, cvec &symbols,int skip_symbols)  {
  ivec ti(1);    
  cvec ci(1);

  qam.set_M(modulation_order);
  codec(codec_ix)->chan_encode(message,transmitted);  

  int f=qam.bits_per_symbol();
  ivec b2s=qam.get_bits2symbols();
  int symbol_index;

  for (int i1=0;i1<transmitted.length()/f;i1++) {

    symbol_index=b2s(bin2dec(transmitted(i1*f,i1*f+f-1)));    
    ti(0)=symbol_index;
    qam.modulate(ti,ci);
    symbols(i1+skip_symbols)=ci(0);
        
  };


};


double AMC::evm(uint32_t modulation_order, bvec &transmitted,  cvec &symbols, int skip_symbols) {
  ivec ti(1);    
  cvec ci(1);
  double evm=0.0;
  ivec b2s=qam.get_bits2symbols();
  int symbol_index;

  qam.set_M(modulation_order);
  int f=qam.bits_per_symbol();

  for (int i1=0;i1<transmitted.length()/f;i1++) {
    symbol_index=b2s(bin2dec(transmitted(i1*f,i1*f+f-1)));
    ti(0)=symbol_index;
    qam.modulate(ti,ci);
    evm+=sqr(symbols(i1+skip_symbols)-ci(0));
    
  };
  evm/=(transmitted.length()/f);
  evm=sqrt(evm);
  return evm;

};



void AMC::demodulate(const cvec &input_symbols, double noise_variance, uint32_t modulation_order, uint32_t codec_ix, bvec &message_hat, bvec &transmitted_hat) 
 {

  qam.set_M(modulation_order);
  for (int i1=1;i1<=log2(modulation_order);i1++) {
    soft_bit(m_code_word_size-1)=0;
  };

  qam.demodulate_soft_bits(input_symbols,noise_variance,soft_bit);

  transmitted_hat=soft_bit<0;

  if (!only_raw_bits) {
    codec(codec_ix)->chan_decode(soft_bit,message_hat);
  };

};

double AMC::noise_variance(const cvec &symbols, uint32_t modulation_order){

  QAM modulation;
  double noise_variance_norm;
  cvec closest_constellation_point(1);
  cvec symbol(1);

  modulation.set_M(modulation_order);  

  // Estimate noise variance in symbols
  noise_variance_norm=0;
  for (int i1=0;i1<symbols.length();i1++) {
    symbol=symbols(i1);
    closest_constellation_point=modulation.modulate_bits(
                      modulation.demodulate_bits(symbol));

       
    noise_variance_norm+=sum(sqr(symbol-closest_constellation_point));
  };
  noise_variance_norm/=symbols.length();
  if (noise_variance_norm<1e-3)
    noise_variance_norm=1e-3;

  return noise_variance_norm;
}


uint32_t AMC::rand_new(uint64_t *next)  {
  uint32_t temp;
  (*next) = (*next) * 1103515245 + 12345;
  temp = (*next) >> 16;
  return (temp % 32768);
};

void AMC::randb_thread_safe(int size, bvec &out, uint64_t *seed)  {
  randb_thread_safe(0,size,out,seed);
};


void AMC::randb_thread_safe(int start_ix, int length, bvec &out, uint64_t *seed)  {
  for (int i1=start_ix;i1<start_ix+length;i1++) {
    out(i1)=rand_new(seed) & 1;
  }; 
};


void AMC::run_through_all_mcs(void)  {
   #define NO_OF_FRAMES 1000
  
   ivec modulation_orders="4,16,64,256";

   vec sigmas="1.5849    1.4125    1.2589    1.1220    1.0000    0.8913  \
               0.7943    0.7079    0.6310    0.5623    0.5012    0.4467  \
               0.3981    0.3548    0.3162    0.2818    0.2512    0.2239  \
               0.1995    0.1778    0.1585    0.1413    0.1259    0.1122  \
               0.1000    0.0891    0.0794    0.0708    0.0631    0.0562  \
               0.0501    0.0447    0.0398    0.0355    0.0316    0.0282 ";


   int modulation_order;
   int l,r;
   double FER, BER, Rate, BERc, BERc_raw, BER_raw;
   mat results(modulation_orders.length()*no_of_codecs(),sigmas.length());
   mat BERs(modulation_orders.length()*no_of_codecs(),sigmas.length());
   mat BERs_raw(modulation_orders.length()*no_of_codecs(),sigmas.length());

   vec Rates(modulation_orders.length()*no_of_codecs());
   bvec input, transmitted, message_hat, transmitted_hat;
   cvec symbols;

   Real_Timer timer;
   double max_time=0;
   BERC berc;
   double s2;
  
   input.set_size(codec_message_size(no_of_codecs()-1));
   message_hat.set_size(codec_message_size(no_of_codecs()-1));
   transmitted.set_size(m_code_word_size);
   transmitted_hat.set_size(m_code_word_size);
   symbols.set_size(m_code_word_size/log2((modulation_orders[0])));
   RNG_reset(10);               
   cout << "Hang on ... this is gonna take a while ...\n";

   for (int i1=0;i1<modulation_orders.length();i1++) {
     modulation_order=modulation_orders[i1];
     cout << "Modulation order=" << modulation_order << endl;
     usleep(1000);

     for (int c=0;c<no_of_codecs();c++) {
       cout << "Code number=" << c << endl;       
       l=codec_message_size(c);

       r=i1+c*modulation_orders.length();
       Rate=(float (l))/((float) m_code_word_size)*
                             log2((float) modulation_order);
       Rates(r)=Rate;
       max_time=0;
  
       for (int i4=0;i4<sigmas.length();i4++) {
         FER=0;
         BERc=0;
         BERc_raw=0;
         usleep(10);


         for (uint32_t f=0;f<NO_OF_FRAMES;f++) {
           
           input.set_subvector(0,randb(l));
           modulate(modulation_order,c,input.mid(0,l),transmitted,symbols,0);
           usleep(10);     

           symbols=symbols+sigmas(i4)*randn_c(symbols.length());
           timer.reset();
           timer.start();

           s2=sigmas(i4)*sigmas(i4);
           s2=noise_variance(symbols.mid(0,m_code_word_size/log2( modulation_order)), modulation_order);

           demodulate(symbols.mid(0,
             m_code_word_size/log2( modulation_order)),
                      s2,
                          modulation_order,c,message_hat,transmitted_hat); 
           
           timer.stop();
           if (timer.get_time()>max_time)
             max_time=timer.get_time();

           berc.clear();
           berc.count(transmitted,transmitted_hat);  
           BER_raw=berc.get_errorrate();

           berc.clear();
           berc.count(input.mid(0,l),message_hat.mid(0,l));
           BER=berc.get_errorrate();
           if (BER>0)
             FER=FER+(1.0/((double) NO_OF_FRAMES));

           BERc=BERc+(BER/((double) NO_OF_FRAMES));
           BERc_raw=BERc_raw+(BER_raw/((double) NO_OF_FRAMES));
           

         };

         std::cout << "res=" << i4 <<"," << s2 << "," << "FER=" << FER << std::endl;
         results(r,i4)=FER;
         BERs(r,i4)=BERc;
         BERs_raw(r,i4)=BERc_raw;
       };
       std::cout << "Max used time=" << max_time << endl;
     };
   }; 


  it_file d1;
  d1.open("results_MCS.it");
  d1<< Name("results") << results;
  d1<< Name("sigmas") << sigmas;
  d1<< Name("Rates") << Rates;
  d1<< Name("BERs") << BERs;
  d1<< Name("BERs_raw") << BERs_raw;
  d1.close();

  // Plot results in matlab
  // plot(-10*log10(sigmas),(1-results).*repmat(Rates,1,23),'-x') %Waterfall
  // plot(-10*log10(sigmas),results,'-x') % FER.


};