openairinterface5g/openair1/SIMULATION/TOOLS/multipath_tv_channel.c

/*******************************************************************************
    OpenAirInterface 
    Copyright(c) 1999 - 2014 Eurecom

    OpenAirInterface 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.


    OpenAirInterface 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 OpenAirInterface.The full GNU General Public License is 
   included in this distribution in the file called "COPYING". If not, 
   see <http://www.gnu.org/licenses/>.

  Contact Information
  OpenAirInterface Admin: openair_admin@eurecom.fr
  OpenAirInterface Tech : openair_tech@eurecom.fr
  OpenAirInterface Dev  : openair4g-devel@eurecom.fr
  
  Address      : Eurecom, Compus SophiaTech 450, route des chappes, 06451 Biot, France.

 *******************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include "defs.h"
#include "SIMULATION/RF/defs.h"
#include <complex.h>

void tv_channel(channel_desc_t *desc,double complex ***H,uint16_t length);
double frand_a_b(double a, double b);
void tv_conv(double complex **h, double complex *x, double complex *y, uint16_t nb_samples, uint8_t nb_taps, int delay);

void multipath_tv_channel(channel_desc_t *desc,
			  double **tx_sig_re, 
			  double **tx_sig_im, 
			  double **rx_sig_re,
			  double **rx_sig_im,
			  uint16_t length,
			  uint8_t keep_channel) {
  
    double complex **tx,**rx,***H_t,*rx_temp;//, *tv_H_t;
  double path_loss = pow(10,desc->path_loss_dB/20);
  int i,j,k,dd;

  dd = abs(desc->channel_offset);

#ifdef DEBUG_CH
  printf("[TV CHANNEL] keep = %d : path_loss = %g (%f), nb_rx %d, nb_tx %d, dd %d, len %d max_doppler %d\n",keep_channel,path_loss,desc->path_loss_dB,desc->nb_rx,desc->nb_tx,dd,desc->channel_length,desc->max_Doppler);
#endif

  tx = (double complex **)malloc(desc->nb_tx*sizeof(double complex));
  rx = (double complex **)malloc(desc->nb_rx*sizeof(double complex));
  
  H_t= (double complex ***) malloc(desc->nb_tx*desc->nb_rx*sizeof(double complex **));
  //  tv_H_t = (double complex *) malloc(length*sizeof(double complex));
  rx_temp= (double complex *) calloc(length,sizeof(double complex));
  
  for(i=0;i<desc->nb_tx;i++){
      tx[i] = (double complex *)calloc(length,sizeof(double complex));
  }
  
  for(i=0;i<desc->nb_rx;i++){
      rx[i] = (double complex *)calloc(length,sizeof(double complex));
  }
  
  for(i=0;i<desc->nb_tx*desc->nb_rx;i++){
      H_t[i] = (double complex **) malloc(length*sizeof(double complex *));
      for(j=0;j<length;j++) {
          H_t[i][j] = (double complex *) calloc (desc->nb_taps,sizeof(double complex));
      }
  }
  
  for (i=0;i<desc->nb_tx;i++){
      for(j=0;j<length;j++) {
          tx[i][j] = tx_sig_re[i][j] +I*tx_sig_im[i][j];   
      }
  }
  
  if (keep_channel) {
      // do nothing - keep channel
  } else {
      tv_channel(desc,H_t,length);
  }  
   
  for(i=0;i<desc->nb_rx;i++){
      for(j=0;j<desc->nb_tx;j++){
          tv_conv(H_t[i+(j*desc->nb_rx)],tx[j],rx_temp,length,desc->nb_taps,dd);
          for(k=0;k<length;k++){
              rx[i][k] += rx_temp[k];
          }
      }	
  }

  for(i=0;i<desc->nb_rx;i++){
      for(j=0;j<length;j++){
          rx_sig_re[i][j] = creal(rx[i][j])*path_loss;
          rx_sig_im[i][j] = cimag(rx[i][j])*path_loss;
      }
  }

  /*  for(k=0;k<length;k++) {
      tv_H_t[k] = H_t[0][k][0];
      }*/

  for(i=0;i<desc->nb_tx;i++){
      free(tx[i]);
  }
  free(tx);
  for(i=0;i<desc->nb_rx;i++) {
      free(rx[i]);
  }
  free(rx);
    
  for(i=0;i<desc->nb_rx*desc->nb_tx;i++){
      for(j=0;j<length;j++){
          free(H_t[i][j]);
      }
      free(H_t[i]);
  }
  free(H_t);
  
  free(rx_temp);
}

