core/simulator.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 "simulator.hpp"
#include "four_multi.hpp"
#include <itpp/itbase.h>
#include <itpp/base/itfile.h>
#include <itpp/itsignal.h>
#include <ctime>


using namespace itpp;

void simulate(std::vector<four_multi_node*> *nodes,channel_model *c){
  std::vector<four_multi_node*> cnodes;
  cnodes=*nodes;
  frame_settings next_frame;
  frame_settings next_frame_all[nodes->size()];
  frame_settings next_next_frame_all[nodes->size()];
  std::complex<int16_t> *pr; // Pointers to receive buffers.
  std::complex<int16_t> *pt; // Pointers to transmit buffers.
  double gain_tx, gain_rx, f_tx, f_rx;
  uint32_t no_ant_tx, no_ant_rx;
  uint32_t buffer_length;
  bool all_are_done;
  void* feedback_buffers[(nodes->size())*(nodes->size())];
  bool feedback_buffers_empty[(nodes->size())*(nodes->size())];
  uint32_t num_nodes;

  num_nodes=nodes->size();


  #if 0
  buffer_length=cnodes[0]->d_buffer_length;
  for (uint32_t node_counter=1;node_counter<nodes->size();node_counter++){
    if (!(buffer_length==cnodes[node_counter]->d_buffer_length)) {
      std::cout << "All nodes must have the same buffer length \n";
      exit(1);
    }
  };
  #endif


  for (uint32_t node_counter1=0;node_counter1<nodes->size();node_counter1++){
      for (uint32_t node_counter2=0;node_counter2<nodes->size();node_counter2++){
        cnodes[node_counter1]->d_simulate=true;
        //feedback_buffers_empty[node_counter1+num_nodes*node_counter2]=new bool;
        feedback_buffers_empty[node_counter1+num_nodes*node_counter2]=true;
      };
  };


  for (uint32_t node_counter=0;node_counter<nodes->size();node_counter++){
        cnodes[node_counter]->feedback_buffers=feedback_buffers;
        cnodes[node_counter]->feedback_buffers_empty=feedback_buffers_empty;
        cnodes[node_counter]->num_nodes=num_nodes;
  };

 
  for (uint32_t node_counter=0;node_counter<nodes->size();node_counter++){
     next_frame=cnodes[node_counter]->node_init();     
     while (IS_CALIB_BURST(next_frame)) {
       next_frame=cnodes[node_counter]->node_process();
     };
     next_frame_all[node_counter]=next_frame;
  };


  all_are_done=true;
  for (uint32_t node_counter=0;node_counter<nodes->size();node_counter++){
    next_frame=next_frame_all[node_counter];
    if (!(next_frame.what==END))
      all_are_done=false;
  };
  


  while (!all_are_done) {


    std::cout << "================================================ \n";

    double time;

    uint32_t node_counter0;
    for (node_counter0=0;node_counter0<nodes->size();node_counter0++){
      if (!(next_frame_all[node_counter0].what==END)) {
        time=(double) next_frame_all[node_counter0].time.get_real_secs();
        break;
      }
    };
    for (uint32_t node_counter=node_counter0;node_counter<nodes->size();node_counter++){

      if (!(next_frame_all[node_counter].what==END)) {
        
        if (!(time==(double) next_frame_all[node_counter].time.get_real_secs())) {
          std::cout << "All nodes must have the same frame times \n";
        exit(1);
        };
      };
    };


    std::cout << "time=" << time << "\n";



    for (uint32_t node_rx=0;node_rx<nodes->size();node_rx++){
      no_ant_rx=cnodes[node_rx]->d_no_ant1+cnodes[node_rx]->d_no_ant2;

      
      if (next_frame_all[node_rx].what==RX) {

        for (uint32_t i1=0;i1<no_ant_rx;i1++) {

          if (i1<next_frame_all[node_rx].frequency.size()) {
            f_rx=next_frame_all[node_rx].frequency[i1];
          } else {
            f_rx=next_frame_all[node_rx].frequency.back();
          };
         
          
          if (i1<next_frame_all[node_rx].gain_rx.size()) {
            gain_rx=next_frame_all[node_rx].gain_rx[i1];
          }
          else {
            gain_rx=next_frame_all[node_rx].gain_rx.back();
          };


          pr=(std::complex<short int>*)(cnodes[node_rx]->d_rx_buffers[i1]);
          for (uint32_t i2=0;i2<cnodes[node_rx]->d_buffer_length;i2++) {
            pr[i2]=0;
          };      

          for (uint32_t node_tx=0;node_tx<nodes->size();node_tx++){


            if (next_frame_all[node_tx].what==TX) {

              
              no_ant_tx=cnodes[node_tx]->d_no_ant1+cnodes[node_tx]->d_no_ant2;

              for (uint32_t i2=0;i2<no_ant_tx;i2++) {


                if (i2<next_frame_all[node_tx].frequency.size())
                  f_tx=next_frame_all[node_tx].frequency[i2];
                else
                  f_tx=next_frame_all[node_tx].frequency.back();                
                               
                if (i2<next_frame_all[node_tx].gain_tx.size())
                  gain_tx=next_frame_all[node_tx].gain_tx[i2];
                else
                  gain_tx=next_frame_all[node_tx].gain_tx.back();

                
                pt=(std::complex<short int>*) cnodes[node_tx]->d_tx_buffers[i2];                

                buffer_length=cnodes[node_tx]->d_buffer_length;
                if ((cnodes[node_rx]->d_buffer_length)<buffer_length)
                  buffer_length=cnodes[node_rx]->d_buffer_length;

                c->sub_channel(node_tx,node_rx,i2,i1,f_tx,f_rx,time,gain_tx,
                               gain_rx,pt,pr,buffer_length);
              };   
            };
          };
          // Add noise
          buffer_length=cnodes[node_rx]->d_buffer_length;
          c->add_noise(gain_rx,pr,buffer_length);
        };       
        // Process the received data at the receivers

        if (cnodes[node_rx]->save_input_on_file) {
          cnodes[node_rx]->log_raw_buffer(cnodes[0]->d_buffer_length);
        };     
        next_frame=cnodes[node_rx]->node_process();
        next_next_frame_all[node_rx]=next_frame;
      }; 
      
    };


    // Do all the transmitters
    for (uint32_t node_tx=0;node_tx<nodes->size();node_tx++){
      if (next_frame_all[node_tx].what==TX) {   
        next_frame=cnodes[node_tx]->node_process();
        next_next_frame_all[node_tx]=next_frame;
     };
    };


    // Those who have ended
    for (uint32_t node_tx=0;node_tx<nodes->size();node_tx++){
      if (next_frame_all[node_tx].what==END) {
        next_next_frame_all[node_tx].what=END;
      };
    };
 

    for (uint32_t node_counter=0;node_counter<nodes->size();node_counter++){
      next_frame_all[node_counter]=next_next_frame_all[node_counter];
    };



    all_are_done=true;
    for (uint32_t node_counter=0;node_counter<nodes->size();node_counter++){
      next_frame=next_frame_all[node_counter];
      if (!(next_frame.what==END))
        all_are_done=false;
     };
    

  };
  for (uint32_t node_ix=0;node_ix<nodes->size();node_ix++){
    if (cnodes[node_ix]->save_input_on_file) {
      cnodes[node_ix]->save_logged_raw_data();
    };  
    cnodes[node_ix]->end_of_run();
  };
  std::cout << "Simulation done! \n";

};

