multipath_tv_channel.c 6.22 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21
/*
 * Licensed to the OpenAirInterface (OAI) Software Alliance under one or more
 * contributor license agreements.  See the NOTICE file distributed with
 * this work for additional information regarding copyright ownership.
 * The OpenAirInterface Software Alliance licenses this file to You under
 * the OAI Public License, Version 1.0  (the "License"); you may not use this file
 * except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *      http://www.openairinterface.org/?page_id=698
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 *-------------------------------------------------------------------------------
 * For more information about the OpenAirInterface (OAI) Software Alliance:
 *      contact@openairinterface.org
 */

22 23 24 25 26 27 28 29
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include "defs.h"
#include "SIMULATION/RF/defs.h"
#include <complex.h>

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

void multipath_tv_channel(channel_desc_t *desc,
35 36 37 38 39 40 41 42 43
                          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;
44 45 46 47 48 49
  double path_loss = pow(10,desc->path_loss_dB/20);
  int i,j,k,dd;

  dd = abs(desc->channel_offset);

#ifdef DEBUG_CH
50 51
  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);
52 53 54 55
#endif

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

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

  for(i=0; i<desc->nb_tx; i++) {
    tx[i] = (double complex *)calloc(length,sizeof(double complex));
63
  }
64 65 66

  for(i=0; i<desc->nb_rx; i++) {
    rx[i] = (double complex *)calloc(length,sizeof(double complex));
67
  }
68 69 70 71 72 73 74

  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));
    }
75
  }
76 77 78 79 80

  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];
    }
81
  }
82

83
  if (keep_channel) {
84
    // do nothing - keep channel
85
  } else {
86
    tv_channel(desc,H_t,length);
87 88
  }

89 90 91 92 93 94
  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];
95
      }
96 97 98 99 100 101 102 103
    }
  }

  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;
    }
104 105 106 107 108 109
  }

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

110 111
  for(i=0; i<desc->nb_tx; i++) {
    free(tx[i]);
112
  }
113

114
  free(tx);
115 116 117

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

120
  free(rx);
121 122 123 124 125 126 127

  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]);
128
  }
129

130
  free(H_t);
131

132 133 134 135
  free(rx_temp);
}

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

139 140
  int i,j,p,l,k;
  double *alpha,*phi_rad,pi=acos(-1),*w_Hz;
141

142 143 144
  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));
145 146 147 148

  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);
149
  }
150 151 152 153 154 155 156 157 158 159 160

  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)));
    }
161
  }
162

163 164 165 166
  /*
  // This is the code when we only consider a SISO case
  for(i=0;i<length;i++)
  {
167 168 169 170
  for(j=0;j<desc->nb_taps;j++)
     {
    for(p=0;p<desc->nb_paths;p++)
       {
Florian Kaltenberger's avatar
Florian Kaltenberger committed
171
         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]));
172 173
       }
       }
174
   }
175
  for(j=0;j<desc->nb_paths;j++)
176
   {
Florian Kaltenberger's avatar
Florian Kaltenberger committed
177
  phi_rad[j] = fmod(2*pi*w_Hz[j]*(length-1)*(1/desc->sampling_rate)+phi_rad[j],2*pi);
178 179 180 181
   }
  */

  // if MIMO
182 183 184 185 186 187 188
  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
189
            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]));
190
          }
191 192 193 194
        }
      }

      for(j=0; j<desc->nb_paths; j++) {
Florian Kaltenberger's avatar
Florian Kaltenberger committed
195
        phi_rad[j] = fmod(2*pi*w_Hz[j]*(length-1)*(1/desc->sampling_rate)+phi_rad[j],2*pi);
196
      }
197
    }
198 199 200 201 202 203 204
  }

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

205 206 207 208
// 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)
{

209
  int i,j;
210 211 212 213 214 215 216

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

    }
217 218 219
  }
}

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