/*---------------------------------------------------------------------------*
 *                                   IT++			             *
 *---------------------------------------------------------------------------*
 * Copyright (c) 1995-2004 by Tony Ottosson, Thomas Eriksson, Pål Frenger,   *
 * Tobias Ringström, and Jonas Samuelsson.                                   *
 *                                                                           *
 * Permission to use, copy, modify, and distribute this software and its     *
 * documentation under the terms of the GNU General Public License is hereby *
 * granted. No representations are made about the suitability of this        *
 * software for any purpose. It is provided "as is" without expressed or     *
 * implied warranty. See the GNU General Public License for more details.    *
 *---------------------------------------------------------------------------*/

/*! 
  \file 
  \brief Implementation of Communication Channel classes and functions.
  \author Tony Ottosson and Pål frenger
  
  1.33
  
  2004/10/28 12:59:55
*/

#include "comm/channel.h"

#include "base/binary.h"
#include "base/stat.h"
#include "base/matfunc.h"
#include "base/scalfunc.h"
#include "base/elmatfunc.h"
#include "base/specmat.h"
#include "base/filter.h"
#include "base/transforms.h"
#include "base/bessel.h"

using std::complex;

namespace itpp {

  vec jake_filter(double NormFDopp, int order)
  {
    int i, L = (int)std::floor(double(order)/2);
    double t, x0;
    vec x_pos(L), x_neg(L), x(2*L+1), h(2*L+1);
    for (i=1; i<L+1; i++) {	
      t = double(i);
      x_pos(i-1) = besselj(0.25, 2*pi*NormFDopp*t) / std::sqrt(std::sqrt(t));
    }
    x0 = double(1.468813) * std::sqrt(std::sqrt(NormFDopp));
    x_neg = reverse(x_pos);
    x = concat( concat(x_neg, x0), x_pos);
    h = elem_mult(hamming(2*L+1), x);
    h /= norm(h);
    return h;
  }



  // ------------------------------------------------------------------------------------------------------------------
  Fading_Generator::Fading_Generator(const double norm_doppler, const Doppler_Spectrum spectrum)
  {
    set_norm_doppler(norm_doppler);
    set_doppler_spectrum(spectrum);
    los_power = 0.0; // no LOS component
    los_dopp = 0.0;
  }

  void Fading_Generator::set_norm_doppler(const double norm_doppler)
  {
    it_assert(norm_doppler >=0 && norm_doppler <1.0, "Fading_Generator: Normalized Doppler must be >=0 and <1.");
    n_dopp = norm_doppler; 
    init_flag = false;
  }
  
  void Fading_Generator::set_doppler_spectrum(const Doppler_Spectrum spectrum)
  {
    dopp_spectrum = spectrum;
    if (spectrum != Rice) { // no LOS component
      los_power = 0.0;
      los_dopp = 0.0;
    }
    init_flag = false;
  }

  void Fading_Generator::set_LOS(const double relative_power, const double relative_doppler)
  {
    it_assert(relative_doppler >=0 && relative_doppler <=1.0, "Relative Doppler must be >=0 and <=1.");
    it_assert(relative_power >= 0.0, "Rice factor need to be >= 0.0");
    it_assert(dopp_spectrum == Rice, "Can only set LOS component if Rice spectrum");
    los_power = relative_power;
    los_dopp = relative_doppler;
    init_flag = false;
  }

  cvec Fading_Generator::generate(const int no_samples)
  {
    cvec output;
    this->generate(no_samples, output);
    return output;
  }

  void Fading_Generator::generate_zero_doppler(const int no_samples, cvec &output)
  {
    output = randn_c(no_samples);
    if(los_power > 0.0) {
      double diffuse = std::sqrt(1.0/(1.0+los_power));
      double direct = diffuse*std::sqrt(los_power);
      for (int i=0; i<no_samples; i++)
	output(i) = diffuse*output(i) + direct*complex<double>(std::cos(2*pi*los_dopp*n_dopp*(i+time_offset)),std::sin(2*pi*los_dopp*n_dopp*(i+time_offset)));
    }
    time_offset += no_samples; 
  }

  // ------------------------------------------------------------------------------------------------------------------
  Rice_Fading_Generator::Rice_Fading_Generator(const double norm_doppler, const Doppler_Spectrum spectrum, const int no_freq, const Rice_Method method)
    : Fading_Generator(norm_doppler, spectrum)
  {
    set_no_frequencies(no_freq);
    set_method(method);
  }

  void Rice_Fading_Generator::set_no_frequencies(const int no_freq)
  {
    it_assert(no_freq >= 7, "Rice_Fading_Generator: number of doppler frequencies should be at least 7");
    Ni = no_freq;
    init_flag = false;
  }