void block_fading::sub_channel(uint32_t from_node_ix
,uint32_t to_node_ix, uint32_t from_antenna_ix,  uint32_t to_antenna_ix,
double freq_tx, double freq_rx,  double time, double gain_tx, double gain_rx, 
                  std::complex<short int> const* signal_in, std::complex<short int>* signal_out,uint32_t number_of_samples)
{


  // Save the old seed and create a random
  ivec seed_state;
  RNG_get_state(seed_state);
  if (seed==0)
    RNG_randomize();
  else {
    RNG_reset(seed);
    seed++;
  };
  // Randomize a single coefficient
  std::complex<double> h;
  itpp::cvec s_in, s_out;
  
  // Set length of ITPP vectors
  s_in.set_length(number_of_samples);
  s_out.set_length(number_of_samples);

  // Convert input to cvec
  for (uint32_t i1=0;i1<number_of_samples;i1++) {
    s_in[i1]=signal_in[i1];
  };

  // Convert receiver buffer to cvec
  for (uint32_t i1=0;i1<number_of_samples;i1++) {
    s_out[i1]=signal_out[i1];
  };


  // Do the channel model
  mat temp;
  double t;
  temp=old_time[from_node_ix][to_node_ix];
  t=temp.get(from_antenna_ix,to_antenna_ix);

  for (uint32_t i3=0;i3<num_of_taps;i3++) {

    if (time>t) {
      h=itpp::randn_c();
      h=forgetting_factor*(old_h[from_node_ix][to_node_ix][i3]).
        get(from_antenna_ix,to_antenna_ix)+
      sqrt(1-forgetting_factor*forgetting_factor)*h;
     old_h[from_node_ix][to_node_ix][i3].set(from_antenna_ix,to_antenna_ix,h);
    } else {
      h=(old_h[from_node_ix][to_node_ix][i3]).get(from_antenna_ix,                                                to_antenna_ix);
    };

    if (is_awgn) {
      if (i3==0)
        h=1;
      else
        h=0;
    };
    if (is_diagonal) {
      if (from_antenna_ix!=to_antenna_ix)  {
        h=0; 
      };
    };
 
    h=path_gain[from_node_ix][to_node_ix]*scaling_signal*h;
    h=h/sqrt(num_of_taps);
    for (int i4=i3;i4<s_in.length();i4++) {
      s_out(i4)+=h*s_in(i4-i3);
    };
  };
  // Convert output to  std::complex<short int>
  for (uint32_t i1=0;i1<number_of_samples;i1++) {
    signal_out[i1]=s_out[i1];
  };
  
  // Reset the seed
  RNG_set_state(seed_state);

};

