modem/modem_OFDM1.cpp

/*
 * Copyright 2012-2013, Per Zetterberg, KTH Royal Institute of Technology
 *
 * This program is free software: you can redistribut 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/>.
 */

/* This code should implemen, more or less, the same functions as 
 * 
 * mod_OFDM1.m and demod_OFDM1.m
 *
 *
 */

#include <itpp/itcomm.h>
#include <itpp/stat/misc_stat.h>
#include <math.h>
#include "modem_OFDM1.hpp"
#include <itpp/itbase.h>
#include <complex>
#include <itpp/base/sort.h>


using namespace itpp;
using namespace std;

OFDM1::OFDM1(void){
  pilot_power_factor=1;
  estimate_SINR=false;
  pulse_shape="0.0546,-0.0898,-0.0587,0.3134,0.5612,0.3133,-0.0587,-0.0898,0.0546";
  pulse_shape_delay=4;
};


uint32_t OFDM1::number_of_modulation_symbols(void){
  return  ((Ns-known_symbol_pos.length())*length_ix);
};

uint32_t OFDM1::waveform_length(void){
  return (max(d_reorder)*Ns*(Nfft+Np));
};


void OFDM1::init(bool use_pilot_carriers, bool prepend_training,uint32_t prefix_length,uint32_t number_of_symbols, ivec known_symbol_postion_index, ivec reorder){
    
    pz_init_fft();
    Ns=number_of_symbols;
    known_symbol_pos=known_symbol_postion_index;
    d_reorder=reorder;

    do_use_pilot_carriers=use_pilot_carriers;
    do_prepend_training=prepend_training;

    Np=prefix_length;

    if (do_use_pilot_carriers) 
      length_ix=36;
    else
      length_ix=38;

    length_ix_all=38;

    pilot_symbol=1.0+1.0i;
    h.set_size(Nfft);
    
    
    string temp_str;
    temp_str="1,1,-1,1,1,1,1,-1,1,1,1,-1,1,1,1,-1,-1,-1,1,1,1,1, ";
    temp_str+="-1,-1,-1,-1,-1,1,-1,-1,1,-1,1,1,-1,1,1,1,-1,1,1,1,1,-1,-1,-1, ";
    temp_str+="-1,1,1,1,1,1,1,-1,-1,-1,1,1,-1,1,-1,-1,-1,-1,-1,-1,1,-1,-1,1, ";
    temp_str+="1,1,1,-1,1,-1,-1,-1,-1,1,1,-1,1,-1,-1,-1,1,1,-1,1,-1,-1,-1, ";
    temp_str+="-1,1,1,-1,-1,1,-1,1,-1,1,1,1,1,1,-1,1,1,1,-1,1,-1,-1,1,-1,-1, ";
    temp_str+="1,1,-1,1,1,1,1,1,1,-1,1,-1,-1,-1,1,-1,-1,1,-1,1 ";
    train=temp_str.c_str();
    ix.set_size(length_ix);



    { uint32_t counter=0;
      for (uint32_t i1=1;i1<pilot_carrier1;i1++) {
        ix(counter++)=i1;
      };

      for (uint32_t i1=pilot_carrier1+1;i1<=rightmost_carrier;i1++) {
        ix(counter++)=i1;
      };
      rightmost_carrier_sub=counter-1;
      leftmost_carrier_sub=counter;

      for (uint32_t i1=leftmost_carrier;i1<pilot_carrier2;i1++) {
        ix(counter++)=i1;
      };
      
      for (uint32_t i1=pilot_carrier2+1;i1<Nfft;i1++){
        ix(counter++)=i1;
      };
    };
    
    ix_all.set_size(length_ix_all);
    {uint32_t counter=0;

      for (uint32_t i1=1;i1<=rightmost_carrier;i1++) {
        ix_all(counter++)=i1;
      };
      for (uint32_t i1=leftmost_carrier;i1<Nfft;i1++) {
        ix_all(counter++)=i1;
      };

    };


    ix_sub.set_size(length_ix);
    {int counter=0;
      for (uint32_t i1=0;i1<(length_ix+2);i1++) {
        for (uint32_t i2=0;i2<length_ix;i2++) {
          if (ix_all(i1)==ix(i2)) {
            ix_sub(counter)=i1;
            counter++;
          };
        };
      };
    };
  
    for (uint32_t i1=0;i1<length_ix+2;i1++) {
      if (ix_all(i1)==(int32_t) pilot_carrier1) {
        pilot_carrier1_sub=i1;
      };

      if (ix_all(i1)==(int32_t) pilot_carrier2) {
        pilot_carrier2_sub=i1;
      };
    };

    if (!do_use_pilot_carriers) {
      ix=ix_all;
      for (uint32_t i1=0;i1<length_ix;i1++)
        ix_sub(i1)=i1;
    };

    int no_interleavers=8;
    symbol_interleave_matrix.set_size(no_interleavers,length_ix);
    symbol_de_interleave_matrix.set_size(no_interleavers,length_ix);
    RNG_reset(101);
    ivec temp1(length_ix), temp2(length_ix);
    Sort<double> my_sort_f;
    Sort<int> my_sort_i;

    for (int i1=0;i1<no_interleavers;i1++) {
      temp1=my_sort_f.sort_index(0,(int) (length_ix-1),randu(length_ix));
      temp2=my_sort_i.sort_index(0,(int) (length_ix-1),temp1);
      
      symbol_interleave_matrix.set_row(i1,temp1);
      symbol_de_interleave_matrix.set_row(i1,temp2);      
    };


      

  } ;

