rfsimulator: enable and correct NTN LEO doppler simulation

openair1/SIMULATION/TOOLS/random_channel.c

@@ -1663,7 +1663,7 @@ channel_desc_t *new_channel_desc_scm(uint8_t nb_tx,
       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
+      chan_desc->enable_dynamic_Doppler = true;
       fill_channel_desc(chan_desc,nb_tx,
                         nb_rx,
                         nb_taps,
@@ -1762,6 +1762,11 @@ void set_channeldesc_name(channel_desc_t *cdesc,char *modelname) {
   cdesc->model_name=strdup(modelname);
 }
 
+void set_channeldesc_direction(channel_desc_t *cdesc, bool is_uplink)
+{
+  cdesc->is_uplink = is_uplink;
+}
+
 #ifdef DEBUG_CH_POWER
 double accumulated_ch_power = 0;
 int ch_power_count = 0;

openair1/SIMULATION/TOOLS/sim.h

@@ -139,6 +139,8 @@ typedef struct {
   float Doppler_phase_inc;
   /// current Doppler phase of each RX antenna (for continuous phase from one block to the next)
   float *Doppler_phase_cur;
+  /// flag indicating if channel direction is UL or DL
+  bool is_uplink;
 } channel_desc_t;
 
 typedef struct {
@@ -358,6 +360,7 @@ void free_channel_desc_scm(channel_desc_t *ch);
 \param module_id identifies the channel model. should be define as a macro in simu.h
 */
 void set_channeldesc_owner(channel_desc_t *cdesc, channelmod_moduleid_t module_id);
+
 /**
 \brief This function set a model name to a model descriptor, can be later used to identify a allocated channel model
 \param cdesc points to the model descriptor
@@ -365,6 +368,13 @@ void set_channeldesc_owner(channel_desc_t *cdesc, channelmod_moduleid_t module_i
 */
 void set_channeldesc_name(channel_desc_t *cdesc,char *modelname);
 
+/**
+\brief This function set a channel model direction to either uplink or downlink
+\param cdesc points to the model descriptor
+\param is_uplink indicates if this channel is applied in uplink (not downlink) direction
+*/
+void set_channeldesc_direction(channel_desc_t *cdesc, bool is_uplink);
+
 /** \fn void get_cexp_doppler(struct complexd *cexp_doppler, channel_desc_t *chan_desc, const uint32_t length)
 \brief This routine generates the complex exponential to apply the Doppler shift
 \param cexp_doppler Output with the complex exponential of Doppler shift

radio/rfsimulator/apply_channelmod.c

@@ -50,6 +50,8 @@ void rxAddInput(const c16_t *input_sig,
                 uint32_t CirSize,
                 bool add_noise)
 {
+  static uint64_t last_TS = 0;
+
   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.
@@ -82,42 +84,76 @@ void rxAddInput(const c16_t *input_sig,
     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;
+    if (channelDesc->is_uplink) {
+      const double dir_ue_sat_x = pos_sat_x - pos_ue_x;
+      const double dir_ue_sat_y = pos_sat_y - pos_ue_y;
+      const double dir_ue_sat_z = pos_sat_z - pos_ue_z;
 
-    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);
+      const double dist_ue_sat = sqrt(dir_ue_sat_x * dir_ue_sat_x + dir_ue_sat_y * dir_ue_sat_y + dir_ue_sat_z * dir_ue_sat_z);
+      const double vel_ue_sat = (vel_sat_x * dir_ue_sat_x + vel_sat_y * dir_ue_sat_y + vel_sat_z * dir_ue_sat_z) / dist_ue_sat;
+
+      double dist_sat_gnb = 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_sat_gnb_x = pos_gnb_x - pos_sat_x;
+        const double dir_sat_gnb_y = pos_gnb_y - pos_sat_y;
+        const double dir_sat_gnb_z = pos_gnb_z - pos_sat_z;
+
+        dist_sat_gnb = sqrt(dir_sat_gnb_x * dir_sat_gnb_x + dir_sat_gnb_y * dir_sat_gnb_y + dir_sat_gnb_z * dir_sat_gnb_z);
+      }
+
+      const double prop_delay = (dist_ue_sat + dist_sat_gnb) / c;
+      if (channelDesc->enable_dynamic_delay)
+        channelDesc->channel_offset = prop_delay * channelDesc->sampling_rate;
+
+      const double f_Doppler_shift_ue_sat = (-vel_ue_sat / c) * channelDesc->center_freq;
+      if (channelDesc->enable_dynamic_Doppler)
+        channelDesc->Doppler_phase_inc = 2 * M_PI * f_Doppler_shift_ue_sat / channelDesc->sampling_rate;
+
+      if(TS - last_TS >= channelDesc->sampling_rate) {
+        last_TS = TS;
+        LOG_D(HW, "Satellite orbit: time %f s, Position = (%f, %f, %f), Velocity = (%f, %f, %f)\n", t, pos_sat_x, pos_sat_y, pos_sat_z, vel_sat_x, vel_sat_y, vel_sat_z);
+        LOG_D(HW, "Uplink delay %f ms, Doppler shift UE->SAT %f kHz\n", prop_delay * 1000, f_Doppler_shift_ue_sat / 1000);
+      }
+    } else {
+      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;
+      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);
+      }
+
+      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;
+      if (channelDesc->enable_dynamic_Doppler)
+        channelDesc->Doppler_phase_inc = 2 * M_PI * f_Doppler_shift_sat_ue / channelDesc->sampling_rate;
+
+      if(TS - last_TS >= channelDesc->sampling_rate) {
+        last_TS = TS;
+        LOG_D(HW, "Satellite orbit: time %f s, Position = (%f, %f, %f), Velocity = (%f, %f, %f)\n", t, pos_sat_x, pos_sat_y, pos_sat_z, vel_sat_x, vel_sat_y, vel_sat_z);
+        LOG_D(HW, "Downlink delay %f ms, Doppler shift SAT->UE %f kHz\n", prop_delay * 1000, f_Doppler_shift_sat_ue / 1000);
+      }
     }
   }
 
@@ -131,6 +167,8 @@ void rxAddInput(const c16_t *input_sig,
   const double noise_per_sample = add_noise ? pow(10,channelDesc->noise_power_dB/10.0) * 256 : 0;
   const uint64_t dd = channelDesc->channel_offset;
   const int nbTx=channelDesc->nb_tx;
+  double Doppler_phase_cur = channelDesc->Doppler_phase_cur[rxAnt];
+  Doppler_phase_cur -= 2 * M_PI * round(Doppler_phase_cur / (2 * M_PI));
 
   for (int i=0; i<nbSamples; i++) {
     cf_t *out_ptr = after_channel_sig + i;
@@ -161,10 +199,10 @@ void rxAddInput(const c16_t *input_sig,
 #else
       double complex in = rx_tmp.r + rx_tmp.i * I;
 #endif
-      double complex out = in * cexp(channelDesc->Doppler_phase_cur[rxAnt] * I);
+      double complex out = in * cexp(Doppler_phase_cur * I);
       rx_tmp.r = creal(out);
       rx_tmp.i = cimag(out);
-      channelDesc->Doppler_phase_cur[rxAnt] += channelDesc->Doppler_phase_inc;
+      Doppler_phase_cur += channelDesc->Doppler_phase_inc;
     }
 
     out_ptr->r += rx_tmp.r * pathLossLinear + noise_per_sample * gaussZiggurat(0.0, 1.0);
@@ -172,6 +210,8 @@ void rxAddInput(const c16_t *input_sig,
     out_ptr++;
   }
 
+  channelDesc->Doppler_phase_cur[rxAnt] = Doppler_phase_cur;
+
   if ( (TS*nbTx)%CirSize+nbSamples <= CirSize )
     // Cast to a wrong type for compatibility !
     LOG_D(HW,"Input power %f, output power: %f, channel path loss %f, noise coeff: %f \n",

radio/rfsimulator/simulator.c

@@ -274,6 +274,7 @@ static int allocCirBuf(rfsimulator_state_t *bridge, int sock)
     }
 
     set_channeldesc_owner(ptr->channel_model, RFSIMU_MODULEID);
+    set_channeldesc_direction(ptr->channel_model, bridge->role == SIMU_ROLE_SERVER);
     random_channel(ptr->channel_model,false);
     LOG_I(HW, "Random channel %s in rfsimulator activated\n", modelname);
   }
@@ -467,6 +468,7 @@ static int rfsimu_setchanmod_cmd(char *buff, int debug, telnet_printfunc_t prnt,
                                                           t->chan_pathloss,
                                                           0); // noise_power
           set_channeldesc_owner(newmodel, RFSIMU_MODULEID);
+          set_channeldesc_direction(newmodel, t->role == SIMU_ROLE_SERVER);
           set_channeldesc_name(newmodel,modelname);
           random_channel(newmodel,false);
           channel_desc_t *oldmodel=b->channel_model;