cvec prePost (cvec data_time, int pre, int post){

   data_time.ins(0, data_time.right(pre));
   data_time.ins(data_time.length(), data_time.mid(pre,post));
   return(data_time);

}

ivec Ant1Relay::bpsk_nc_ml_demapping(Array<complex<double> > &miso_out, cvec symbols) {

	ivec res_label;
	res_label.set_size(miso_out.size());

	for(int i=0; i<miso_out.size(); i++) {

		complex<double> Y = miso_out(i);
		cvec X = zeros_c(2);

		double min_norm = numeric_limits<double>::max();
		int min_pt = -1;


		for(int q=0; q<symbols.size(); q++) {
			X(0) = symbols(q);
			for(int r=0; r<symbols.size(); r++) {
				X(1) = symbols(r);
				complex<double> hx = channel * X;
				double calc_norm = norm(Y-hx);
				if(calc_norm < min_norm) {
					min_norm = calc_norm;
					min_pt = (q^r) & 0x01; // xor
				}
			}
		}

		res_label(i) = min_pt;
	}

	return res_label;
}

void roots(const vec &p, cvec &r)
{
  int n = p.size(), m, l;
  ivec f = find(p != 0.0);
  m = f.size();
  vec v = p;
  mat A;

  if (m > 0 && n > 1) {
    v = v(f(0), f(m - 1));
    l = v.size();

    if (l > 1) {

      A = diag(ones(l - 2), -1);
      A.set_row(0, -v(1, l - 1) / v(0));
      r = eig(A);
      cvec d;
      cmat V;
      eig(A, d , V);

      if (f(m - 1) < n)
        r = concat(r, zeros_c(n - f(m - 1) - 1));
    }
    else {
      r.set_size(n - f(m - 1) - 1, false);
      r.zeros();
    }
  }
  else
    r.set_size(0, false);
}

cvec fir_x::process(bvec ce, cvec x)
{
 cvec y;	
 int N;

	#if (DEBUG_LEVEL==3)
	cout << "***** fir_x::proc *****" << endl;	
	cout << "ce=" << ce << endl;
	cout << "x=" << x << endl;
	sleep(1000);
	#endif

	N=ce.length();
	if (x.length()!=N) { 
		throw sci_exception("fir_x::process - ce.size <> x.size", x.length() );
	}	
	y.set_size(N);
	for (int i=0; i<N; i++) {
		if (bool(ce[i])) y0=update(x[i]);				
		y[i]=y0;
	}
	

	#if (DEBUG_LEVEL==3)
	cout << "y=" << y << endl;
	cout << "+++++ fir_x::proc +++++" << endl;	
	sleep(1000);
	#endif

	return (y);
}

static
void assert_cvec(const cvec &ref, const cvec &act)
{
  ASSERT_EQ(ref.length(), act.length());
  for (int n = 0; n < ref.length(); ++n) {
    ASSERT_NEAR(ref[n].real(), act[n].real(), tol);
    ASSERT_NEAR(ref[n].imag(), act[n].imag(), tol);
  }
}

static
void assert_cvec_p(const cvec &ref, const cvec &act, int line)
{
  ASSERT_EQ(ref.length(), act.length()) << line;
  for (int n = 0; n < ref.length(); ++n) {
    ASSERT_NEAR(ref(n).real(), act(n).real(), tol) << line;
    ASSERT_NEAR(ref(n).imag(), act(n).imag(), tol) << line;
  }
}

cvec operator/(const cvec &v, const double &s)
{
  it_assert_debug(v.size() > 0, "operator/(): Vector of zero length");

  cvec temp = v;
  for (int i = 0;i < v.size();i++) {
    temp(i) /= (double)s;
  }
  return temp;
}

void zero_pad_back(cvec& vec_in, int m_in)
{
  int quotient = vec_in.length() / m_in;
  if (quotient*m_in == vec_in.length()) {
    return;
  }
  else {
    int n_remainder = (quotient + 1) * m_in - vec_in.length();
    vec_in = concat(vec_in, zeros_c(n_remainder));
  }
}

	cvec descrambling(cvec data,cvec code)
	{
		cvec result;
		result.set_length(data.length());
		result.zeros();
		for(int i=0;i<data.length();i++)
		{
			result[i]=data(i)*code(i);
		}
		return result;
	}

void poly(const cvec &r, cvec &p)
{
  int n = r.size();

  p.set_size(n + 1, false);
  p.zeros();
  p(0) = 1.0;

  for (int i = 0; i < n; i++)
    p.set_subvector(1, p(1, i + 1) - r(i)*p(0, i));
}

void auto_correlation(const cvec &a, cvec &x)
{
    int k = a.size();
    x.set_size(k);
    for (int i = 0; i < k; ++i) {
        x[i] = 0;
        for (int j = 0; j < k - i; ++j)
            x[i] += a[j] * conj(a[j + i]);
        for (int j = k - i; j < k; ++j)
            x[i] += a[j] * conj(a[j - k + i]);
    }
}

cvec operator*(const cmat &m, const cvec &v)
{
  it_assert_debug(m.no_cols == v.size(), "cmat::operator*(): Wrong sizes");
  cvec r(m.no_rows);
  std::complex<double> alpha = std::complex<double>(1.0);
  std::complex<double> beta = std::complex<double>(0.0);
  char trans = 'n';
  int incr = 1;
  blas::zgemv_(&trans, &m.no_rows, &m.no_cols, &alpha, m.data, &m.no_rows,
               v._data(), &incr, &beta, r._data(), &incr);
  return r;
}