void OFDM1::pilot_pattern(cvec &burst, cvec precoder) {


  burst(0)=0.0+0.0i;  burst(1)=-1.0-1.0i; burst(2)=-1.0-1.0i; 
  burst(3)=1.0+1.0i; burst(4)=-1.0-1.0i; 
  burst(5)=1.0+1.0i; burst(6)=-1.0-1.0i; burst(7)=1.0+1.0i; 
  burst(8)=1.0+1.0i; burst(9)=-1.0-1.0i; 
  burst(10)=-1.0-1.0i; burst(11)=-1.0+1.0i; 
  burst(12)=-1.0+1.0i; burst(13)=1.0+1.0i; 
  burst(14)=-1.0+1.0i;  burst(15)=-1.0+1.0i; 
  burst(16)=-1.0+1.0i; burst(17)=1.0+1.0i; 
  burst(18)=1.0+1.0i; burst(19)=1.0+1.0i; 
  burst(20)=0.0+0.0i; burst(21)=0.0+0.0i; 
  burst(22)=0.0+0.0i; burst(23)=0.0+0.0i; 
  burst(24)=0.0+0.0i; burst(25)=0.0+0.0i; 
  burst(26)=0.0+0.0i; burst(27)=0.0+0.0i; 
  burst(28)=0.0+0.0i; burst(29)=0.0+0.0i; 
  burst(30)=0.0+0.0i; burst(31)=0.0+0.0i; 
  burst(32)=0.0+0.0i; burst(33)=0.0+0.0i; 
  burst(34)=0.0+0.0i; burst(35)=0.0+0.0i; 
  burst(36)=0.0+0.0i; burst(37)=0.0+0.0i; 
  burst(38)=0.0+0.0i; burst(39)=0.0+0.0i; 
  burst(40)=0.0+0.0i; burst(41)=0.0+0.0i; 
  burst(42)=0.0+0.0i; burst(43)=0.0+0.0i; 
  burst(44)=0.0+0.0i; 
  burst(45)=0.0+0.0i; burst(46)=0.0+0.0i; 
  burst(47)=0.0+0.0i; burst(48)=0.0+0.0i; 
  burst(49)=0.0+0.0i; burst(50)=0.0+0.0i; 
  burst(51)=0.0+0.0i; burst(52)=0.0+0.0i; 
  burst(53)=0.0+0.0i; burst(54)=0.0+0.0i; 
  burst(55)=0.0+0.0i; burst(56)=0.0+0.0i; 
  burst(57)=0.0+0.0i; burst(58)=0.0+0.0i; 
  burst(59)=0.0+0.0i; burst(60)=0.0+0.0i; 
  burst(61)=-1.0+1.0i; burst(62)=1.0+1.0i; 
  burst(63)=-1.0-1.0i; burst(64)=-1.0+1.0i; 
  burst(65)=1.0+1.0i; burst(66)=-1.0-1.0i; 
  burst(67)=-1.0-1.0i; burst(68)=1.0+1.0i; 
  burst(69)=-1.0+1.0i;  burst(70)=-1.0-1.0i; burst(71)=-1.0+1.0i; 
  burst(72)=-1.0-1.0i; burst(73)=-1.0+1.0i; burst(74)=-1.0+1.0i; 
  burst(75)=-1.0-1.0i; burst(76)=1.0+1.0i; 
  burst(77)=-1.0-1.0i; burst(78)=-1.0+1.0i; 
  burst(79)=1.0+1.0i;

  for (uint32_t i1=0;i1<Nfft;i1++) {
    burst(i1)=burst(i1)*pilot_power_factor/sqrt(2);
  };


  for (uint32_t i1=0;i1<length_ix_all;i1++) {
    burst(ix_all(i1))=burst(ix_all(i1))*precoder(i1);
  };


}; 