  void Rice_Fading_Generator::set_method(const Rice_Method method)
  {
    // check if this method works for the given spectrum
    rice_method = method;
    init_flag = false;
  }

  int Rice_Fading_Generator::get_no_frequencies()
  {
    return Ni;
  }

  Rice_Method Rice_Fading_Generator::get_method()
  {
    return rice_method;
  }

  void Rice_Fading_Generator::init()
  {
    switch(rice_method) {
    case MEDS: // Method of Exact Doppler Spread (MEDS)
      init_MEDS();
      break;

    default:
      it_error("Rice_Fading_Generator: this method is not implemented");
    };

    time_offset = 0; // state in the process
    init_flag = true; // generator ready to use
  }

  void Rice_Fading_Generator::set_time_offset(const int offset)
  {
    it_assert(offset >= 0, "Rice_Fading_Generator: time offset need to be >=0");
    time_offset = offset;
  }
  
  int Rice_Fading_Generator::set_time_offset()
  {
    return time_offset;
  }

  void Rice_Fading_Generator::generate(const int no_samples, cvec &output)
  {
    if (init_flag == false)
      init();

    if (n_dopp == 0.0)
      generate_zero_doppler(no_samples, output);
    else {
      output.set_size(no_samples, false);

      for (int i=0; i<no_samples; i++) {
	output(i) = complex<double>( sum( elem_mult( c1, cos(2*pi*f1*n_dopp*(i+time_offset)+th1) ) ), 
				     sum( elem_mult( c2, cos(2*pi*f2*n_dopp*(i+time_offset)+th2) ) ) );
      }
      
      if(los_power > 0.0) { // LOS component exist
	double diffuse = std::sqrt(1.0/(1.0+los_power));
	double direct = diffuse*std::sqrt(los_power);
	for (int i=0; i<no_samples; i++)
	  output(i) = diffuse*output(i) + direct*complex<double>(std::cos(2*pi*los_dopp*n_dopp*(i+time_offset)),std::sin(2*pi*los_dopp*n_dopp*(i+time_offset)));
      }
      time_offset += no_samples; 
    }
  }


  void Rice_Fading_Generator::init_MEDS()
  {
    vec n;

    switch(dopp_spectrum) {
    case Jakes: // Jakes (Clark) spectrum
    case Rice: // Rice also have a Jakes spectrum component
      n = linspace(1,Ni,Ni);
      f1 = sin(pi/(2*Ni)*(n-0.5));
      c1 = std::sqrt(1.0/Ni)*ones(Ni);
      th1 = randu(Ni)*2*pi;

      n = linspace(1,Ni+1,Ni+1);
      f2 = sin(pi/(2*(Ni+1))*(n-0.5));
      c2 = std::sqrt(1.0/(Ni+1))*ones(Ni+1);
      th2 = randu(Ni+1)*2*pi;
      break;

    default:
      it_error("Rice_Fading_Generator: this spectrum is not implemented for the MEDS Rice fading generator");
    };
  }

  // ------------------------------------------------------------------------------------------------------------------

  FIR_Fading_Generator::FIR_Fading_Generator(const double norm_doppler, const Doppler_Spectrum spectrum, const int filter_length)
    : Fading_Generator(norm_doppler, spectrum)
  {
    set_filter_length(filter_length);
  }


  void FIR_Fading_Generator::set_filter_length(const int filter_length)
  {
    it_assert(filter_length >= 50, "FIR_Fading_Generator: filter length should be at least 50");
    fir_length = filter_length;
    init_flag = false;
  }

  int FIR_Fading_Generator::get_filter_length()
  {
    return fir_length;
  }

  void FIR_Fading_Generator::init()
  {
    double norm_dopp = n_dopp;
    if (n_dopp > 0.0) {
      switch(dopp_spectrum) {
      case Jakes:
      case Rice: // Rice also has a Jakes spectrum component
	upsample_rate = 1; // upsampling rate      
	while (norm_dopp < 0.1) { // calculate a reasonable upsample rate so that normalized doppler is > 0.1
	  norm_dopp *= 2;
	  upsample_rate *= 2;
	}
	fir_filter.set_coeffs(jake_filter(norm_dopp, fir_length));
	break;

      default:
	it_error("FIR_Fading_Generator: doppler spectrum is not implemented");
      };

      // fill filter with dummy data
      cvec dummy = fir_filter(randn_c(fir_length));
      
      left_overs.set_size(0, false);
    }
    init_flag = true; // generator ready to use
  }