void block_fading::add_noise(double gain_rx, 
        std::complex<short int>* signal, 
                                 uint32_t number_of_samples) {

  itpp::cvec s, n;

  // Save the old seed and create a random
  ivec seed_state;
  RNG_get_state(seed_state);
  RNG_randomize();


  // Set length of ITPP vectors
  s.set_length(number_of_samples);

  // Convert input to cvec
  for (uint32_t i1=0;i1<number_of_samples;i1++) {
    s[i1]=signal[i1];
  };

  // Generate noise
  n=itpp::randn_c(number_of_samples);
  n=scaling_noise*n;


  // Add noise
  s=s+n;

  // Convert output to  std::complex<short int>
  for (uint32_t i1=0;i1<number_of_samples;i1++) {
     signal[i1]=s[i1];
  };

  // Reset the seed
  RNG_set_state(seed_state);

  

}

block_fading::block_fading(void) {

  scaling_signal=0.45;
  scaling_noise=4.5;//21.7;
  seed=0;
  forgetting_factor=0.0;
  is_awgn=false;
  is_diagonal=false;
  num_of_taps=1;
  // Save the old seed and create a random
  ivec seed_state;
  RNG_get_state(seed_state);
  if (seed==0)
    RNG_randomize();
  else {
    RNG_reset(seed);
    seed=seed+1;
  };


  for (uint32_t i1=0;i1<max_num_nodes;i1++) {
    for (uint32_t i2=0;i2<max_num_nodes;i2++) {
      old_time[i1][i2].set_size(max_num_ant,max_num_ant);
      old_time[i1][i2].zeros();
      path_gain[i1][i2]=1.0;
      for (uint32_t i3=0;i3<max_num_taps;i3++) {
        old_h[i1][i2][i3].set_size(max_num_ant,max_num_ant);
        old_h[i1][i2][i3]=randn_c(max_num_ant,max_num_ant);
      };
    };
  };
  // Reset the seed
  RNG_set_state(seed_state);


};