ivec OFDM1::start_of_known_symbol(void){
  ivec temp;
  temp.set_length(known_symbol_pos.length());
  for (int i1=0;i1<known_symbol_pos.length();i1++) {
    temp(i1)=d_reorder.get(known_symbol_pos(i1))*(Nfft+Np);
  };
  return temp;
};



void OFDM1::modulate(cvec &waveform,const cvec &symbols_in,  
                     const uint32_t offset,
double &rms_power,const cvec &precoder)

{

    cvec symbol(Nfft), symbol_temp(Nfft);
    cvec waveform_temp(Nfft);
    double power=0,temp;
    bool the_symbol_is_known;
    int i4; // Counter from symbols_in

    i4=0;
    symbol.zeros();

    for (int i1=0;i1<(int) Ns;i1++) {
      
      the_symbol_is_known=false;
      for (int i2=0;i2<known_symbol_pos.length();i2++) {
         if (i1==known_symbol_pos(i2))
           the_symbol_is_known=true;
      };

          

      if (the_symbol_is_known) {
        
         pilot_pattern(symbol,precoder);

         
      } else {  

        for (uint32_t i3=0;i3<length_ix;i3++) {
          symbol(ix(i3))=symbols_in(i4)*precoder(ix_sub(i3)); 
          i4++;
        };
        if (do_use_pilot_carriers) {
                symbol[pilot_carrier1]=(1.0+1.0i)/sqrt(2);
                symbol[pilot_carrier2]=(1.0+1.0i)/sqrt(2);
                symbol[pilot_carrier1]*=precoder(pilot_carrier1_sub);
                symbol[pilot_carrier2]*=precoder(pilot_carrier2_sub);
        };
      };          
        
      
      //waveform_temp=ifft(symbol,Nfft);
      pz_ifft(symbol,waveform_temp);
      waveform_temp=concat(waveform_temp.mid(Nfft-Np,Np),waveform_temp);
      waveform_temp=filter(pulse_shape,1,waveform_temp);

      if (rms_power!=-10) {
        temp=norm(waveform_temp.mid(Np,Nfft));
        power=power+temp*temp;
        temp=norm(waveform_temp);
        power=power+temp*temp;
      };

      int replace_pos=d_reorder(i1)*(Nfft+Np)+offset-pulse_shape_delay;

      if (replace_pos<0) {
        std::cout << "writing outside of waveform in OFDM1::modulate \n";
        exit(1);
      };
      
      for (int i2=0;i2<waveform_temp.size();i2++) {
         waveform(i2+replace_pos)+=waveform_temp(i2);
      };

      

    }; 
    if (rms_power!=-10) {
      power=power/(Ns*Nfft);
      rms_power=sqrt(power);
    };
    

};




void OFDM1::insert_pilot(cvec &waveform, uint32_t symbol_position, uint32_t offset, 
cvec precoder) {
  cvec symbol(Nfft);
  cvec waveform_temp(Nfft);
  symbol.zeros();
  pilot_pattern(symbol,precoder);
  //waveform_temp=ifft(symbol,Nfft);
  pz_ifft(symbol,waveform_temp);

  waveform.replace_mid(symbol_position*(Nfft+Np)+Np+offset,waveform_temp);
  waveform.replace_mid(symbol_position*(Nfft+Np)+offset,waveform_temp.mid(Nfft-Np,Np));
};


