refactor (co)sinus calculation out of phase_noise() into openair1/PHY/TOOLS/get_sin_cos.c

CMakeLists.txt

@@ -823,6 +823,7 @@ set(PHY_SRC_COMMON
   ${OPENAIR1_DIR}/PHY/TOOLS/sqrt.c
   ${OPENAIR1_DIR}/PHY/TOOLS/lut.c
   ${OPENAIR1_DIR}/PHY/TOOLS/simde_operations.c
+  ${OPENAIR1_DIR}/PHY/TOOLS/get_sin_cos.c
   )
 
 set(PHY_SRC

openair1/PHY/TOOLS/get_sin_cos.c

@@ -0,0 +1,70 @@
+/*
+ * 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
+ *
+ * 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 <stdint.h>
+#include <math.h>
+#include "tools_defs.h"
+
+#define LUT_SIN_RESOL 128
+// look-up table for the sine (cosine) function
+static uint16_t LUTSin[LUT_SIN_RESOL + 1];
+
+void InitSinLUT(void)
+{
+  for (int i = 0; i < (LUT_SIN_RESOL + 1); i++) {
+    LUTSin[i] = round(sin((M_PI * i) / (2 * LUT_SIN_RESOL)) * (1 << 14)); // Format: Q14
+  }
+}
+
+// fast calculation of e^{i*2*pi*phase}
+// ret.r == cosinus << 14
+// ret.i == sinus << 14
+c16_t get_sin_cos(double phase)
+{
+  int index;
+  if (phase < 0) {
+    index = (int)(-phase * 4 * LUT_SIN_RESOL + 0.5) % (4 * LUT_SIN_RESOL);
+    index = 4 * LUT_SIN_RESOL - index;
+  } else {
+    index = (int)(phase * 4 * LUT_SIN_RESOL + 0.5) % (4 * LUT_SIN_RESOL);
+  }
+
+  c16_t ret;
+  if (index < 2 * LUT_SIN_RESOL) { // check for 1st and 2nd Quadrant
+    if (index < LUT_SIN_RESOL) { // 1st Quadrant
+      ret.r = LUTSin[LUT_SIN_RESOL - index];
+      ret.i = LUTSin[index];
+    } else { // 2nd Quadrant
+      ret.r = -LUTSin[index - LUT_SIN_RESOL];
+      ret.i = LUTSin[2 * LUT_SIN_RESOL - index];
+    }
+  } else { // 3rd and 4th Quadrant
+    if (index < 3 * LUT_SIN_RESOL) { // 3rd Quadrant
+      ret.r = -LUTSin[3 * LUT_SIN_RESOL - index];
+      ret.i = -LUTSin[index - 2 * LUT_SIN_RESOL];
+    } else { // 4th Quadrant
+      ret.r = LUTSin[index - 3 * LUT_SIN_RESOL];
+      ret.i = -LUTSin[4 * LUT_SIN_RESOL - index];
+    }
+  }
+
+  return ret;
+}

openair1/PHY/TOOLS/tools_defs.h

@@ -792,6 +792,13 @@ double interp(double x, double *xs, double *ys, int count);
 void simde_mm128_separate_real_imag_parts(simde__m128i *out_re, simde__m128i *out_im, simde__m128i in0, simde__m128i in1);
 void simde_mm256_separate_real_imag_parts(simde__m256i *out_re, simde__m256i *out_im, simde__m256i in0, simde__m256i in1);
 
+void InitSinLUT(void);
+
+// fast calculation of e^{i*2*pi*phase}
+// ret.r == cosinus << 14
+// ret.i == sinus << 14
+c16_t get_sin_cos(double phase);
+
 #ifdef __cplusplus
 }
 #endif

openair1/SIMULATION/TOOLS/phase_noise.c

@@ -22,49 +22,18 @@
 #include <math.h>
 #include "sim.h"
 
-#define ResolSinCos 100
-static uint16_t LUTSin[ResolSinCos+1];
-
-/* Initialisation function for SIN table values */
-void InitSinLUT( void ) {
-  for ( int i=0; i<(ResolSinCos+1); i++ ) {
-    LUTSin[i] = sin((double)(M_PI*i)/(2*ResolSinCos)) * (1<<14); //Format: Q14
 
 /* linear phase noise model */
-void phase_noise(double ts, int16_t * InRe, int16_t * InIm)
+void phase_noise(double ts, int16_t *InRe, int16_t *InIm)
 {
   static uint64_t i=0;
-  int32_t x=0 ,y=0;
   double fd = 300;//0.01*30000
-  int16_t SinValue = 0, CosValue= 0;
-  double IdxDouble = (double)(i*fd * ts * ResolSinCos * 4);
-  int16_t IdxModulo = ((int32_t)(IdxDouble>0 ? IdxDouble+0.5 : IdxDouble-0.5)) % (ResolSinCos*4);
-  IdxModulo = IdxModulo<0 ? IdxModulo+ResolSinCos*4 : IdxModulo;
-
-  if(IdxModulo<2*ResolSinCos) {//< 2 check for 1st and 2nd
-    if(IdxModulo < ResolSinCos) {// 1st Quadrant
-      SinValue = LUTSin[IdxModulo];
-      CosValue = LUTSin[ResolSinCos-IdxModulo];
-    }
-    else {// 2nd Quadrant
-      SinValue = LUTSin[2*ResolSinCos-IdxModulo];
-      CosValue = -LUTSin[IdxModulo-ResolSinCos];
-    }
-  }
-  else {// 3rd and 4th Quadrant
-    if(IdxModulo < 3*ResolSinCos) {// 3rd Quadrant
-      SinValue = -LUTSin[IdxModulo-2*ResolSinCos];
-      CosValue = -LUTSin[3*ResolSinCos-IdxModulo];
-    }
-    else {//4th Quadrant
-      SinValue = -LUTSin[4*ResolSinCos-IdxModulo];
-      CosValue = LUTSin[IdxModulo-3*ResolSinCos];
-    }
-  }
-  x = ( ((int32_t)InRe[0] * CosValue) - ((int32_t)InIm[0] * SinValue ));
-  y = ( ((int32_t)InIm[0] * CosValue) + ((int32_t)InRe[0] * SinValue ));
-  InRe[0]= (int16_t)(x>>14);
-  InIm[0]= (int16_t)(y>>14);
+  double phase = (double)(i * fd * ts);
+  c16_t rot = get_sin_cos(phase);
+  c16_t val = {.r = *InRe, .i = *InIm};
+  val = c16mulShift(val, rot, 14);
+  *InRe = val.r;
+  *InIm = val.i;
   i++;
 }
 

openair1/SIMULATION/TOOLS/sim.h

@@ -543,8 +543,6 @@ int load_channellist(uint8_t nb_tx, uint8_t nb_rx, double sampling_rate, uint64_
 double N_RB2sampling_rate(uint16_t N_RB);
 double N_RB2channel_bandwidth(uint16_t N_RB);
 
-void InitSinLUT(void);
-
 /* Linear phase noise model */
 /**
   \brief This function produce phase noise and add to input signal