  void FIR_Fading_Generator::generate(const int no_samples, cvec &output)
  {
    if (init_flag == false)
      init();

    if (n_dopp == 0.0)
      generate_zero_doppler(no_samples, output);
    else {
      int no_upsamples = ceil_i(double(no_samples-left_overs.size())/double(upsample_rate)) + 1; // calculate number of samples before upsampling
      
      // should make a smarter interpolation here!!!
      lininterp(fir_filter(randn_c(no_upsamples)), upsample_rate, output);
      output = concat(left_overs, output); // add left-overs from previous filtering
      left_overs = output.right(output.size()-no_samples); // save left-overs for next round of filtering
      output.set_size(no_samples, true);

      if(los_power > 0.0) { // LOS component exist
	double diffuse = std::sqrt(1.0/(1.0+los_power));
	double direct = diffuse*std::sqrt(los_power);
	for (int i=0; i<no_samples; i++)
	  output(i) = diffuse*output(i) + direct*complex<double>(std::cos(2*pi*los_dopp*n_dopp*(i+time_offset)),std::sin(2*pi*los_dopp*n_dopp*(i+time_offset)));
      }
      time_offset += no_samples; 
    }
  }

  // ------------------------------------------------------------------------------------------------------------------
  IFFT_Fading_Generator::IFFT_Fading_Generator(const double norm_doppler, const Doppler_Spectrum spectrum)
    : Fading_Generator(norm_doppler, spectrum) { }


  void IFFT_Fading_Generator::init()
  {
    init_flag = true; // generator ready to use
  }


  void IFFT_Fading_Generator::generate(const int no_samples, cvec &output)
  {
    if (init_flag == false)
      init();

    if (n_dopp == 0.0)
      generate_zero_doppler(no_samples, output);
    else {
      switch(dopp_spectrum) {

      case Jakes:
      case Rice: // Rice also has a Jakes spectrum component
	generate_Jakes(no_samples, output);
	break;
      default:
	it_error("IFFT_Fading_Generator: doppler spectrum is not implemented");
      };

      if(los_power > 0.0) { // LOS component exist
	double diffuse = std::sqrt(1.0/(1.0+los_power));
	double direct = diffuse*std::sqrt(los_power);
	for (int i=0; i<no_samples; i++)
	  output(i) = diffuse*output(i) + direct*complex<double>(std::cos(2*pi*los_dopp*n_dopp*(i+time_offset)),std::sin(2*pi*los_dopp*n_dopp*(i+time_offset)));
      }
      time_offset += no_samples; 
    }
  }

  void IFFT_Fading_Generator::generate_Jakes(const int no_samples, cvec &output)
  {
    int Nfft = pow2i(needed_bits(no_samples));
    double df = 1.0/Nfft;
    int noisesamp = (int)std::ceil(n_dopp/df);
		
    while (noisesamp <= 10) { // if too few samples, increase the FFT size
      Nfft *= 2;
      df = 1.0/double(Nfft);
      noisesamp = (int)std::ceil(n_dopp/df);
    }

    vec Fpos = linspace(0,0.5,Nfft/2+1);
    vec F = concat(Fpos, reverse(-Fpos(1,Nfft/2-1)));
    vec S = zeros(Nfft);
	
    for (int i=0; i<F.size(); i++) {
      if (std::abs(F(i)) < n_dopp)
	S(i) = std::sqrt(1.5 / ( pi*n_dopp*std::sqrt(1-std::pow(F(i)/n_dopp,2))));
      else if (std::abs(F(i)) == n_dopp)
	S(i) = 1000000;
    }

    S /= norm(S,2); S *= Nfft;

    cvec x(Nfft);

    // lots of zeros. Not necessary to do the multiplication for all elements!!!
    // S is converted. Not necessary???
    x = ifft( elem_mult(to_cvec(S), concat(randn_c(noisesamp), zeros_c(Nfft-2*noisesamp), randn_c(noisesamp)) ) ); 
    output = x.mid(0, no_samples);   
  }

  // ------------------------------------------------------------------------------------------------------------------
  Channel_Specification::Channel_Specification(const vec &avg_power_dB, const vec &delay_prof)
  {
    set_channel_profile(avg_power_dB, delay_prof);
  }

  Channel_Specification::Channel_Specification(const Channel_Profile profile)
  {
    set_channel_profile(profile);
  }

