multipath_tv_channel.c 6.46 KB
Newer Older
ghaddab's avatar
ghaddab committed
1
/*******************************************************************************
2
    OpenAirInterface
ghaddab's avatar
ghaddab committed
3 4 5 6 7 8 9 10 11 12 13 14 15 16
    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
17 18
    along with OpenAirInterface.The full GNU General Public License is
   included in this distribution in the file called "COPYING". If not,
ghaddab's avatar
ghaddab committed
19 20 21 22 23
   see <http://www.gnu.org/licenses/>.

  Contact Information
  OpenAirInterface Admin: openair_admin@eurecom.fr
  OpenAirInterface Tech : openair_tech@eurecom.fr
24
  OpenAirInterface Dev  : openair4g-devel@lists.eurecom.fr
25

ghaddab's avatar
ghaddab committed
26
  Address      : Eurecom, Campus SophiaTech, 450 Route des Chappes, CS 50193 - 06904 Biot Sophia Antipolis cedex, FRANCE
ghaddab's avatar
ghaddab committed
27 28

 *******************************************************************************/
29 30 31 32 33 34 35 36
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include "defs.h"
#include "SIMULATION/RF/defs.h"
#include <complex.h>

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

void multipath_tv_channel(channel_desc_t *desc,
42 43 44 45 46 47 48 49 50
                          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;
51 52 53 54 55 56
  double path_loss = pow(10,desc->path_loss_dB/20);
  int i,j,k,dd;

  dd = abs(desc->channel_offset);

#ifdef DEBUG_CH
57 58
  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);
59 60 61 62
#endif

  tx = (double complex **)malloc(desc->nb_tx*sizeof(double complex));
  rx = (double complex **)malloc(desc->nb_rx*sizeof(double complex));
63

64 65 66
  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));
67 68 69

  for(i=0; i<desc->nb_tx; i++) {
    tx[i] = (double complex *)calloc(length,sizeof(double complex));
70
  }
71 72 73

  for(i=0; i<desc->nb_rx; i++) {
    rx[i] = (double complex *)calloc(length,sizeof(double complex));
74
  }
75 76 77 78 79 80 81

  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));
    }
82
  }
83 84 85 86 87

  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];
    }
88
  }
89

90
  if (keep_channel) {
91
    // do nothing - keep channel
92
  } else {
93
    tv_channel(desc,H_t,length);
94 95
  }

96 97 98 99 100 101
  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];
102
      }
103 104 105 106 107 108 109 110
    }
  }

  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;
    }
111 112 113 114 115 116
  }

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

117 118
  for(i=0; i<desc->nb_tx; i++) {
    free(tx[i]);
119
  }
120

121
  free(tx);
122 123 124

  for(i=0; i<desc->nb_rx; i++) {
    free(rx[i]);
125
  }
126

127
  free(rx);
128 129 130 131 132 133 134

  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]);
135
  }
136

137
  free(H_t);
138

139 140 141 142
  free(rx_temp);
}

//TODO: make phi_rad a parameter of this function
143 144 145
void tv_channel(channel_desc_t *desc,double complex ***H,uint16_t length)
{

146 147
  int i,j,p,l,k;
  double *alpha,*phi_rad,pi=acos(-1),*w_Hz;
148

149 150 151
  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));
152 153 154 155

  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);
156
  }
157 158 159 160 161 162 163 164 165 166 167

  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)));
    }
168
  }
169

170 171 172 173
  /*
  // This is the code when we only consider a SISO case
  for(i=0;i<length;i++)
  {
174 175 176 177
  for(j=0;j<desc->nb_taps;j++)
     {
    for(p=0;p<desc->nb_paths;p++)
       {
Florian Kaltenberger's avatar
Florian Kaltenberger committed
178
         H[i][j] += sqrt(desc->amps[j]/2)*alpha[p]*cexp(-I*(2*pi*w_Hz[p]*i*(1/(desc->sampling_rate*1e6))+phi_rad[p]));
179 180
       }
       }
181
   }
182
  for(j=0;j<desc->nb_paths;j++)
183
   {
Florian Kaltenberger's avatar
Florian Kaltenberger committed
184
  phi_rad[j] = fmod(2*pi*w_Hz[j]*(length-1)*(1/desc->sampling_rate)+phi_rad[j],2*pi);
185 186 187 188
   }
  */

  // if MIMO
189 190 191 192 193 194 195
  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++) {
Florian Kaltenberger's avatar
Florian Kaltenberger committed
196
            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->sampling_rate*1e6))+phi_rad[p]));
197
          }
198 199 200 201
        }
      }

      for(j=0; j<desc->nb_paths; j++) {
Florian Kaltenberger's avatar
Florian Kaltenberger committed
202
        phi_rad[j] = fmod(2*pi*w_Hz[j]*(length-1)*(1/desc->sampling_rate)+phi_rad[j],2*pi);
203
      }
204
    }
205 206 207 208 209 210 211
  }

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

212 213 214 215
// 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)
{

216
  int i,j;
217 218 219 220 221 222 223

  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]));

    }
224 225 226
  }
}

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