void OFDM1::estimate_channel(cvec waveform,double  &noise_variance_normalized, 
                             cvec &h, uint32_t start_sample) {
  double noise_variance;
  estimate_channel(waveform,noise_variance_normalized,h,start_sample,noise_variance);

};



void OFDM1::estimate_channel(cvec waveform,double  &noise_variance_normalized,
                             cvec &h, uint32_t start_sample, double &noise_variance) {

    cvec pilot(Nfft);
    int i1,counter;
    complex <double> hhat, actual_meas, ideal_meas, error;
    double error_div_h;
    cvec signal(Nfft);
    cvec h_rotated(Nfft);
    complex<double> rotation;



    //signal=fft(waveform.mid(0,Nfft));
    pz_fft(waveform.mid(start_sample,Nfft),signal);
    pilot.set_length(Nfft);
    pilot_pattern(pilot,ones_c(length_ix_all));
    h=signal(ix_all);    
    h/=pilot(ix_all);



    noise_variance_normalized=0;
    noise_variance=0;
    counter=0;

    for (i1=1;i1<ix_all.length()-1;i1++) {

      if (((ix_all(i1)-ix_all(i1-1))==1) && ((ix_all(i1+1)-ix_all(i1))==1)) {

        hhat=0.5*abs(h(i1-1))+0.5*abs(h(i1+1));
        complex<double> t;
        t=h(i1+1)/h(i1-1);
        double a=arg(h(i1-1))+0.5*arg(t);
        hhat=polar(abs(hhat),a);
        ideal_meas=hhat*pilot(ix_all(i1));

        actual_meas=signal(ix_all(i1)); 
        error=ideal_meas-actual_meas;
                
        noise_variance+=abs(error)*abs(error);
        error_div_h=abs(error)/abs(h(i1));
        noise_variance_normalized+=error_div_h*error_div_h;
        counter++;
      }
    };
    noise_variance_normalized/=counter;
    noise_variance/=counter;



 };
 


void OFDM1::init_multi_antenna(uint32_t no_ant, int32_t buffer_size, ivec &interferers_known_pos_in) {
 
  for (uint32_t i1=0;i1<Ns;i1++) {
     x[i1].set_size(no_ant,length_ix_all);
  };
  he.set_size(no_ant,length_ix_all);
  interferers_known_pos=interferers_known_pos_in;
  for (int i1=0;i1<interferers_known_pos.length();i1++) {
    hi[i1].set_size(no_ant,length_ix_all);
  };

  int num_payload_symbols=0;
  bool the_symbol_is_known;
  for (int i1=0;i1<(int) Ns;i1++) {
    the_symbol_is_known=false;
    for (int i2=0;i2<known_symbol_pos.length();i2++) {
         if (i1==known_symbol_pos(i2))
           the_symbol_is_known=true;
    };

    if (!(the_symbol_is_known)) {
        num_payload_symbols++;
    };    
  };

  x_c.set_size(length_ix_all,num_payload_symbols);


}


