Replace gaussdouble() by gaussZiggurat() in random_channel() and add_noise()

openair1/SIMULATION/TOOLS/multipath_channel.c

@@ -160,8 +160,8 @@ void add_noise(c16_t **rxdata,
   for (int i = 0; i < length; i++) {
     for (int ap = 0; ap < nb_antennas_rx; ap++) {
       c16_t *rxd = &rxdata[ap][slot_offset + i + delay];
-      rxd->r = r_re[ap][i] + sqrt(sigma / 2) * gaussdouble(0.0, 1.0); // convert to fixed point
-      rxd->i = r_im[ap][i] + sqrt(sigma / 2) * gaussdouble(0.0, 1.0);
+      rxd->r = r_re[ap][i] + sqrt(sigma / 2) * gaussZiggurat(0.0, 1.0); // convert to fixed point
+      rxd->i = r_im[ap][i] + sqrt(sigma / 2) * gaussZiggurat(0.0, 1.0);
       /* Add phase noise if enabled */
       if (pdu_bit_map & PUSCH_PDU_BITMAP_PUSCH_PTRS) {
         phase_noise(ts, &rxdata[ap][slot_offset + i + delay].r, &rxdata[ap][slot_offset + i + delay].i);

openair1/SIMULATION/TOOLS/random_channel.c

@@ -467,8 +467,9 @@ double get_normalization_ch_factor(channel_desc_t *desc)
     for (int l = 0; l < (int)desc->nb_taps; l++) {
       for (int aarx = 0; aarx < desc->nb_rx; aarx++) {
         for (int aatx = 0; aatx < desc->nb_tx; aatx++) {
-          anew[aarx + (aatx * desc->nb_rx)].r = sqrt(desc->ricean_factor * desc->amps[l] / 2) * gaussdouble(0.0, 1.0);
-          anew[aarx + (aatx * desc->nb_rx)].i = sqrt(desc->ricean_factor * desc->amps[l] / 2) * gaussdouble(0.0, 1.0);
+          struct complexd *anewp = &anew[aarx + (aatx * desc->nb_rx)];
+          anewp->r = sqrt(desc->ricean_factor * desc->amps[l] / 2) * gaussZiggurat(0.0, 1.0);
+          anewp->i = sqrt(desc->ricean_factor * desc->amps[l] / 2) * gaussZiggurat(0.0, 1.0);
           if ((l == 0) && (desc->ricean_factor != 1.0)) {
             anew[aarx + (aatx * desc->nb_rx)].r += sqrt((1.0 - desc->ricean_factor) / 2);
             anew[aarx + (aatx * desc->nb_rx)].i += sqrt((1.0 - desc->ricean_factor) / 2);
@@ -527,6 +528,15 @@ channel_desc_t *new_channel_desc_scm(uint8_t nb_tx,
                                      int32_t channel_offset,
                                      double path_loss_dB,
                                      float noise_power_dB) {
+
+  // To create tables for normal distribution
+  uint64_t rand;
+  FILE *h = fopen("/dev/random", "r");
+  if (fread(&rand, sizeof(rand), 1, h) != 1) {
+    LOG_W(HW, "Simulator can't read /dev/random\n");
+  }
+  tableNor(rand);
+
   channel_desc_t *chan_desc = (channel_desc_t *)calloc(1,sizeof(channel_desc_t));
 
   for(int i=0; i<max_chan; i++) {
@@ -1712,8 +1722,9 @@ int random_channel(channel_desc_t *desc, uint8_t abstraction_flag) {
     for (aarx=0; aarx<desc->nb_rx; aarx++) {
       for (aatx=0; aatx<desc->nb_tx; aatx++) {
 
-        anew[aarx + (aatx * desc->nb_rx)].r = sqrt(desc->ricean_factor * desc->amps[i] / 2) * gaussdouble(0.0, 1.0) * desc->normalization_ch_factor;
-        anew[aarx + (aatx * desc->nb_rx)].i = sqrt(desc->ricean_factor * desc->amps[i] / 2) * gaussdouble(0.0, 1.0) * desc->normalization_ch_factor;
+        struct complexd *anewp = &anew[aarx + (aatx * desc->nb_rx)];
+        anewp->r = sqrt(desc->ricean_factor * desc->amps[i] / 2) * gaussZiggurat(0.0, 1.0) * desc->normalization_ch_factor;
+        anewp->i = sqrt(desc->ricean_factor * desc->amps[i] / 2) * gaussZiggurat(0.0, 1.0) * desc->normalization_ch_factor;
 
         if ((i==0) && (desc->ricean_factor != 1.0)) {
           if (desc->random_aoa==1) {

openair1/SIMULATION/TOOLS/rangen_double.c

@@ -109,7 +109,78 @@ double gaussdouble(double mean, double variance)
   }
 }
 
+// Ziggurat
+static double wn[128], fn[128];
+static uint32_t iz, jz, jsr = 123456789, kn[128];
+static int32_t hz;
+#define SHR3 (jz = jsr, jsr ^= (jsr << 13), jsr ^= (jsr >> 17), jsr ^= (jsr << 5), jz + jsr)
+#define UNI (0.5 + (signed)SHR3 * 0.2328306e-9)
 
+double nfix(void)
+{
+  const double r = 3.442620;
+  static double x, y;
+
+  for (;;) {
+    x = hz * wn[iz];
+
+    if (iz == 0) {
+      do {
+        x = -0.2904764 * log(UNI);
+        y = -log(UNI);
+      } while (y + y < x * x);
+
+      return (hz > 0) ? r + x : -r - x;
+    }
+
+    if (fn[iz] + UNI * (fn[iz - 1] - fn[iz]) < exp(-0.5 * x * x)) {
+      return x;
+    }
+
+    hz = SHR3;
+    iz = hz & 127;
+
+    if (abs(hz) < kn[iz]) {
+      return ((hz)*wn[iz]);
+    }
+  }
+}
+
+/*!\Procedure to create tables for normal distribution kn,wn and fn. */
+void tableNor(unsigned long seed)
+{
+  jsr = seed;
+  double dn = 3.442619855899;
+  int i;
+  const double m1 = 2147483648.0;
+  double q;
+  double tn = 3.442619855899;
+  const double vn = 9.91256303526217E-03;
+  q = vn / exp(-0.5 * dn * dn);
+  kn[0] = ((dn / q) * m1);
+  kn[1] = 0;
+  wn[0] = (q / m1);
+  wn[127] = (dn / m1);
+  fn[0] = 1.0;
+  fn[127] = (exp(-0.5 * dn * dn));
+
+  for (i = 126; 1 <= i; i--) {
+    dn = sqrt(-2.0 * log(vn / dn + exp(-0.5 * dn * dn)));
+    kn[i + 1] = ((dn / tn) * m1);
+    tn = dn;
+    fn[i] = (exp(-0.5 * dn * dn));
+    wn[i] = (dn / m1);
+  }
+
+  return;
+}
+
+double gaussZiggurat(double mean, double variance)
+{
+  hz = SHR3;
+  iz = hz & 127;
+  return abs(hz) < kn[iz] ? hz * wn[iz] : nfix();
+}
 
 #ifdef MAIN
 main(int argc,char **argv)

openair1/SIMULATION/TOOLS/sim.h

@@ -513,6 +513,8 @@ int gauss(unsigned int *gauss_LUT,unsigned char Nbits);
 double gaussdouble(double,double);
 void randominit(unsigned int seed_init);
 double uniformrandom(void);
+double gaussZiggurat(double mean, double variance);
+void tableNor(unsigned long seed);
 int freq_channel(channel_desc_t *desc,uint16_t nb_rb, int16_t n_samples,int scs);
 int init_freq_channel(channel_desc_t *desc,uint16_t nb_rb,int16_t n_samples,int scs);
 void term_freq_channel(void);