From e44a99b3c1724d9f28ba4c5cd8014d7a5b511e45 Mon Sep 17 00:00:00 2001
From: Roberto Louro Magueta <rmagueta@allbesmart.pt>
Date: Wed, 9 Nov 2022 15:06:01 +0000
Subject: [PATCH] Implementation of get_cexp_doppler() function to compute the
 Doppler Effect

---
 openair1/SIMULATION/TOOLS/random_channel.c | 38 ++++++++++++++++++++++
 openair1/SIMULATION/TOOLS/sim.h            |  2 ++
 2 files changed, 40 insertions(+)

diff --git a/openair1/SIMULATION/TOOLS/random_channel.c b/openair1/SIMULATION/TOOLS/random_channel.c
index 39c5e5f564..8774108437 100644
--- a/openair1/SIMULATION/TOOLS/random_channel.c
+++ b/openair1/SIMULATION/TOOLS/random_channel.c
@@ -24,6 +24,7 @@
 #include <stdio.h>
 #include <stdlib.h>
 #include <string.h>
+#include <complex.h>
 
 
 #include "PHY/TOOLS/tools_defs.h"
@@ -37,6 +38,7 @@
 
 //#define DEBUG_CH
 //#define DEBUG_CH_POWER
+//#define DOPPLER_DEBUG
 
 #include "assertions.h"
 
@@ -451,6 +453,42 @@ void tdlModel(int  tdl_paths, double *tdl_delays, double *tdl_amps_dB, double DS
   }
 }
 
+void get_cexp_doppler(struct complexd *cexp_doppler, channel_desc_t *chan_desc)
+{
+  // TS 38.104 - Table G.3-1
+  uint16_t Dmin = 2;
+  uint16_t Ds = 300;
+  double c = 299792458;
+  double v = chan_desc->max_Doppler * (c / (double)chan_desc->center_freq);
+
+#ifdef DOPPLER_DEBUG
+  printf("v = %f\n", v);
+#endif
+
+  double phase0 = 2 * M_PI * uniformrandom();
+  double cos_theta[chan_desc->channel_length];
+  double fs[chan_desc->channel_length];
+  for (int t_idx = 0; t_idx < chan_desc->channel_length; t_idx++) {
+    double t = t_idx / (chan_desc->sampling_rate * 1e6);
+    if (t >= 0 && t <= Ds / v) {
+      cos_theta[t_idx] = (Ds / 2 - v * t) / sqrt(Dmin * Dmin + (Ds / 2 - v * t) * (Ds / 2 - v * t));
+    } else if (t > Ds / v && t <= 2 * Ds / v) {
+      cos_theta[t_idx] = (-1.5 * Ds + v * t) / sqrt(Dmin * Dmin + (-1.5 * Ds + v * t) * (-1.5 * Ds + v * t));
+    } else {
+      cos_theta[t_idx] = cos(fmod(t, 2 * Ds / v));
+    }
+    fs[t_idx] = chan_desc->max_Doppler * cos_theta[t_idx];
+
+    double complex tmp_cexp_doppler = cexp(I * (2 * M_PI * fs[t_idx] * t + phase0));
+    cexp_doppler[t_idx].r = creal(tmp_cexp_doppler);
+    cexp_doppler[t_idx].i = cimag(tmp_cexp_doppler);
+
+#ifdef DOPPLER_DEBUG
+    printf("(%2i) t_us = %f, cos_theta = %f, fs = %f, cexp_doppler = (%f, %f)\n", t_idx, t * 1e6, cos_theta[t_idx], fs[t_idx], cexp_doppler[t_idx].r, cexp_doppler[t_idx].i);
+#endif
+  }
+}
+
 double get_normalization_ch_factor(channel_desc_t *desc)
 {
   if (!(desc->channel_length > 1 && desc->modelid >= TDL_A && desc->modelid <= TDL_E)) {
diff --git a/openair1/SIMULATION/TOOLS/sim.h b/openair1/SIMULATION/TOOLS/sim.h
index 1e2602650f..0704bae02a 100644
--- a/openair1/SIMULATION/TOOLS/sim.h
+++ b/openair1/SIMULATION/TOOLS/sim.h
@@ -80,6 +80,8 @@ typedef struct {
   struct complexd **chF;
   ///Maximum path delay in mus.
   double Td;
+  ///Carrier center frequency
+  uint64_t center_freq;
   ///Channel bandwidth in MHz.
   double channel_bandwidth;
   ///System sampling rate in Msps.