  void Channel_Specification::set_channel_profile(const vec &avg_power_dB, const vec &delay_prof)
  {
    it_assert(min(delay_prof) == 0, "Minimum relative delay must be 0.");
    it_assert(avg_power_dB.size() == delay_prof.size(), "Power and delay vectors must be of equal length!");
    it_assert(delay_prof(0) == 0, "First tap must be at zero delay");


    N_taps = delay_prof.size();
    // taps should be sorted and unique
    for (int i=1; i<N_taps; i++) {
      it_assert(delay_prof(i) > delay_prof(i-1), "Delays should be sorted and unique");
    }

    a_prof_dB = avg_power_dB;
    d_prof = delay_prof;

    // set doppler spectrum to Jakes per default
    tap_doppler_spectrum.set_size(N_taps, false);
    for (int i=0; i<N_taps; i++)
      tap_doppler_spectrum(i) = Jakes;

    discrete = false;
  }

  void Channel_Specification::set_channel_profile(const Channel_Profile profile)
  {
    switch (profile) {
      // -------------- ITU Channel models -----------------
    case ITU_Vehicular_A:
      set_channel_profile( vec("0 -1 -9 -10 -15 -20"), ivec("0 310 710 1090 1730 2510") * 1e-9 );
      break;

    case ITU_Vehicular_B:
      set_channel_profile( vec("-2.5 0 -12.8 -10 -25.2 -16"), ivec("0 300 8900 12900 17100 20000") * 1e-9 );
      break;
      
    case ITU_Pedestrian_A:
      set_channel_profile( vec("0 -9.7 -19.2 -22.8"), ivec("0 110 190 410") * 1e-9 );
      break;
 
    case ITU_Pedestrian_B:
      set_channel_profile( vec("0 -0.9 -4.9 -8 -7.8 -23.9"), ivec("0 200 800 1200 2300 3700") * 1e-9 );
      break;

      // -------------- COST259 Channel models -----------------
    case COST259_TUx:
      set_channel_profile( vec("-5.7 -7.6 -10.1 -10.2 -10.2 -11.5 -13.4 -16.3 -16.9 -17.1 -17.4 -19 -19 -19.8 -21.5 -21.6 -22.1 -22.6 -23.5 -24.3"), 
			   ivec("217 512 514 517 674 882 1230 1287 1311 1349 1533 1535 1622 1818 1836 1884 1943 2048 2140") * 1e-9 );
      break;

    case COST259_RAx:
      set_channel_profile( vec("-5.2 -6.4 -8.4 -9.3 -10 -13.1 -15.3 -18.5 -20.4 -22.4"), ivec("0 42 101 129 149 245 312 410 469 528") * 1e-9 );
      set_doppler_spectrum(0, Rice);
      set_LOS( sqr(0.91/0.41), 0.7); // What should the rice factor be??? Not sure from report!
      break;

    case COST259_HTx:
      set_channel_profile( vec("-3.6 -8.9 -10.2 -11.5 -11.8 -12.7 -13.0 -16.2 -17.3 -17.7 -17.6 -22.7 -24.1 -25.8 -25.8 -26.2 -29 -29.9 -30 -30.7"), 
			   ivec("356 441 528 546 609 625 842 916 941 15000 16172 16492 16876 16882 16978 17615 17827 17849 18016") * 1e-9 );
      break;

      // -------------- COST207 Channel models -----------------
    case COST207_RA:
      set_channel_profile( vec("0 -2 -10 -20"), ivec("0 200 400 600") * 1e-9 );
      set_doppler_spectrum(0, Rice);
      set_LOS( sqr(0.91/0.41), 0.7);
      break;

    case COST207_RA6:
      set_channel_profile( vec("0 -4 -8 -12 -16 -20"), ivec("0 100 200 300 400 500") * 1e-9 );
      set_doppler_spectrum(0, Rice);
      set_LOS( sqr(0.91/0.41), 0.7);
      break;

    case COST207_TU:
      set_channel_profile( vec("-3 0 -2 -6 -8 -10"), ivec("0 200 600 1600 2400 5000") * 1e-9 );
      set_doppler_spectrum(2, GaussI);
      set_doppler_spectrum(3, GaussI);
      set_doppler_spectrum(4, GaussII);
      set_doppler_spectrum(5, GaussII);
      break;

    case COST207_TU6alt:
      set_channel_profile( vec("-3 0 -2 -6 -8 -10"), ivec("0 200 500 1600 2300 5000") * 1e-9 );
      set_doppler_spectrum(3, GaussI);
      set_doppler_spectrum(4, GaussII);
      set_doppler_spectrum(5, GaussII);
      break;

    case COST207_TU12:
      set_channel_profile( vec("-4 -3 0 -2 -3 -5 -7 -5 -6 -9 -11 -10"), ivec("0 200 400 600 800 1200 1400 1800 2400 3000 3200 5000") * 1e-9 );
      set_doppler_spectrum(3, GaussI);
      set_doppler_spectrum(4, GaussI);
      set_doppler_spectrum(5, GaussI);
      set_doppler_spectrum(6, GaussI);
      set_doppler_spectrum(7, GaussI);
      set_doppler_spectrum(8, GaussII);
      set_doppler_spectrum(9, GaussII);
      set_doppler_spectrum(10, GaussII);
      set_doppler_spectrum(11, GaussII);
      break;

    case COST207_TU12alt:
      set_channel_profile( vec("-4 -3 0 -2.6 -3 -5 -7 -5 -6.5 -8.6 -11 -10"), ivec("0 200 400 600 800 1200 1400 1800 2400 3000 3200 5000") * 1e-9 );
      set_doppler_spectrum(4, GaussI);
      set_doppler_spectrum(5, GaussI);
      set_doppler_spectrum(6, GaussI);
      set_doppler_spectrum(7, GaussI);
      set_doppler_spectrum(8, GaussII);
      set_doppler_spectrum(9, GaussII);
      set_doppler_spectrum(10, GaussII);
      set_doppler_spectrum(11, GaussII);
      break;

    case COST207_BU:
      set_channel_profile( vec("-3 0 -3 -5 -2 -4"), ivec("0 400 1000 1600 5000 6600") * 1e-9 );
      set_doppler_spectrum(2, GaussI);
      set_doppler_spectrum(3, GaussI);
      set_doppler_spectrum(4, GaussII);
      set_doppler_spectrum(5, GaussII);
      break;

    case COST207_BU6alt:
      set_channel_profile( vec("-2.5 0 -3 -5 -2 -4"), ivec("0 300 1000 1600 5000 6600") * 1e-9 );
      set_doppler_spectrum(2, GaussI);
      set_doppler_spectrum(3, GaussI);
      set_doppler_spectrum(4, GaussII);
      set_doppler_spectrum(5, GaussII);
      break;

    case COST207_BU12:
      set_channel_profile( vec("-7 -3 -1 0 -2 -6 -7 -1 -2 -7 -10 -15"), ivec("0 200 400 800 1600 2200 3200 5000 6000 7200 8200 10000") * 1e-9 );
      set_doppler_spectrum(3, GaussI);
      set_doppler_spectrum(4, GaussI);
      set_doppler_spectrum(5, GaussII);
      set_doppler_spectrum(6, GaussII);
      set_doppler_spectrum(7, GaussII);
      set_doppler_spectrum(8, GaussII);
      set_doppler_spectrum(9, GaussII);
      set_doppler_spectrum(10, GaussII);
      set_doppler_spectrum(11, GaussII);
      break;

    case COST207_BU12alt:
      set_channel_profile( vec("-7.7 -3.4 -1.3 0 -2.3 -5.6 -7.4 -1.4 -1.6 -6.7 -9.8 -15.1"), ivec("0 100 300 700 1600 2200 3100 5000 6000 7200 8100 10000") * 1e-9 );
      set_doppler_spectrum(3, GaussI);
      set_doppler_spectrum(4, GaussI);
      set_doppler_spectrum(5, GaussII);
      set_doppler_spectrum(6, GaussII);
      set_doppler_spectrum(7, GaussII);
      set_doppler_spectrum(8, GaussII);
      set_doppler_spectrum(9, GaussII);
      set_doppler_spectrum(10, GaussII);
      set_doppler_spectrum(11, GaussII);
      break;


    case COST207_HT:
      set_channel_profile( vec("0 -2 -4 -7 -6 -12"), ivec("0 200 400 600 15000 17200") * 1e-9 );
      set_doppler_spectrum(4, GaussII);
      set_doppler_spectrum(5, GaussII);
      break;

    case COST207_HT6alt:
      set_channel_profile( vec("0 -1.5 -4.5 -7.5 -8 -17.7"), ivec("0 100 300 500 15000 17200") * 1e-9 );
      set_doppler_spectrum(4, GaussII);
      set_doppler_spectrum(5, GaussII);
      break;

    case COST207_HT12:
      set_channel_profile( vec("-10 -8 -6 -4 0 0 -4 -8 -9 -10 -12 -14"), ivec("0 200 400 600 800 2000 2400 15000 15200 15800 17200 20000") * 1e-9 );
      set_doppler_spectrum(3, GaussI);
      set_doppler_spectrum(4, GaussI);
      set_doppler_spectrum(5, GaussI);
      set_doppler_spectrum(6, GaussII);
      set_doppler_spectrum(7, GaussII);
      set_doppler_spectrum(8, GaussII);
      set_doppler_spectrum(9, GaussII);
      set_doppler_spectrum(10, GaussII);
      set_doppler_spectrum(11, GaussII);
      break;

    case COST207_HT12alt:
      set_channel_profile( vec("-10 -8 -6 -4 0 0 -4 -8 -9 -10 -12 -14"), ivec("0 100 300 500 700 1000 1300 15000 15200 15700 17200 20000") * 1e-9 );
      set_doppler_spectrum(4, GaussI);
      set_doppler_spectrum(5, GaussI);
      set_doppler_spectrum(6, GaussI);
      set_doppler_spectrum(7, GaussII);
      set_doppler_spectrum(8, GaussII);
      set_doppler_spectrum(9, GaussII);
      set_doppler_spectrum(10, GaussII);
      set_doppler_spectrum(11, GaussII);
      break;
    };
  }