void OFDM1::demodulate_multi_antenna(const cmat &waveform,vec &soft_bit, uint32_t start_sample, int modulation_order){

    cvec waveform_temp(Nfft);
    double noise_variance, noise_variance_norm;
    vec noise_variance_vec; // Noise variance estmimated during pilot.
    vec noise_variance_SINR; // Noise variance estmimated when there is
                             // no signal. Only for benchmarking
    vec soft_temp;
    cvec symbolsc(Nfft);
    bool the_symbol_is_known;
    int32_t num_payload_symbols;
    double D=0.0, DI=0.0;
    
    noise_variance_vec.set_length(waveform.rows());
    // Estimate noise variance for analysis purposes

    if (estimate_SINR) {
      noise_variance_SINR.set_length(waveform.rows());
      int Nsymb;
      int start_ix=noise_estimation_start_ix;
      int stop_ix;
      
      for (int i2=0;i2<waveform.rows();i2++) {
        noise_variance_SINR(i2)=0;
        start_ix=noise_estimation_start_ix;
        stop_ix=start_ix+Nfft;
        Nsymb=0;
        while (stop_ix<noise_estimation_stop_ix) {
          waveform_temp=(waveform.get_row(i2)).mid(start_ix,Nfft);
          pz_fft(waveform_temp,symbolsc);
          noise_variance_SINR(i2)+=sum(sqr(symbolsc(ix_all)))/length_ix_all;
          start_ix=stop_ix+1;
          stop_ix=start_ix+Nfft;
          Nsymb++;
        };
        noise_variance_SINR(i2)=noise_variance_SINR(i2)/Nsymb;
      };
    };

    // Do FFT on all signals except pilot
    int i10=0;
    for (int i1=0;i1<(int) Ns;i1++) {
      the_symbol_is_known=false;
      for (int i2=0;i2<known_symbol_pos.length();i2++) {         
        if (i1==known_symbol_pos(i2))
           the_symbol_is_known=true;
      };

      if (!(the_symbol_is_known)) {
        for (int i2=0;i2<waveform.rows();i2++) {
           waveform_temp=(waveform.get_row(i2)).mid((Nfft+Np)*d_reorder(i1)+start_sample,Nfft);
           //symbolsc=fft(waveform_temp);
           pz_fft(waveform_temp,symbolsc);
           x[i10].set_row(i2,symbolsc(ix_all));
        };
        i10++;
      };          
    };
    num_payload_symbols=i10;

    if (modulation_order==0) return;

    // Estimate the vector channel of desired signal and interferers
    for (int i2=0;i2<he.rows();i2++) {
       double dummy;

       int s=start_of_known_symbol()(0)+start_sample;
       waveform_temp=(waveform.get_row(i2)).mid(s,Nfft);
       estimate_channel(waveform_temp,noise_variance_norm,
                        h,0,noise_variance);

       noise_variance_vec(i2)=noise_variance;
       he.set_row(i2,h);
       for (int i1=0;i1<interferers_known_pos.length();i1++) {
         s=start_sample+interferers_known_pos(i1)*(Nfft+Np);
         waveform_temp=(waveform.get_row(i2)).mid(s,Nfft);
         estimate_channel(waveform_temp,dummy,h,0);
         hi[i1].set_row(i2,h);
       };
    };

    noise_variance=mean(noise_variance_vec);

    if (modulation_order==1) return;


    cmat Rs; // Structured covariance matrix
    cmat Ru; // Un-structured covariance matrix
    cmat temp;
    cvec temp1;
    complex<double> hc;
    cvec w(waveform.rows());
    

    Rs.set_size(waveform.rows(),waveform.rows());
    Ru.set_size(waveform.rows(),waveform.rows());

    temp.set_size(waveform.rows(),waveform.rows());
    temp1.set_length(waveform.rows());

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

      Ru.zeros();
      for (int i1=0;i1<num_payload_symbols;i1++) {
        temp1=(x[i1]).get_col(i2);
        Ru+=outer_product(temp1,temp1,true);
      };
      Rs=outer_product(he.get_col(i2),he.get_col(i2),true);
      for (int i1=0;i1<interferers_known_pos.length();i1++) {
        Rs=Rs+outer_product(hi[i1].get_col(i2),hi[i1].get_col(i2),true);
      };
      Rs=Rs+diag(noise_variance_vec);
      w=inv(Rs)*he.get_col(i2);
      w=w/norm(w,2); 

      if (estimate_SINR) {
        double d;
        DI+=abs(elem_mult_sum(conj(w),(Rs+
                         diag(noise_variance_SINR-noise_variance_vec))*w));
        d=abs(elem_mult_sum(conj(w),he.get_col(i2)));
        D+=d*d;
      };
      

      h(i2)=elem_mult_sum(he.get_col(i2),conj(w));
      for (int i1=0;i1<num_payload_symbols;i1++) {
        x_c(i2,i1)=elem_mult_sum(x[i1].get_col(i2),conj(w));
      };

    };

    if (estimate_SINR) {
       SINR=D/(DI-D);
    };

    if (modulation_order==2) return;

    #if 0    
    cvec symbol;
    symbol.set_length(length_ix_all);
    soft_temp.set_length(length_ix);
    int bits_per_symbol=log2(modulation_order);

    // Re-estimate noise variance
    double noise_variance_temp;
    noise_variance=0;
    for (int i1=0;i1<num_payload_symbols;i1++) {
       symbol=x_c.get_col(i1);
       noise_variance_temp=0;
       soft_output(symbol,noise_variance_temp,modulation_order,soft_temp);
       noise_variance+=noise_variance_temp;
    };
    noise_variance/=num_payload_symbols;


    // Soft-value extraction
    for (int i1=0;i1<num_payload_symbols;i1++) {
       symbol=x_c.get_col(i1);
       soft_output(symbol,noise_variance,modulation_order,soft__temp);
       soft_bit.replace_mid(i1*length_ix*bits_per_symbol,soft_temp);
    };
    #endif
};

