From 255ad2441053df7fc88ad1b44a42372fc346f796 Mon Sep 17 00:00:00 2001
From: Roberto Louro Magueta <rmagueta@allbesmart.pt>
Date: Wed, 2 Nov 2022 12:30:22 +0000
Subject: [PATCH] Replace gaussdouble() by gaussZiggurat() in random_channel()
and add_noise()
---
CMakeLists.txt | 1 -
openair1/SIMULATION/TOOLS/multipath_channel.c | 4 +-
openair1/SIMULATION/TOOLS/random_channel.c | 19 ++-
openair1/SIMULATION/TOOLS/rangen_double.c | 71 ++++++++++++
openair1/SIMULATION/TOOLS/sim.h | 2 +
radio/rfsimulator/new_channel_sim.c | 108 ------------------
radio/rfsimulator/rfsimulator.h | 2 -
7 files changed, 90 insertions(+), 117 deletions(-)
delete mode 100644 radio/rfsimulator/new_channel_sim.c
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 8559320d514..4bb1fa6545a 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -2546,7 +2546,6 @@ target_link_libraries(nrscope ${XFORMS_LIBRARIES})
add_library(rfsimulator MODULE
${OPENAIR_DIR}/radio/rfsimulator/simulator.c
${OPENAIR_DIR}/radio/rfsimulator/apply_channelmod.c
- ${OPENAIR_DIR}/radio/rfsimulator/new_channel_sim.c
${OPENAIR1_DIR}/PHY/TOOLS/signal_energy.c
)
target_link_libraries(rfsimulator SIMU_COMMON ${ATLAS_LIBRARIES})
diff --git a/openair1/SIMULATION/TOOLS/multipath_channel.c b/openair1/SIMULATION/TOOLS/multipath_channel.c
index 6adc33b6462..c1bc2689168 100644
--- a/openair1/SIMULATION/TOOLS/multipath_channel.c
+++ b/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);
diff --git a/openair1/SIMULATION/TOOLS/random_channel.c b/openair1/SIMULATION/TOOLS/random_channel.c
index 8f7cc8ea0e5..b11e4b7f675 100644
--- a/openair1/SIMULATION/TOOLS/random_channel.c
+++ b/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) {
diff --git a/openair1/SIMULATION/TOOLS/rangen_double.c b/openair1/SIMULATION/TOOLS/rangen_double.c
index 0ba9e196b13..20bd9c1ecb0 100644
--- a/openair1/SIMULATION/TOOLS/rangen_double.c
+++ b/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)
diff --git a/openair1/SIMULATION/TOOLS/sim.h b/openair1/SIMULATION/TOOLS/sim.h
index 9b3fd98ae62..1e2602650ff 100644
--- a/openair1/SIMULATION/TOOLS/sim.h
+++ b/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);
diff --git a/radio/rfsimulator/new_channel_sim.c b/radio/rfsimulator/new_channel_sim.c
deleted file mode 100644
index fcfa85265b2..00000000000
--- a/radio/rfsimulator/new_channel_sim.c
+++ /dev/null
@@ -1,108 +0,0 @@
-/*
-* 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.1 (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
-*
-* Author and copyright: Laurent Thomas, open-cells.com
-*
-* 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
-*/
-
-
-#include <stdlib.h>
-#include <stdio.h>
-#include <string.h>
-#include <unistd.h>
-#include <stdbool.h>
-#include <errno.h>
-
-
-#include <common/utils/assertions.h>
-#include <common/utils/LOG/log.h>
-#include <common/config/config_userapi.h>
-#include <openair1/SIMULATION/TOOLS/sim.h>
-#include "rfsimulator.h"
-
-// 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();
-}
diff --git a/radio/rfsimulator/rfsimulator.h b/radio/rfsimulator/rfsimulator.h
index e69d41fa6ce..068250a149f 100644
--- a/radio/rfsimulator/rfsimulator.h
+++ b/radio/rfsimulator/rfsimulator.h
@@ -23,8 +23,6 @@
#ifndef __RFSIMULATOR_H
#define __RFSIMULATOR_H
-double gaussZiggurat(double mean, double variance);
-void tableNor(unsigned long seed);
void rxAddInput( const c16_t *input_sig,
c16_t *after_channel_sig,
int rxAnt,
--
GitLab