  void Channel_Specification::set_doppler_spectrum(Doppler_Spectrum *tap_spectrum)
  {
    for (int i=0; i<N_taps; i++)
      tap_doppler_spectrum(i) = tap_spectrum[i];
  }

  void Channel_Specification::set_doppler_spectrum(const int tap_number, Doppler_Spectrum tap_spectrum)
  {
    tap_doppler_spectrum(tap_number) = tap_spectrum;
  }

  void Channel_Specification::set_LOS(const double relative_power, const double norm_doppler)
  {
    it_assert(N_taps >= 1, "Cannot set LOS component if not set channel profile");
    it_assert(norm_doppler >=0 && norm_doppler <=1.0, "Normalized Doppler must be >=0 and <=1.");
    it_assert(relative_power >= 0.0, "Rice factor need to be >= 0.0");
    it_assert(tap_doppler_spectrum(0) == Rice, "Can only set LOS component if Rice spectrum");

    los_power = relative_power;
    los_dopp = norm_doppler;
  }

  void Channel_Specification::get_channel_profile(vec &avg_power_dB, vec &delay_prof)
  {
    avg_power_dB = a_prof_dB;
    delay_prof = d_prof;
  }

  vec Channel_Specification::get_avg_power_dB()
  {
    return a_prof_dB;
  }