void OFDM1::extract_multi_antenna(const cmat &waveform,
                                  uint32_t start_sample){

  // When calling this function the following members will be set
  // x_c : A matrix which contains the combined signal (#subcarrier,#symbol)
  // h : Channel estimate including combiner #subcarriers
  // x[#symbols]: The raw complex samples. Each matrix (#subcarrier,#symbol).
  // he: Channel estimate for the desired signal (#antenna,#subcarrier)
  // hi[#interferers]: Channel estimate of interfering signals.
  
  vec dummy;

  demodulate_multi_antenna(waveform,dummy,start_sample,2);
};


void OFDM1::extract_multi_antenna(const cmat &waveform,
                                  uint32_t start_sample, cvec &symbols_out){

  // As the short version but also including equalization and estimation of the noise
  // in the equalized samples
  
  vec dummy;
  cvec signal(length_ix_all);
  std::complex<double> cpec;

  demodulate_multi_antenna(waveform,dummy,start_sample,2);


  uint32_t i3=0;  
  //std::cout << "h=" << h << "; \n";
  //std::cout << "x_c=" << x_c << "; \n";

  for (int i1=0;i1<x_c.cols();i1++) {

    elem_div_out( x_c.get_col(i1),h,signal);

    if (do_use_pilot_carriers) {
          // Common phase error correction (CPEC)
          //cpec=pilot_symbol/(signal(pilot_carrier1_sub)/h(pilot_carrier1_sub)+
          //         signal(pilot_carrier2_sub)/h(pilot_carrier2_sub));

        cpec=pilot_symbol/signal(pilot_carrier1_sub)+
          pilot_symbol/signal(pilot_carrier2_sub);      
        cpec=cpec/abs(cpec);
        signal=signal*cpec;

    };

    for (uint32_t i2=0;i2<length_ix;i2++) {
      symbols_out(i3)=signal(ix_sub(i2));
      i3=i3+1;
    };
  };

  
};

#if 0
void OFDM1::extract_multi_antenna(const cmat &waveform,
                                  uint32_t start_sample, cvec &symbols_out){

  // As the short version but also including equalization and estimation of the noise
  // in the equalized samples
  
  //QAM modulation;
  vec dummy;
  //cvec symbol(1);
  //cvec closest_constellation_point(1);

  //modulation.set_M(modulation_order);  
  demodulate_multi_antenna(waveform,dummy,start_sample,2);


  uint32_t i3=0;
  cvec temp(length_ix);

  for (int i1=0;i1<x_c.cols();i1++) {
      elem_div_out( x_c.get_col(i1),h(ix_sub),temp);
      for (uint32_t i2=0;i2<length_ix;i2++) {
        symbols_out(i3)=temp(i2);
        i3=i3+1;
      };
  };


};
#endif