//TODO: make phi_rad a parameter of this function
void tv_channel(channel_desc_t *desc,double complex ***H,uint16_t length){
  
  int i,j,p,l,k;
  double *alpha,*phi_rad,pi=acos(-1),*w_Hz;
  
  alpha = (double *)calloc(desc->nb_paths,sizeof(double));
  phi_rad = (double *)calloc(desc->nb_paths,sizeof(double));
  w_Hz = (double *)calloc(desc->nb_paths,sizeof(double));
  
  for(i=0;i<desc->nb_paths;i++) {
      w_Hz[i]=desc->max_Doppler*cos(frand_a_b(0,2*pi));
      phi_rad[i]=frand_a_b(0,2*pi);
  }
  
  if(desc->ricean_factor == 1){
      for(i=0;i<desc->nb_paths;i++) {
          alpha[i]=1/sqrt(desc->nb_paths);
      }
  }
  else {
      alpha[0]=sqrt(desc->ricean_factor/(desc->ricean_factor+1));
      for(i=1;i<desc->nb_paths;i++) {
          alpha[i] = (1/sqrt(desc->nb_paths-1))*(sqrt(1/(desc->ricean_factor+1)));
      }
  }
  /*
  // This is the code when we only consider a SISO case
  for(i=0;i<length;i++)
  {
	for(j=0;j<desc->nb_taps;j++)
	   {
		for(p=0;p<desc->nb_paths;p++)
		   {
		     H[i][j] += sqrt(desc->amps[j]/2)*alpha[p]*cexp(-I*(2*pi*w_Hz[p]*i*(1/(desc->BW*1e6))+phi_rad[p]));
		   }
   	   }
   }
for(j=0;j<desc->nb_paths;j++)
   {
	phi_rad[j] = fmod(2*pi*w_Hz[j]*(length-1)*(1/desc->BW)+phi_rad[j],2*pi);
   }
  */

  // if MIMO
  for (i=0;i<desc->nb_rx;i++){
      for(j=0;j<desc->nb_tx;j++){
          for(k=0;k<length;k++){
              for(l=0;l<desc->nb_taps;l++){
                  H[i+(j*desc->nb_rx)][k][l] = 0;
                  for(p=0;p<desc->nb_paths;p++){
                      H[i+(j*desc->nb_rx)][k][l] += sqrt(desc->amps[l]/2)*alpha[p]*cexp(I*(2*pi*w_Hz[p]*k*(1/(desc->BW*1e6))+phi_rad[p]));
                  }
              }
          }
          for(j=0;j<desc->nb_paths;j++){
              phi_rad[j] = fmod(2*pi*w_Hz[j]*(length-1)*(1/desc->BW)+phi_rad[j],2*pi);
          }
      }
  }

  free(alpha);
  free(w_Hz);
  free(phi_rad);
}

// time varying convolution 
void tv_conv(double complex **h, double complex *x, double complex *y, uint16_t nb_samples, uint8_t nb_taps, int dd) {
  
  int i,j;
  for(i=0;i<((int)nb_samples-dd);i++) {
      for(j=0;j<nb_taps;j++) {
          if(i>j)
              y[i+dd] += creal(h[i][j])*creal(x[i-j])-cimag(h[i][j])*cimag(x[i-j]) + I*(creal(h[i][j])*cimag(x[i-j])+cimag(h[i][j])*creal(x[i-j]));
          
      }
  }
}

double frand_a_b(double a, double b){
    return (rand()/(double)RAND_MAX)*(b-a)+a;
}