/*****************************************************************************
 * Author: Christos Komninakis                                               *
 * Date: January 4, 2000.                                                    *
 * Content: C++ code file for the CLASS "flat_rayleigh".                     *
 *****************************************************************************/

#include "flat_rayleigh.h"

#define PI 3.14159265358979

/**************** Constructor ***************************************/
flat_rayleigh::flat_rayleigh(int sed, double fD)
{
  chan_seed = sed; //-23;
  chan_val = Complex(0.0, 0.0);
  K = 7; H = 7; H2 = 2*H;
  I = (int) (0.2 / fD);
  last_i = 0; IP = 0;
  int i, t, j, k;
  
  a = new double [K];
  b = new double [K];
  c = new double [K];
  d = new double [K];
  double A_o;
  st = new Complex* [3];
  for (i=0; i<3; i++)
    st[i] = new Complex [K];
  sinc_matrix = new double* [I];
  for (i=0; i<I; i++)
    sinc_matrix[i] = new double [H2];
  buff_f = new Complex [H2];

  /* Follwing are the values of a designed ROOT DOPPLER filter prototype,
     consisting of K=7 biquads, with fD = 0.2 */

  a[0] = 0.98323218491754;    b[0] = 0.13953807910136;
  c[0] = -0.60134822489663;   d[0] = 0.48938674026073;
  a[1] = -0.58056307877395;   b[1] = 0.99977164701498;
  c[1] = -0.42846077703117;   d[1] = 0.24612477470991;
  a[2] = -0.50498097209321;   b[2] = 0.99986531836969;
  c[2] = -0.63050157648166;   d[2] = 0.9350218659307;
  a[3] = 0.20534422769486;    b[3] = 0.98691526068222;
  c[3] = -0.096351029308223;  d[3] = -0.060828833103004;
  a[4] = 1.2667005604552;     b[4] = 0.93699896920708;
  c[4] = -0.64949215100022;   d[4] = 0.80427856161307;
  a[5] = -0.60545728027338;   b[5] = 0.99998818741754;
  c[5] = -0.61939958040555;   d[5] = 0.99895172368533;
  a[6] = -0.29077222501287;   b[6] = 0.99832970528904;
  c[6] = -0.62124137243814;   d[6] = 0.98111826315356;
  A_o = 0.028241759549181;

  /* now really start */
  Ao = A_o / sqrt(2.0);
  /* produce the table storing all the values of the interpolating function,
     windowed and in reversed order, so that filtering is straightforward */
  double temp;
  for (i=0; i<I; i++)
    for (j=0; j<H2; j++)
      sinc_matrix[i][j] = sin(PI*(j-H)+(PI/I)*i) / (PI*(j-H)+(PI/I)*i) *
                          (0.5 - 0.5*cos(2.0*PI*(i+I*j)/(H2*I))); // windowed, Hanning
                          //(0.54 - 0.46*cos(2.0*PI*(i+I*j)/(H2*I))); // windowed, Hamming
                          //(0.42 - 0.5*cos(2.0*PI*(i+I*j)/(H2*I)) + 0.08*cos(4.0*PI*(i+I*j)/(H2*I))); // windowed, Blackman
  sinc_matrix[0][H] = 1.0;
  for (i=0; i<I; i++)
    for (j=0; j<H; j++) {
      temp = sinc_matrix[i][j];
      sinc_matrix[i][j] = sinc_matrix[i][H2-1-j];
      sinc_matrix[i][H2-1-j] = temp;
    }

  //ofstream out("channel.txt");
  /* work the channel for some time, to overcome transients */
  for (t=0; t<I*2*(4*K+H2); t++) {
    if (last_i == I-1) {
      last_i = 0;

      /* now produce new element, and insert it into the buffer: buff_f[IP] */
      st[0][0] = Complex(Gaussian(&chan_seed), Gaussian(&chan_seed))*Ao -
                 st[1][0]*c[0] - st[2][0]*d[0];
      for (k=1; k<K; k++)
        st[0][k] = st[0][k-1] + st[1][k-1]*a[k-1] + st[2][k-1]*b[k-1] -
                                st[1][k]*c[k] - st[2][k]*d[k];
      buff_f[IP] = st[0][K-1] + st[1][K-1]*a[K-1] + st[2][K-1]*b[K-1];
      for (k=0; k<K; k++) {
        st[2][k] = st[1][k];
        st[1][k] = st[0][k];
      }
      /* done, produced buff_f[IP], now move IP to next position */
      IP = (IP + 1) % H2;
    } else
      last_i++;

    buff_sinc = sinc_matrix[last_i];
    chan_val = Complex(0.0, 0.0);
    for (j=0; j<H2; j++)
      chan_val += buff_f[(IP + j) % H2] * buff_sinc[j];
      //out << chan_val << endl;
  }
  //out.close;
}