void block_fading::reset(void) {

  ivec seed_state;

  RNG_get_state(seed_state);
  if (seed==0)
    RNG_randomize();
  else {
    RNG_reset(seed);
    seed++;
  };

  for (uint32_t i1=0;i1<max_num_nodes;i1++) {
    for (uint32_t i2=0;i2<max_num_nodes;i2++) {
      old_time[i1][i2].zeros();
      for (uint32_t i3=0;i3<max_num_taps;i3++) {        
        old_h[i1][i2][i3]=randn_c(max_num_ant,max_num_ant);
      };
    };
  };

  //std::cout << "old_h=" << old_h[0][0][0](0,0) << "\n";

  // Reset the seed
  RNG_set_state(seed_state);


};

from_file::from_file(std::string directory_name) {
  raw_data_dir=directory_name;
  for (uint32_t i1=0;i1<max_num_nodes;i1++) {
    for (uint32_t i2=0;i2<max_num_nodes;i2++) {
      old_time[i1][i2].set_size(max_num_ant,max_num_ant);
      old_time[i1][i2]=-ones(max_num_ant,max_num_ant);
      frame_ix[i1][i2].set_size(max_num_ant,max_num_ant);
      frame_ix[i1][i2]=-ones_i(max_num_ant,max_num_ant);

    };
  };
};







void from_file::sub_channel(uint32_t from_node_ix,uint32_t to_node_ix, 
                         uint32_t from_antenna_ix, uint32_t to_antenna_ix 
                            , double freq_tx, double freq_rx,
                                 double time, double gain_tx, double gain_rx, 
                            std::complex<short int> const* signal_in, std::complex<short int>* signal_out,uint32_t number_of_samples) {


  mat temp;
  double t;
  int f_ix;
  std::string filename;
  std::stringstream temp_ss;
  std::ifstream s1; 

  temp=old_time[from_node_ix][to_node_ix];
  t=temp.get(from_antenna_ix,to_antenna_ix);

  f_ix=frame_ix[from_node_ix][to_node_ix](from_antenna_ix,to_antenna_ix);

  if (time>t) {
    f_ix++;
    (frame_ix[from_node_ix][to_node_ix])(from_antenna_ix,to_antenna_ix)=f_ix;
    (old_time[from_node_ix][to_node_ix])(from_antenna_ix,to_antenna_ix)=time;
  };

  temp_ss << to_node_ix;
  filename=raw_data_dir+"node="+temp_ss.str()+"_";
  temp_ss.str("");
  temp_ss << f_ix;
  filename+="frame="+temp_ss.str()+".dat";
  temp_ss.str("");
  //std::cout << "filename=" << filename << std::endl;

  s1.open(filename.c_str(), std::ios::binary);
  for (uint32_t i1=0;i1<to_antenna_ix+1;i1++) {
    s1.read((char *) signal_out,number_of_samples*4); 
  };
  s1.close(); 

  
};