  vec Channel_Specification::get_delay_prof()
  {
    return d_prof;
  }

  Doppler_Spectrum Channel_Specification::get_doppler_spectrum(const int index)
  {
    it_assert( (index >=0) && (index < N_taps), "Channel_Specification: index of of range");
    return tap_doppler_spectrum(index);
  }

  double Channel_Specification::get_LOS_power()
  {
    return los_power;
  }

  double Channel_Specification::get_LOS_doppler()
  {
    return los_dopp;
  }

  double Channel_Specification::calc_mean_excess_delay()
  {
    vec a_prof = inv_dB(a_prof_dB);
    vec delay_prof = d_prof;

    if (discrete) {
      delay_prof *= discrete_Ts;
    }

    return ( a_prof*delay_prof / sum(a_prof));
  }

  double Channel_Specification::calc_rms_delay_spread()
  {
    vec a_prof = inv_dB(a_prof_dB);
    vec delay_prof = d_prof;

    if (discrete) {
      delay_prof *= discrete_Ts;
    }

    double a = ( a_prof*delay_prof / sum(a_prof));
    double b = ( a_prof*sqr(delay_prof) / sum(a_prof) );

    return ( std::sqrt(b-a*a) );
  }

  void Channel_Specification::discretize(const double Ts)
  {
    it_assert(N_taps > 0, "Channel_Specification::discretize: no channel profile specified");
    vec delay_prof(N_taps);
    vec power(N_taps);

    int j = 0, no_taps, j_delay = 0;

    vec a_prof = inv_dB(a_prof_dB); // Convert power profile


    it_assert(d_prof(0) == 0, "Channel_Specification: first tap should be at zero delay");
    delay_prof(0) = d_prof(0);
    power(0) = a_prof(0);

    // Taps within ( (j-0.5)Ts,(j+0.5)Ts] are included in the jth tap
    // 
    for(int i=1; i<N_taps; i++) {
      if( d_prof(i) > (j_delay+0.5)*Ts ) {
	while( d_prof(i) > (j_delay+0.5)*Ts ) { j_delay++; } // first skip empty taps
	// create a new tap at (j+1)Ts
	j++;
	delay_prof(j) = j_delay;
	power(j) = a_prof(i);
      } else {
	power(j) += a_prof(i);
      }
    }

    no_taps = j+1; // number of taps found
    power.set_size(no_taps, true);
    delay_prof.set_size(no_taps, true);
    
    // write over existing channel profile with the discretized version
    a_prof_dB = dB(power);
    d_prof = delay_prof;
    N_taps = no_taps;

    discrete = true;    
    discrete_Ts = Ts;

  }

  // ------------------------------------------------------------------------------------------------------------------
  TDL_Channel::TDL_Channel(const double norm_doppler, const vec &avg_power_dB, const ivec &delay_prof)
  {
    set_channel_profile(avg_power_dB, delay_prof);
    set_norm_doppler(norm_doppler);
  }