/******* Function pass_through, without providing any CSI ************/
void flat_rayleigh::pass_through(int length, Complex *inp, Complex *outp)
{
  register int k, j, t;

  
  for (t=0; t<length; t++) {
    if (last_i == I-1) {
      last_i = 0;

      /* now produce new element, and insert it into the buffer: buff_f[IP] */
      st[0][0] = Complex( Gaussian(&chan_seed), Gaussian(&chan_seed) )*Ao -
                 st[1][0]*c[0] - st[2][0]*d[0];
      for (k=1; k<K; k++)
        st[0][k] = st[0][k-1] + st[1][k-1]*a[k-1] + st[2][k-1]*b[k-1] -
                                st[1][k]*c[k] - st[2][k]*d[k];
      buff_f[IP] = st[0][K-1] + st[1][K-1]*a[K-1] + st[2][K-1]*b[K-1];
      for (k=0; k<K; k++) {
        st[2][k] = st[1][k];
        st[1][k] = st[0][k];
      }
      /* done, produced buff_f[IP], now move IP to next position */
      IP = (IP + 1) % H2;
    } else
      last_i++;

    buff_sinc = sinc_matrix[last_i];
    chan_val = Complex(0.0, 0.0);
    for (j=0; j<H2; j++)
      chan_val += buff_f[(IP + j) % H2] * buff_sinc[j];
    outp[t] = inp[t] * chan_val;
  }
}



/******* Overloaded Function pass_through, providing CSI also ************/
void flat_rayleigh::pass_through(int length, Complex *inp,
                                    Complex *outp, Complex *csi)
{
  register int k, j, t;

  
  for (t=0; t<length; t++) {
    if (last_i == I-1) {
      last_i = 0;

      /* now produce new element, and insert it into the buffer: buff_f[IP] */
      st[0][0] = Complex( Gaussian(&chan_seed), Gaussian(&chan_seed) )*Ao -
                 st[1][0]*c[0] - st[2][0]*d[0];
      for (k=1; k<K; k++)
        st[0][k] = st[0][k-1] + st[1][k-1]*a[k-1] + st[2][k-1]*b[k-1] -
                                st[1][k]*c[k] - st[2][k]*d[k];
      buff_f[IP] = st[0][K-1] + st[1][K-1]*a[K-1] + st[2][K-1]*b[K-1];
      for (k=0; k<K; k++) {
        st[2][k] = st[1][k];
        st[1][k] = st[0][k];
      }
      /* done, produced buff_f[IP], now move IP to next position */
      IP = (IP + 1) % H2;
    } else
      last_i++;

    buff_sinc = sinc_matrix[last_i];
    chan_val = Complex(0.0, 0.0);
    for (j=0; j<H2; j++)
      chan_val += buff_f[(IP + j) % H2] * buff_sinc[j];
    csi[t] = chan_val;
    outp[t] = inp[t] * chan_val;
  }
}



/***** Overloaded Function pass_through, providing AMPLITUDE CSI only *****/
void flat_rayleigh::pass_through(int length, Complex *inp,
                                 Complex *outp, double *amp_csi)
{
  register int k, j, t;
  
  for (t=0; t<length; t++) {
    if (last_i == I-1) {
      last_i = 0;

      /* now produce new element, and insert it into the buffer: buff_f[IP] */
      st[0][0] = Complex( Gaussian(&chan_seed), Gaussian(&chan_seed) )*Ao -
                 st[1][0]*c[0] - st[2][0]*d[0];
      for (k=1; k<K; k++)
        st[0][k] = st[0][k-1] + st[1][k-1]*a[k-1] + st[2][k-1]*b[k-1] -
                                st[1][k]*c[k] - st[2][k]*d[k];
      buff_f[IP] = st[0][K-1] + st[1][K-1]*a[K-1] + st[2][K-1]*b[K-1];
      for (k=0; k<K; k++) {
        st[2][k] = st[1][k];
        st[1][k] = st[0][k];
      }
      /* done, produced buff_f[IP], now move IP to next position */
      IP = (IP + 1) % H2;
    } else
      last_i++;

    buff_sinc = sinc_matrix[last_i];
    chan_val = Complex(0.0, 0.0);
    for (j=0; j<H2; j++)
      chan_val += buff_f[(IP + j) % H2] * buff_sinc[j];
    amp_csi[t] = sqrt(abs(chan_val));
    outp[t] = inp[t] * chan_val;
  }
}