OFDM_channel_from_file::OFDM_channel_from_file(std::string filename,
ivec num_ant,int Nfft,int Nc, 
             ivec ix_all ,int tx_start_sample, int rx_start_sample,
            double signal_scaling, double EVM,double EVM_RX, double CPE) {

  d_EVM=EVM;
  d_EVM_RX=EVM_RX;
  d_CPE=CPE;
  d_filename=filename;
  d_num_ant=num_ant;
  d_num_nodes=d_num_ant.length();
  d_Nfft=Nfft;
  d_signal_scaling=signal_scaling;
  d_Nc=Nc;
  d_tx_start_sample=tx_start_sample;
  d_rx_start_sample=rx_start_sample;
  d_frame_ix=0;
  t_old=0;
  d_ix_all=ix_all;
  counter=0;
  rand_cpe.set_size(max_num_ant,max_num_nodes);

  for (int i1=0;i1<d_num_nodes;i1++) {
    for (int i2=0;i2<d_num_nodes;i2++) {
      for (int i3=0;i3<d_num_ant(i1);i3++) {
        (H[i1][i2][i3]).set_size(d_num_ant(i2),ix_all.length());
      };
    };
  };


  std::stringstream temp_ss;
  temp_ss << d_frame_ix;  
  load_file(d_filename+"_frame_ix="+temp_ss.str()+".dat");

};

void OFDM_channel_from_file::reset(void) {
  d_frame_ix=0;
  t_old=0;
  counter=0;
  std::stringstream temp_ss;
  temp_ss << d_frame_ix;  
  load_file(d_filename+"_frame_ix="+temp_ss.str()+".dat");
};

void OFDM_channel_from_file::load_file(std::string filename) {


 std::ifstream s1;   
 s1.open(filename.c_str(), std::ios::binary);
 double r,i;
 
 for (int i2=0;i2<d_num_nodes;i2++) {
   for (int i1=0;i1<d_num_nodes;i1++) {
     

    for (int i5=0;i5<d_ix_all.length();i5++) { 
      for (int i4=0;i4<d_num_ant(i1);i4++) {
         for (int i3=0;i3<d_num_ant(i2);i3++) {  
           

           s1.read((char *) &r,sizeof(double));
           real(H[i1][i2][i4](i3,i5))=r;

           
           
         };
       };
     };


    for (int i5=0;i5<d_ix_all.length();i5++) { 
      for (int i4=0;i4<d_num_ant(i1);i4++) {
         for (int i3=0;i3<d_num_ant(i2);i3++) {  
           s1.read((char *) &i,sizeof(double));
           imag(H[i1][i2][i4](i3,i5))=i;

         };
       };
     };

   };
 };
 s1.close();

 double x,y;
 vec tempn;
 //rand_cpe=d_CPE*randn(max_num_ant,max_num_nodes);
 for (int i1=0;i1<(int) max_num_ant;i1++) {
   for (int i2=0;i2<(int) max_num_nodes;i2++) {
     tempn=randu(2);
     x=tempn(0);
     if (x<0.5)
       y=(1/sqrt(2))*log(2.0*x);
     else
       y=(1/sqrt(2))*log(2.0*(1-x));

     if (x>0.5){
       y=-y;
     };     
     rand_cpe(i1,i2)=d_CPE*y;
   };
 };


};