  TDL_Channel::TDL_Channel(Channel_Specification &channel_spec)
  {
    set_channel_profile(channel_spec);
    
    // set doppler spectrum
    tap_doppler_spectrum.set_size(N_taps, false);
    for (int i=0; i<N_taps; i++) {
      tap_doppler_spectrum(i) = channel_spec.get_doppler_spectrum(i);
    }

    if (tap_doppler_spectrum(0) == Rice) // set LOS component
      set_LOS(channel_spec.get_LOS_power(), channel_spec.get_LOS_doppler());
  }

  void TDL_Channel::set_channel_profile(const vec &avg_power_dB, const ivec &delay_prof)
  {
    it_assert(min(delay_prof) == 0, "Minimum relative delay must be 0.");
    it_assert(avg_power_dB.size() == delay_prof.size(), "Power and delay vectors must be of equal length!");
    it_assert(delay_prof(0) == 0, "First tap must be at zero delay");

    N_taps = delay_prof.size();
    // taps should be sorted and unique
    for (int i=1; i<N_taps; i++) {
      it_assert(delay_prof(i) > delay_prof(i-1), "Delays should be sorted and unique");
    }

    a_prof = pow(10.0,avg_power_dB/20.0); // Convert power profile to ampiltude profile
    a_prof /= norm(a_prof); // Normalize

    d_prof = delay_prof;

    init_flag = false;
  }

  void TDL_Channel::set_norm_doppler(const double norm_doppler)
  {
    it_assert(norm_doppler >=0 && norm_doppler <=1.0, "TDL_Channel: Normalized Doppler must be >=0 and <=1.");
    n_dopp = norm_doppler;
    init_flag = false;
  }

  void TDL_Channel::set_channel_profile_uniform(const int no_taps)
  {
    it_assert(no_taps >= 1, "Minimum number of taps is 1.");

    vec avg_power_dB = zeros(no_taps);
    ivec delay_prof(no_taps);

    for (int i=0; i<no_taps; i++)
      delay_prof(i) = i;

    set_channel_profile(avg_power_dB, delay_prof);
  }

  // not implemented
  void TDL_Channel::set_channel_profile_exponential(const double delay_spread)
  {
    it_assert(delay_spread > 0.0, "Delay spread must be larger than 0.");


    //a_prof /= norm(a_prof); // Normalize
  }


  void TDL_Channel::set_channel_profile(Channel_Specification &channel_spec)
  {
    vec avg_power_dB;
    vec delay_profile;
    
    it_assert(channel_spec.is_discrete() == true, "TDL_Channel: Channel_Specification has not been discretized");

    channel_spec.get_channel_profile(avg_power_dB, delay_profile);
    set_channel_profile(avg_power_dB, to_ivec(delay_profile));
  }


  void TDL_Channel::set_doppler_spectrum(Doppler_Spectrum *tap_spectrum)
  {
    it_assert(N_taps > 0, "TDL_Channel:: set_doppler_spectrum: channel profile not defined yet");
    tap_doppler_spectrum.set_size(N_taps, false);

    for (int i=0; i<N_taps; i++)
      tap_doppler_spectrum = tap_spectrum[i];

    init_flag = false;
  }

  void TDL_Channel::set_generation_method(const Fading_Generation_Method generation_method)
  {
    method = generation_method;
    init_flag = false;
  }

  void TDL_Channel::set_LOS(const double relative_power, const double norm_doppler)
  {
    it_assert(N_taps >= 1, "Cannot set LOS component if not set channel profile");
    it_assert(norm_doppler >=0 && norm_doppler <=1.0, "Normalized Doppler must be >=0 and <=1.");
    it_assert(relative_power >= 0.0, "Rice factor need to be >= 0.0");
    it_assert(tap_doppler_spectrum(0) == Rice, "Can only set LOS component if Rice spectrum");

    los_power = relative_power;
    los_dopp = norm_doppler;
    init_flag = false;
  }

  void TDL_Channel::get_channel_profile(vec &avg_power_dB, ivec &delay_prof)
  {
    avg_power_dB = 20 * log10(a_prof);
    delay_prof = d_prof;
  }


  vec TDL_Channel::get_avg_power_dB()
  {
    return ( 20 * log10(a_prof) );
  }

  ivec TDL_Channel::get_delay_prof()
  {
    return d_prof;
  }


  double TDL_Channel::get_norm_doppler()
  {
    return n_dopp;
  }

  double TDL_Channel::get_LOS_power()
  {
    return los_power;
  }

  double TDL_Channel::get_LOS_doppler()
  {
    return los_dopp;
  }

