rfsimulator: add delay and Doppler simulation for circular orbit and add according channel models SAT_LEO_TRANS and SAT_LEO_REGEN

Delay and Doppler computation is done as in Matlab function dopplerShiftCircularOrbit:
https://de.mathworks.com/help/satcom/ref/dopplershiftcircularorbit.html

openair1/SIMULATION/TOOLS/random_channel.c

@@ -128,6 +128,7 @@ void fill_channel_desc(channel_desc_t *chan_desc,
   chan_desc->first_run                  = 1;
   chan_desc->ip                         = 0.0;
   chan_desc->max_Doppler                = max_Doppler;
+  chan_desc->Doppler_phase_cur          = calloc(nb_rx, sizeof(double));
   chan_desc->ch                         = calloc(nb_tx*nb_rx, sizeof(struct complexd *));
   chan_desc->chF                        = calloc(nb_tx*nb_rx, sizeof(struct complexd *));
   chan_desc->a                          = calloc(nb_taps, sizeof(struct complexd *));
@@ -1653,6 +1654,37 @@ channel_desc_t *new_channel_desc_scm(uint8_t nb_tx,
                         0);
       break;
 
+    case SAT_LEO_TRANS:
+    case SAT_LEO_REGEN:
+      nb_taps = 1;
+      Td = 0;
+      channel_length = 1;
+      ricean_factor = 0.0;
+      aoa = 0.0;
+      maxDoppler = 0;
+      chan_desc->sat_height = 600e3;
+      chan_desc->enable_dynamic_delay = true;
+      chan_desc->enable_dynamic_Doppler = false; // TODO: requires UE to support continuous Doppler estimation, compensation and pre-compensation
+      fill_channel_desc(chan_desc,nb_tx,
+                        nb_rx,
+                        nb_taps,
+                        channel_length,
+                        default_amp_lin,
+                        NULL,
+                        NULL,
+                        Td,
+                        sampling_rate,
+                        channel_bandwidth,
+                        ricean_factor,
+                        aoa,
+                        forgetting_factor,
+                        maxDoppler,
+                        channel_offset,
+                        path_loss_dB,
+                        0);
+      printf("%s: satellite orbit height %f km\n", map_int_to_str(channelmod_names, channel_model), chan_desc->sat_height / 1000);
+      break;
+
     default:
       LOG_W(OCM,"channel model not yet supported\n");
       free(chan_desc);
@@ -1711,6 +1743,7 @@ void free_channel_desc_scm(channel_desc_t *ch) {
       free(ch->R_sqrt[i]);
 
   free(ch->R_sqrt);
+  free(ch->Doppler_phase_cur);
   free(ch->ch);
   free(ch->chF);
   free(ch->a);
@@ -1743,9 +1776,9 @@ int random_channel(channel_desc_t *desc, uint8_t abstraction_flag) {
   struct complexd phase, alpha, beta;
   start_meas(&desc->random_channel);
 
-  // For AWGN channel, the received signal (Srx) is equal to transmitted signal (Stx) plus noise (N), i.e., Srx = Stx + N,
+  // For AWGN and SAT_LEO_* channels, the received signal (Srx) is equal to transmitted signal (Stx) plus noise (N), i.e., Srx = Stx + N,
   //  therefore, the channel matrix is the identity matrix.
-  if (desc->modelid == AWGN) {
+  if (desc->modelid == AWGN || desc->modelid == SAT_LEO_TRANS || desc->modelid == SAT_LEO_REGEN) {
     for (aarx=0; aarx<desc->nb_rx; aarx++) {
       for (aatx = 0; aatx < desc->nb_tx; aatx++) {
         desc->ch[aarx+(aatx*desc->nb_rx)][0].r = aarx%desc->nb_tx == aatx ? 1.0 : 0.0;
@@ -2071,6 +2104,8 @@ static void display_channelmodel(channel_desc_t *cd,int debug, telnet_printfunc_
        cd->forgetting_factor);
   prnt("Initial phase: %lf   nb_path: %i \n",
        cd->ip, cd->nb_paths);
+  if (cd->modelid == SAT_LEO_TRANS || cd->modelid == SAT_LEO_REGEN)
+    prnt("satellite orbit height: %f\n", cd->sat_height);
 
   for (int i=0; i<cd->nb_taps ; i++) {
     prnt("taps: %i   lin. ampli. : %lf    delay: %lf \n",i,cd->amps[i], cd->delays[i]);

openair1/SIMULATION/TOOLS/sim.h

@@ -126,6 +126,18 @@ typedef struct {
   char *model_name;  
   /// flags to properly trigger memory free
   unsigned int free_flags;
+  /// time stamp when the time varying channel emulation starts (when client connected)
+  uint64_t start_TS;
+  /// height of LEO satellite
+  float sat_height;
+  /// flag to enable dynamic delay simulation for LEO satellite
+  bool enable_dynamic_delay;
+  /// flag to enable dynamic Doppler simulation for LEO satellite
+  bool enable_dynamic_Doppler;
+  /// Doppler phase increment (might vary over time, e.g. for LEO satellite)
+  float Doppler_phase_inc;
+  /// current Doppler phase of each RX antenna (for continuous phase from one block to the next)
+  float *Doppler_phase_cur;
 } channel_desc_t;
 
 typedef struct {
@@ -218,6 +230,8 @@ typedef enum {
   EPA_low,
   EPA_medium,
   EPA_high,
+  SAT_LEO_TRANS,
+  SAT_LEO_REGEN,
 } SCM_t;
 #define CHANNELMOD_MAP_INIT \
   {"custom",custom},\
@@ -253,6 +267,8 @@ typedef enum {
   {"EPA_low",EPA_low},\
   {"EPA_medium",EPA_medium},\
   {"EPA_high",EPA_high},\
+  {"SAT_LEO_TRANS",SAT_LEO_TRANS},\
+  {"SAT_LEO_REGEN",SAT_LEO_REGEN},\
   {NULL, -1}
 
 #define CONFIG_HLP_SNR     "Set average SNR in dB (for --siml1 option)\n"

radio/rfsimulator/apply_channelmod.c

@@ -28,7 +28,7 @@
 #include <unistd.h>
 #include <stdbool.h>
 #include <errno.h>
-
+#include <complex.h>
 
 #include <common/utils/assertions.h>
 #include <common/utils/LOG/log.h>
@@ -52,14 +52,85 @@
   either we regenerate the channel (call again random_channel(desc,0)), or we keep it over subframes
   legacy: we regenerate each sub frame in UL, and each frame only in DL
 */
-void rxAddInput( const c16_t *input_sig,
+void rxAddInput(const c16_t *input_sig,
-		 c16_t *after_channel_sig,
+                c16_t *after_channel_sig,
-                 int rxAnt,
+                int rxAnt,
-                 channel_desc_t *channelDesc,
+                channel_desc_t *channelDesc,
-                 int nbSamples,
+                int nbSamples,
-                 uint64_t TS,
+                uint64_t TS,
-                 uint32_t CirSize
+                uint32_t CirSize)
-               ) {
+{
+  if ((channelDesc->sat_height > 0) && (channelDesc->enable_dynamic_delay || channelDesc->enable_dynamic_Doppler)) { // model for transparent satellite on circular orbit
+    /* assumptions:
+       - The Earth is spherical, the ground station is static, and that the Earth does not rotate.
+       - An access or link is possible from the satellite to the ground station at all times.
+       - The ground station is located at the North Pole (positive Zaxis), and the satellite starts from the initial elevation angle 0° in the second quadrant of the YZplane.
+       - Satellite moves in the clockwise direction in its circular orbit.
+    */
+    const double radius_earth = 6371e3; // m
+    const double radius_sat = radius_earth + channelDesc->sat_height;
+    const double GM_earth = 3.986e14; // m^3/s^2
+    const double w_sat = sqrt(GM_earth / (radius_sat * radius_sat * radius_sat)); // rad/s
+
+    // start_time and end_time are when the pos_sat_z == pos_gnb_z (elevation angle == 0 and 180 degree)
+    const double start_phase = -acos(radius_earth / radius_sat); // SAT is just rising above the horizon
+    const double end_phase = -start_phase; // SAT is just falling behind the horizon
+    const double start_time = start_phase / w_sat; // in seconds
+    const double end_time = end_phase / w_sat; // in seconds
+    const uint64_t duration_samples = (end_time - start_time) * channelDesc->sampling_rate;
+    const double t = start_time + ((TS - channelDesc->start_TS) % duration_samples) / channelDesc->sampling_rate;
+
+    const double pos_sat_x = 0;
+    const double pos_sat_y = radius_sat * sin(w_sat * t);
+    const double pos_sat_z = radius_sat * cos(w_sat * t);
+
+    const double vel_sat_x = 0;
+    const double vel_sat_y =  w_sat * radius_sat * cos(w_sat * t);
+    const double vel_sat_z = -w_sat * radius_sat * sin(w_sat * t);
+
+    const double pos_ue_x = 0;
+    const double pos_ue_y = 0;
+    const double pos_ue_z = radius_earth;
+
+    const double dir_sat_ue_x = pos_ue_x - pos_sat_x;
+    const double dir_sat_ue_y = pos_ue_y - pos_sat_y;
+    const double dir_sat_ue_z = pos_ue_z - pos_sat_z;
+
+    const double dist_sat_ue = sqrt(dir_sat_ue_x * dir_sat_ue_x + dir_sat_ue_y * dir_sat_ue_y + dir_sat_ue_z * dir_sat_ue_z);
+    const double vel_sat_ue = (vel_sat_x * dir_sat_ue_x + vel_sat_y * dir_sat_ue_y + vel_sat_z * dir_sat_ue_z) / dist_sat_ue;
+
+    double dist_gnb_sat = 0;
+    double vel_gnb_sat = 0;
+    if (channelDesc->modelid == SAT_LEO_TRANS) {
+      const double pos_gnb_x = 0;
+      const double pos_gnb_y = 0;
+      const double pos_gnb_z = radius_earth;
+
+      const double dir_gnb_sat_x = pos_sat_x - pos_gnb_x;
+      const double dir_gnb_sat_y = pos_sat_y - pos_gnb_y;
+      const double dir_gnb_sat_z = pos_sat_z - pos_gnb_z;
+
+      dist_gnb_sat = sqrt(dir_gnb_sat_x * dir_gnb_sat_x + dir_gnb_sat_y * dir_gnb_sat_y + dir_gnb_sat_z * dir_gnb_sat_z);
+      vel_gnb_sat = (vel_sat_x * dir_gnb_sat_x + vel_sat_y * dir_gnb_sat_y + vel_sat_z * dir_gnb_sat_z) / dist_gnb_sat;
+    }
+
+    const double c = 299792458; // m/s
+    const double prop_delay = (dist_gnb_sat + dist_sat_ue) / c;
+    if (channelDesc->enable_dynamic_delay)
+      channelDesc->channel_offset = prop_delay * channelDesc->sampling_rate;
+
+    const double f_Doppler_shift_sat_ue = (vel_sat_ue / (c - vel_sat_ue)) * channelDesc->center_freq;
+    const double f_Doppler_shift_gnb_sat = (-vel_gnb_sat / c) * channelDesc->center_freq;
+    if (channelDesc->enable_dynamic_Doppler)
+      channelDesc->Doppler_phase_inc = 2 * M_PI * (f_Doppler_shift_gnb_sat + f_Doppler_shift_sat_ue) / channelDesc->sampling_rate;
+
+    static uint64_t last_TS = 0;
+    if(TS - last_TS >= channelDesc->sampling_rate) {
+      last_TS = TS;
+      LOG_I(HW, "Satellite orbit: time %f s, Doppler: gNB->SAT %f kHz, SAT->UE %f kHz, Delay %f ms\n", t, f_Doppler_shift_gnb_sat / 1000, f_Doppler_shift_sat_ue / 1000, prop_delay * 1000);
+    }
+  }
+
   // channelDesc->path_loss_dB should contain the total path gain
   // so, in actual RF: tx gain + path loss + rx gain (+antenna gain, ...)
   // UE and NB gain control to be added
@@ -94,9 +165,16 @@ void rxAddInput( const c16_t *input_sig,
       } //l
     }
 
-    // Fixme: lround(), rount(), ... is detected by valgrind as error, not found why
+    if (channelDesc->Doppler_phase_inc != 0.0) {
-    out_ptr->r += lround(rx_tmp.r*pathLossLinear + noise_per_sample*gaussZiggurat(0.0,1.0));
+      double complex in = CMPLX(rx_tmp.r, rx_tmp.i);
-    out_ptr->i += lround(rx_tmp.i*pathLossLinear + noise_per_sample*gaussZiggurat(0.0,1.0));
+      double complex out = in * cexp(channelDesc->Doppler_phase_cur[rxAnt] * I);
+      rx_tmp.r = creal(out);
+      rx_tmp.i = cimag(out);
+      channelDesc->Doppler_phase_cur[rxAnt] += channelDesc->Doppler_phase_inc;
+    }
+
+    out_ptr->r = lround(rx_tmp.r*pathLossLinear + noise_per_sample*gaussZiggurat(0.0,1.0));
+    out_ptr->i = lround(rx_tmp.i*pathLossLinear + noise_per_sample*gaussZiggurat(0.0,1.0));
     out_ptr++;
   }
 

radio/rfsimulator/simulator.c

@@ -709,6 +709,10 @@ static bool flushInput(rfsimulator_state_t *t, int timeout, int nsamps_for_initi
         samplesVoid[i]=(void *)&v;
 
       rfsimulator_write_internal(t, t->lastWroteTS > 1 ? t->lastWroteTS - 1 : 0, samplesVoid, 1, t->tx_num_channels, 1, false);
+
+      buffer_t *b = &t->buf[conn_sock];
+      if (b->channel_model)
+        b->channel_model->start_TS = t->lastWroteTS;
     } else {
       if ( events[nbEv].events & (EPOLLHUP | EPOLLERR | EPOLLRDHUP) ) {
         socketError(t,fd);
@@ -761,6 +765,8 @@ static bool flushInput(rfsimulator_state_t *t, int timeout, int nsamps_for_initi
                             nsamps_for_initial;
           LOG_D(HW, "UE got first timestamp: starting at %lu\n", t->nextRxTstamp);
           b->trashingPacket=true;
+          if (b->channel_model)
+            b->channel_model->start_TS = t->nextRxTstamp;
         } else if (b->lastReceivedTS < b->th.timestamp) {
           int nbAnt= b->th.nbAnt;