void fir_x::set_state(cvec cv)
{
int k,N;

	N=get_size();
	if (cv.length()!=N ) {
		throw sci_exception("fir::set_state v.lenght <> fir_size", cv.length() );
	}	

	for (k=0; k < N; k++) {
		z[k]=cv[k];
	}
}

void Plot::Plot_data(cvec ydata, int plot_number){

    double x[ydata.length()];
    double y[ydata.length()];
    if(plot_number==2){
        curve2.setStyle(QwtPlotCurve::CurveStyle(4));
        for (int i=0;i<ydata.length();i++){
            x[i]=real(ydata.get(i));
            y[i]=imag(ydata.get(i));
        }
            curve2.setData(x,y,ydata.length());
            gui->qwtPlot_2->replot();
    }
}

cvec polyval(const cvec &p, const vec &x)
{
  it_error_if(p.size() == 0, "polyval: size of polynomial is zero");
  it_error_if(x.size() == 0, "polyval: size of input value vector is zero");

  cvec out(x.size());

  out = p(0);

  for (int i = 1; i < p.size(); i++)
    out = std::complex<double>(p(i)) + elem_mult(to_cvec(x), out);

  return out;
}

/** Return the channel FIR estimation based on trainSeq and fill phiHat, Ahat, delay, sigmaSqrNoise
 *
 * @pre:
 *    - dataC: cvec of the received data
 *    - trainC: pointer to array of the training sequence
 *             !!!  length assume to be < than dataC.length()!!!
 *    - delay: pointer to a int with synchornization point
 *
 * @post:
 *    - Estimators are now filled!
 *    - returned the channel FIR estimation
 */
void synchCatchChannel(cvec dataC, cvec trainCUp , int *delay ){

  if(dataC.length()<trainCUp.length()){
    cerr << "The training sequence is too big compare to datas!\n";
    exit(1);
  }

  //std::cout<<"Data="<<dataC<<"\n";
  //std::cout<<"Train="<<trainCUp<<"\n";
  
  //Computes crosscorrelation:

  int timeLag=1000;

  cvec dataCSmall=dataC(1,timeLag);

  cvec crossCorr=itpp::xcorr(dataCSmall, trainCUp, timeLag);

  
    //std::cout<<"Cross="<<crossCorr<<"\n";

  //Search for the maximum
  double max=0.0;
  int index=0, count=0;
  double tmp[crossCorr.length()-timeLag+1];
  for(int i=timeLag-1;i<crossCorr.length();i++){
    tmp[count]=abs(crossCorr[i]);
   
    //DispVal(tmp[count]);
    //DispVal(count);
    if (tmp[count]> max ){
      max=tmp[count];
      index=i;
     
    }
     count++;
  }
  
  index=index-timeLag+1;
  *delay=index;


  // Save data to file
     std::ofstream ofs( "xcorr.dat" , std::ifstream::out );
     ofs.write((char * ) tmp, (crossCorr.length()-timeLag+1)*sizeof(double));
     ofs.close();

}

void Matlab_Engine::get(cvec &v, const char *name){
  mxArray *T;
  if((T = engGetArray(e, name)) == NULL)
    cout << "No variable with the given name exists in the Matlab workspace!\n";  
  else{
    const int *dims = mxGetDimensions(T);
    const int ndims = mxGetNumberOfDimensions(T);
    if((ndims>2) || !((dims[0]==1)||(dims[1]==1)))
      cout << "The requested variable is not a vector!\n";    
    if(!mxIsComplex(T))
      cout << "The requested variable is real-valued!\n";    
    double *pt_r = (double*) mxGetPr(T);
    double *pt_i = (double*) mxGetPi(T);
    const int N = mxGetNumberOfElements(T);
    v.set_size(N, false);
    if(mxIsComplex(T)) // Copy both real and imaginary part.
      for(int k=0; k<N; k++){
	complex<double> value(*pt_r++, *pt_i++);
	v(k) = value;
      }
    else // Copy only the real part.
      for(int k=0; k<N; k++){
	complex<double> value(*pt_r++, 0);
	v(k) = value;
      }         
  }
}

cvec dataAlloc (cvec data_qam, std::complex<double> data_pilot[], double pilot_pattern[], double data_pattern[], int N){

  int n=0;
  int m=0;
  int u=0;
  std::complex<double> zeroC(0,0);
  cvec data_alloc(N);

  for(int i=0;i<N;i++){;
   if (pilot_pattern[m]-1==i){
     data_alloc.set(i,data_pilot[m]);   
      m++;
     
   }
   else if (data_pattern[n]-1==i){
     data_alloc.set(i,data_qam.get(n));   
      n++;
   }
   else{
     data_alloc.set(i,zeroC);
     u++;
   }
  
 }

  return(data_alloc);

}

cvec xcorr(const cvec &x, const cvec &y, const int max_lag, const std::string scaleopt)
{
  cvec out(2*x.length() - 1); //Initial size does ont matter, it will get adjusted
  xcorr(x, y, out, max_lag, scaleopt, false);

  return out;
}

std::complex<double> operator*(const ivec &a, const cvec &b)
{
  it_assert_debug(a.size() == b.size(), "operator*(): sizes does not match");
  std::complex<double> temp = 0;
  for (int i = 0;i < a.size();i++) {temp += (double)a(i) * b(i);}
  return temp;
}