  Fading_Generation_Method TDL_Channel::get_generation_method()
  {
    return method;
  }

  double TDL_Channel::calc_mean_excess_delay()
  {
    return ( sqr(a_prof)*d_prof / sum_sqr(a_prof));
  }

  double TDL_Channel::calc_rms_delay_spread()
  {
    double a = ( sqr(a_prof)*d_prof / sum_sqr(a_prof));
    double b = ( sqr(a_prof)*sqr(to_vec(d_prof)) / sum_sqr(a_prof) );

    return ( std::sqrt(b-a*a) );
  }

  void TDL_Channel::init()
  {
    it_assert(N_taps > 0, "TDL_Channel:: init: channel profile not defined yet");
    fading_gen.set_size(N_taps, false);

    if (tap_doppler_spectrum.size() == 0) { // dopplerspectrum is not set. Assume Jakes
      tap_doppler_spectrum.set_size(N_taps);
      for(int i=0; i<N_taps; i++)
	tap_doppler_spectrum(i) = Jakes;
    }

    // create all generators and set the parameters
    switch (method) {
    case FIR:
      for(int i=0; i<N_taps; i++) { fading_gen(i) = new FIR_Fading_Generator(n_dopp, tap_doppler_spectrum(i)); }
      break;

    case IFFT:
      for(int i=0; i<N_taps; i++) { fading_gen(i) = new IFFT_Fading_Generator(n_dopp, tap_doppler_spectrum(i)); }
      break;

    case Rice_MEDS:
      // Ni= smallest number of doppler frequencies, increase by 2 for each tap to make taps uncorrelated
      for(int i=0; i<N_taps; i++) { fading_gen(i) = new Rice_Fading_Generator(n_dopp, tap_doppler_spectrum(i), 16+i*2, MEDS); }
      break;
      
    default:
      it_error("TDL_Channel::init(): Generation method is not implemented");
    };

    if (tap_doppler_spectrum(0) == Rice && los_power > 0.0) { // set LOS component
      fading_gen(0)->set_LOS(los_power, los_dopp);
    }

    // initialize all fading generators
    for(int i=0; i<N_taps; i++)
      fading_gen(i)->init();

    init_flag = true;
  }

  void TDL_Channel::generate(const int no_samples, Array<cvec> &channel_coeff)
  {
    if(init_flag == false)
      init();

    channel_coeff.set_size(N_taps, false);

    for (int i=0; i<N_taps; i++)
      channel_coeff(i) = a_prof(i) * fading_gen(i)->generate(no_samples);
  }

  void TDL_Channel::filter_known_channel(const cvec &input, cvec &output, const Array<cvec> &channel_coeff)
  {
    int maxdelay = max(d_prof);

    output.set_size(input.size()+maxdelay, false);
    output.zeros();

    for (int i=0; i<N_taps; i++)
      output += concat( zeros_c(d_prof(i)), elem_mult(input, channel_coeff(i)), zeros_c(maxdelay-d_prof(i)) );
  }

  void TDL_Channel::filter(const cvec &input, cvec &output, Array<cvec> &channel_coeff)
  {
    generate(input.size(), channel_coeff);
    filter_known_channel(input, output, channel_coeff);
  }

  cvec TDL_Channel::filter(const cvec &input, Array<cvec> &channel_coeff)
  {
    cvec output;
    filter(input, output, channel_coeff);
    return output;
  }

  void TDL_Channel::filter(const cvec &input, cvec &output)
  {
    Array<cvec> channel_coeff;
    filter(input, output, channel_coeff);
  }

  cvec TDL_Channel::filter(const cvec &input)
  {
    cvec output;
    filter(input, output);
    return output;
  }

  cvec TDL_Channel::operator()(const cvec &input, Array<cvec> &channel_coeff)
  {
    return filter(input, channel_coeff);
  }

  cvec TDL_Channel::operator()(const cvec &input)
  {
    return filter(input);
  }



  //------------------------------------------------------------------------------
  // Binary Symetric Channel, BSC.
  //------------------------------------------------------------------------------


  bvec BSC::operator()(const bvec &input)
  {
    int i, length = input.length();
    bvec output(length);
	
    for (i=0; i<length; i++) {
      if (u() <= p) {
	output(i) = input(i) + bin(1);
      } else {
	output(i) = input(i);
      }
    }
    return output;
  }



  cvec AWGN_Channel::operator()(const cvec &input)
  {
    cvec output = input + sigma*randn_c(input.size());
    return output;
  }

  vec AWGN_Channel::operator()(const vec &input)
  {
    vec output = input + sigma*randn(input.size());
    return output;
  }


 } //namespace itpp