void OFDM1::demodulate(const cvec &waveform,cvec &symbols_out, uint32_t start_sample, const cvec &h,  double &rms_signal)
{
   double power,temp;
    cvec symbolsc(Nfft);
    vec soft_temp;
    bool the_symbol_is_known;
    cvec signal(Nfft), demodulator_input;
    int i4; // Counter in symbols_out
    std::complex<double> cpec;
    int s;

    demodulator_input.set_length(length_ix);

    i4=0;
    power=0.0;
    cpec=1.0;
    
    for (int i1=0;i1<(int) Ns;i1++) {
        
      
       the_symbol_is_known=false;
       for (int i2=0;i2<known_symbol_pos.length();i2++) {
         if (i1==known_symbol_pos(i2))
           the_symbol_is_known=true;
       };

      if (!the_symbol_is_known) {

        s=(Nfft+Np)*d_reorder(i1)+start_sample;
        
        symbolsc=waveform.mid(s,Nfft);
        temp=norm(symbolsc);
        power=power+temp*temp;

        pz_fft(symbolsc,signal);

        //std::cout << "h1=" << h(pilot_carrier1_sub) << "\n";
        //std::cout << "s1=" << signal(pilot_carrier1) << "\n";
        //std::cout << "h2=" << h(pilot_carrier2_sub) << "\n";
        //std::cout << "s2=" << signal(pilot_carrier2) << "\n";
        //std::cout << "pilot_symbol=" << pilot_symbol << "\n";

        if (do_use_pilot_carriers) {
          // Common phase error correction (CPEC)
          //cpec=pilot_symbol/(signal(pilot_carrier1_sub)/h(pilot_carrier1_sub)+
          //         signal(pilot_carrier2_sub)/h(pilot_carrier2_sub));

          cpec=pilot_symbol/(signal(pilot_carrier1)/h(pilot_carrier1_sub)+
                     signal(pilot_carrier2)/h(pilot_carrier2_sub));

          cpec=cpec/abs(cpec);

        };
            
        elem_div_out(signal(ix_all),h,demodulator_input);


        if (do_use_pilot_carriers) 
          demodulator_input=demodulator_input*cpec;


        for (uint32_t i5=0;i5<length_ix;i5++) {
          symbols_out(i4)=demodulator_input(ix_sub(i5));
          i4++;
        };
      };  
    };
    power=power/((Ns-known_symbol_pos.length())*Nfft);
    rms_signal=sqrt(power);
}


void OFDM1::synchronize(cvec &waveform, uint32_t &start_pos,double &f_offset) {
    cvec de_spread(train.size());
    
    uint32_t Nlags;
    double max_value;
    int start_pos_int, max_index;

    Nlags=waveform.length()-train.length()+1;
    vec FrequencyIx(Nlags), Criterions(Nlags);

    for (uint32_t i1=0; i1<Nlags;i1 ++) {
      de_spread=elem_mult(waveform.mid(i1,train.length()),conj(train));
      max_value=max(abs(fft(de_spread,Nfft_sync)),max_index);
      FrequencyIx(i1)=max_index;
      Criterions(i1)=max_value;
    };

    max_value=max(Criterions,start_pos_int);
    start_pos=start_pos_int;
    f_offset=FrequencyIx(start_pos)/Nfft_sync;
    if (f_offset>0.5)
      f_offset=f_offset-1;

    f_offset=f_offset*25e6;
    
 
}

void OFDM1::synchronize(cvec &waveform, uint32_t &start_pos) {
    cvec de_spread(train.size());
    
    uint32_t Nlags;
    double max_value;
    int start_pos_int;

    Nlags=waveform.length()-train.length()+1;
    vec Criterions(Nlags);

    for (uint32_t i1=0; i1<Nlags;i1 ++) {
      de_spread=elem_mult(waveform.mid(i1,train.length()),conj(train));
      max_value=sum(abs(de_spread));
      Criterions(i1)=max_value;
    };

    max_value=max(Criterions,start_pos_int);
    start_pos=start_pos_int;
    
 
}

void OFDM1::interleave_symbols(cvec &symbols, int index) {

  cvec symbol;
  
  index=index % symbol_interleave_matrix.rows();
  symbol.set_length(length_ix);
  
  for (uint32_t i1=0;i1<Ns-known_symbol_pos.length();i1++) {
    symbol=symbols.mid(i1*length_ix,length_ix);
    symbol=symbol(symbol_interleave_matrix.get_row(index));
    symbols.set_subvector(i1*length_ix,symbol);
  };
 
};

void OFDM1::de_interleave_symbols(cvec &symbols, int index) {

  cvec symbol;
  
  index=index % symbol_de_interleave_matrix.rows();
  symbol.set_length(length_ix);
  
  for (uint32_t i1=0;i1<Ns-known_symbol_pos.length();i1++) {
    symbol=symbols.mid(i1*length_ix,length_ix);
    symbol=symbol(symbol_de_interleave_matrix.get_row(index));
    symbols.set_subvector(i1*length_ix,symbol);
  };
  

};