void OFDM_channel_from_file::sub_channel(uint32_t from_node_ix,
uint32_t to_node_ix, 
           uint32_t from_antenna_ix, uint32_t to_antenna_ix 
                 , double freq_tx, double freq_rx,
                                 double time, double gain_tx, double gain_rx, 
std::complex<short int> const* signal_in, std::complex<short int>* signal_out,
uint32_t number_of_samples) {

  cvec s_in, s_out;
  
  if (time>t_old) {
    d_frame_ix++;
    std::stringstream temp_ss;
    temp_ss << d_frame_ix;  
    load_file(d_filename+"_frame_ix="+temp_ss.str()+".dat");
    t_old=time;
  };


  // Set length of ITPP vectors
  s_in.set_length(number_of_samples);
  s_out.set_length(number_of_samples);

  // Convert input to cvec
  for (uint32_t i1=0;i1<number_of_samples;i1++) {
    s_in[i1]=signal_in[i1];
  };
  s_in/=d_signal_scaling;

  // Convert receiver buffer to cvec
  for (uint32_t i1=0;i1<number_of_samples;i1++) {
    s_out[i1]=signal_out[i1];
  };
  
  // Do the channel model
  int sample_ix=d_tx_start_sample;
  cvec w(d_Nfft), W(d_Nfft);
  cvec Hc(d_Nfft);
  cvec st(d_Nfft), St(d_Nfft);
  bool cond=true;

  Hc.zeros();

  while (cond) {

    for (int i1=0;i1<d_ix_all.length();i1++){
      Hc(d_ix_all(i1))
        =H[from_node_ix][to_node_ix][from_antenna_ix](to_antenna_ix,i1);
    };
    Hc=std::polar(1.0,rand_cpe(from_antenna_ix,from_node_ix)*0.017453)*Hc;

    w=s_in.mid(sample_ix,d_Nfft);
    W=fft(w);    
    St=elem_mult(Hc,W);
    st=ifft(St);
    double c=sum_sqr(abs(st))/((double) d_Nfft);
    
    c=c*(((double) (d_Nfft))/((double) (d_ix_all.length())))*
      (((double) (d_Nfft))/((double) (d_ix_all.length())));
    c=d_EVM*sqrt(c);
    st=st+c*randn_c(d_Nfft);
    
    
    for (int i1=0;i1<d_Nfft;i1++) {
      s_out[sample_ix+d_rx_start_sample-d_tx_start_sample+i1]+=st(i1);
    };
    sample_ix+=d_Nfft+d_Nc;

    cond=sample_ix<((int) number_of_samples-d_Nfft);
    cond=cond && (sample_ix+d_rx_start_sample-d_tx_start_sample+d_Nfft+d_Nc <
                  (int) number_of_samples);

      counter++;

  };


  // Convert output to  std::complex<short int>
  for (int i1=0;i1<(int) number_of_samples;i1++) {
    signal_out[i1]=s_out[i1];
  };   
};


void OFDM_channel_from_file::add_noise(double gain_rx, 
        std::complex<short int>* signal, 
                                 uint32_t number_of_samples) {

  itpp::cvec s, n, n_evm;

  // Save the old seed and create a random
  ivec seed_state;
  RNG_get_state(seed_state);
  RNG_randomize();


  // Set length of ITPP vectors
  s.set_length(number_of_samples);

  // Convert input to cvec
  for (uint32_t i1=0;i1<number_of_samples;i1++) {
    s[i1]=signal[i1];
  };

  // Calculate power for receive EVM generation
  double power=0.0;
  double power_local=0;
  int sample_ix=d_rx_start_sample;
  int no_of_symbols=0;
  bool cond=sample_ix<((int) number_of_samples-d_Nfft);
  while (cond) {
    power_local=0.0;
    for (int i1=0;i1<d_Nfft;i1++) {
      power_local+=abs(s(sample_ix))*abs(s(sample_ix));
    };
    power_local=power_local/((double) d_Nfft);
    if (sqrt(power_local)>10) {
      power+=power_local;
      no_of_symbols++;
    };
    sample_ix+=d_Nfft+d_Nc;
    cond=sample_ix<((int) number_of_samples-d_Nfft);
  };

  power*=(((double) (d_Nfft))/((double) (d_ix_all.length())))*
      (((double) (d_Nfft))/((double) (d_ix_all.length())));
  power=power/((double) no_of_symbols);

  // Generate noise
  n=itpp::randn_c(number_of_samples);
  n=scaling_noise*n;

  // Generate RX EVM
  n_evm=itpp::randn_c(number_of_samples);
  n_evm=d_EVM_RX*sqrt(power)*n_evm;


  // Add noise and EVM
  //s=s+n+n_evm;

  // Convert output to  std::complex<short int>
  for (uint32_t i1=0;i1<number_of_samples;i1++) {
     signal[i1]=s[i1];
  };

  // Reset the seed
  RNG_set_state(seed_state);

  

}