Channel convolution function to repalce multipath_channel

openair1/PHY/TOOLS/tests/CMakeLists.txt

@@ -40,6 +40,18 @@ add_dependencies(tests test_sse_intrinsics)
 add_test(NAME test_sse_intrinsics
   COMMAND ./test_sse_intrinsics)
 
+add_executable(test_channel_convolution test_channel_convolution.cpp)
+target_link_libraries(test_channel_convolution PRIVATE GTest::gtest LOG minimal_lib SIMU shlib_loader)
+add_dependencies(tests test_channel_convolution)
+add_test(NAME test_channel_convolution
+  COMMAND ./test_channel_convolution)
+
+add_executable(benchmark_channel_convolution benchmark_channel_convolution.cpp)
+target_link_libraries(benchmark_channel_convolution PRIVATE benchmark::benchmark minimal_lib SIMU shlib_loader)
+add_dependencies(tests benchmark_channel_convolution)
+add_test(NAME benchmark_channel_convolution
+  COMMAND ./benchmark_channel_convolution)
+
 if(CUDA_ENABLE)
 add_executable(test_multipath test_multipath.c)
 target_link_libraries(test_multipath PRIVATE UTIL SIMU LOG CONFIG_LIB shlib_loader m oai_cuda_lib)

openair1/PHY/TOOLS/tests/benchmark_channel_convolution.cpp

@@ -0,0 +1,117 @@
+/*
+ * 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 <vector>
+#include <algorithm>
+#include <numeric>
+extern "C" {
+#include "openair1/PHY/TOOLS/tools_defs.h"
+#include "common/config/config_userapi.h"
+struct configmodule_interface_s;
+configmodule_interface_t *uniqCfg;
+void exit_function(const char *file, const char *function, const int line, const char *s, const int assert)
+{
+  if (assert) {
+    abort();
+  } else {
+    exit(EXIT_SUCCESS);
+  }
+}
+#include "SIMULATION/TOOLS/sim.h"
+}
+#include <cstdio>
+#include "common/utils/LOG/log.h"
+#include "benchmark/benchmark.h"
+#include "openair1/PHY/TOOLS/phy_test_tools.hpp"
+
+static void BM_channel_convolution_template(benchmark::State &state, void (*func)(channel_desc_t *, c16_t **, c16_t **, uint32_t))
+{
+  const int nb_tx = 2;
+  const int nb_rx = 2;
+  const uint32_t length = state.range(0);
+  const uint32_t channel_length = 5;
+  const uint64_t channel_offset = 10;
+  const double path_loss_dB = -3.0;
+
+  // Setup channel descriptor
+  channel_desc_t desc;
+  memset(&desc, 0, sizeof(channel_desc_t));
+  desc.nb_tx = nb_tx;
+  desc.nb_rx = nb_rx;
+  desc.channel_length = channel_length;
+  desc.channel_offset = channel_offset;
+  desc.path_loss_dB = path_loss_dB;
+
+  // Allocate channel taps
+  std::vector<struct complexd *> ch_ptrs(nb_tx * nb_rx);
+  std::vector<std::vector<struct complexd>> ch_data(nb_tx * nb_rx, std::vector<struct complexd>(channel_length));
+
+  desc.ch = ch_ptrs.data();
+  for (int i = 0; i < nb_tx * nb_rx; i++) {
+    desc.ch[i] = ch_data[i].data();
+    for (uint32_t l = 0; l < channel_length; l++) {
+      desc.ch[i][l].r = (double)rand() / RAND_MAX;
+      desc.ch[i][l].i = (double)rand() / RAND_MAX;
+    }
+  }
+
+  // Input
+  std::vector<c16_t *> input_ptrs(nb_tx);
+  std::vector<AlignedVector512<c16_t>> input_data(nb_tx);
+  for (int i = 0; i < nb_tx; i++) {
+    input_data[i] = generate_random_c16(length + 8);
+    input_ptrs[i] = input_data[i].data();
+  }
+
+  // Output
+  std::vector<c16_t *> output_ptrs(nb_rx);
+  std::vector<AlignedVector512<c16_t>> output_data(nb_rx);
+  for (int i = 0; i < nb_rx; i++) {
+    output_data[i].resize(length);
+    output_ptrs[i] = output_data[i].data();
+  }
+
+  for (auto _ : state) {
+    func(&desc, input_ptrs.data(), output_ptrs.data(), length);
+  }
+}
+
+static void BM_channel_convolution(benchmark::State &state)
+{
+  BM_channel_convolution_template(state, channel_convolution);
+}
+
+static void BM_channel_convolution_avx2(benchmark::State &state)
+{
+  BM_channel_convolution_template(state, channel_convolution_avx2);
+}
+
+static void BM_channel_convolution_avx512(benchmark::State &state)
+{
+  BM_channel_convolution_template(state, channel_convolution_avx512);
+}
+
+BENCHMARK(BM_channel_convolution)->RangeMultiplier(4)->Range(128, 30000);
+BENCHMARK(BM_channel_convolution_avx2)->RangeMultiplier(4)->Range(128, 30000);
+BENCHMARK(BM_channel_convolution_avx512)->RangeMultiplier(4)->Range(128, 30000);
+
+BENCHMARK_MAIN();

openair1/PHY/TOOLS/tests/test_channel_convolution.cpp

@@ -0,0 +1,245 @@
+#include <gtest/gtest.h>
+extern "C" {
+#include "openair1/PHY/TOOLS/tools_defs.h"
+#include "common/utils/LOG/log.h"
+#include "SIMULATION/TOOLS/sim.h"
+}
+#include <cmath>
+#include <cstdlib>
+#include <vector>
+
+void channel_convolution_ref(channel_desc_t *desc, c16_t **input, c16_t **output, uint32_t length)
+{
+  float path_loss = powf(10.0f, (float)desc->path_loss_dB / 20.0f);
+  uint64_t dd = desc->channel_offset;
+
+  for (int i = 0; i < ((int)length - (int)dd); i++) {
+    for (int ii = 0; ii < desc->nb_rx; ii++) {
+      cf_t rx_tmp = {0};
+      for (int j = 0; j < desc->nb_tx; j++) {
+        struct complexd *chan = desc->ch[ii + (j * desc->nb_rx)];
+        for (int l = 0; l < (int)desc->channel_length; l++) {
+          if ((i >= 0) && (i - l) >= 0) {
+            cf_t tx;
+            tx.r = input[j][i - l].r;
+            tx.i = input[j][i - l].i;
+            rx_tmp.r += (tx.r * (float)chan[l].r) - (tx.i * (float)chan[l].i);
+            rx_tmp.i += (tx.i * (float)chan[l].r) + (tx.r * (float)chan[l].i);
+          }
+        } // l
+      } // j
+
+      output[ii][i + dd].r = rx_tmp.r * path_loss;
+      output[ii][i + dd].i = rx_tmp.i * path_loss;
+    } // ii
+  } // i
+}
+
+
+TEST(convolve, compare_against_reference)
+{
+  // Parameters
+  const int nb_tx = 2;
+  const int nb_rx = 2;
+  const uint32_t length = 1000;
+  const uint32_t channel_length = 5;
+  const uint64_t channel_offset = 10;
+  const double path_loss_dB = -3.0;
+
+  // Setup channel descriptor
+  channel_desc_t desc;
+  desc.nb_tx = nb_tx;
+  desc.nb_rx = nb_rx;
+  desc.channel_length = channel_length;
+  desc.channel_offset = channel_offset;
+  desc.path_loss_dB = path_loss_dB;
+
+  // Allocate channel taps
+  std::vector<struct complexd *> ch_ptrs(nb_tx * nb_rx);
+  std::vector<std::vector<struct complexd>> ch_data(nb_tx * nb_rx, std::vector<struct complexd>(channel_length));
+
+  desc.ch = ch_ptrs.data();
+  for (int i = 0; i < nb_tx * nb_rx; i++) {
+    desc.ch[i] = ch_data[i].data();
+    for (uint32_t l = 0; l < channel_length; l++) {
+      desc.ch[i][l].r = (double)rand() / RAND_MAX;
+      desc.ch[i][l].i = (double)rand() / RAND_MAX;
+    }
+  }
+
+  // Input
+  std::vector<c16_t *> input_ptrs(nb_tx);
+  std::vector<std::vector<c16_t>> input_data(nb_tx, std::vector<c16_t>(length));
+  for (int i = 0; i < nb_tx; i++) {
+    input_ptrs[i] = input_data[i].data();
+    for (uint32_t j = 0; j < length; j++) {
+      input_data[i][j].r = rand() % 1000;
+      input_data[i][j].i = rand() % 1000;
+    }
+  }
+
+  // Output
+  std::vector<c16_t *> output_ref_ptrs(nb_rx);
+  std::vector<std::vector<c16_t>> output_ref_data(nb_rx, std::vector<c16_t>(length, {0, 0}));
+
+  std::vector<c16_t *> output_test_ptrs(nb_rx);
+  std::vector<std::vector<c16_t>> output_test_data(nb_rx, std::vector<c16_t>(length, {0, 0}));
+
+  for (int i = 0; i < nb_rx; i++) {
+    output_ref_ptrs[i] = output_ref_data[i].data();
+    output_test_ptrs[i] = output_test_data[i].data();
+  }
+
+  // Run
+  channel_convolution_ref(&desc, input_ptrs.data(), output_ref_ptrs.data(), length);
+  channel_convolution(&desc, input_ptrs.data(), output_test_ptrs.data(), length);
+
+  // Verify
+  for (int i = 0; i < nb_rx; i++) {
+    for (uint32_t j = 0; j < length; j++) {
+      EXPECT_NEAR(output_ref_data[i][j].r, output_test_data[i][j].r, 1) << "Mismatch at rx " << i << " sample " << j;
+      EXPECT_NEAR(output_ref_data[i][j].i, output_test_data[i][j].i, 1) << "Mismatch at rx " << i << " sample " << j;
+    }
+  }
+}
+
+TEST(convolve, compare_simd_against_reference)
+{
+  // Parameters
+  const int nb_tx = 2;
+  const int nb_rx = 2;
+  const uint32_t length = 1000;
+  const uint32_t channel_length = 5;
+  const uint64_t channel_offset = 10;
+  const double path_loss_dB = -3.0;
+
+  // Setup channel descriptor
+  channel_desc_t desc;
+  desc.nb_tx = nb_tx;
+  desc.nb_rx = nb_rx;
+  desc.channel_length = channel_length;
+  desc.channel_offset = channel_offset;
+  desc.path_loss_dB = path_loss_dB;
+
+  // Allocate channel taps
+  std::vector<struct complexd *> ch_ptrs(nb_tx * nb_rx);
+  std::vector<std::vector<struct complexd>> ch_data(nb_tx * nb_rx, std::vector<struct complexd>(channel_length));
+
+  desc.ch = ch_ptrs.data();
+  for (int i = 0; i < nb_tx * nb_rx; i++) {
+    desc.ch[i] = ch_data[i].data();
+    for (uint32_t l = 0; l < channel_length; l++) {
+      desc.ch[i][l].r = (double)rand() / RAND_MAX;
+      desc.ch[i][l].i = (double)rand() / RAND_MAX;
+    }
+  }
+
+  // Input - allocate extra padding for SIMD reads
+  std::vector<c16_t *> input_ptrs(nb_tx);
+  std::vector<std::vector<c16_t>> input_data(nb_tx, std::vector<c16_t>(length + 8));
+  for (int i = 0; i < nb_tx; i++) {
+    input_ptrs[i] = input_data[i].data();
+    for (uint32_t j = 0; j < length; j++) {
+      input_data[i][j].r = rand() % 1000;
+      input_data[i][j].i = rand() % 1000;
+    }
+  }
+
+  // Output
+  std::vector<c16_t *> output_ref_ptrs(nb_rx);
+  std::vector<std::vector<c16_t>> output_ref_data(nb_rx, std::vector<c16_t>(length, {0, 0}));
+
+  std::vector<c16_t *> output_test_ptrs(nb_rx);
+  std::vector<std::vector<c16_t>> output_test_data(nb_rx, std::vector<c16_t>(length, {0, 0}));
+
+  for (int i = 0; i < nb_rx; i++) {
+    output_ref_ptrs[i] = output_ref_data[i].data();
+    output_test_ptrs[i] = output_test_data[i].data();
+  }
+
+  // Run
+  channel_convolution_ref(&desc, input_ptrs.data(), output_ref_ptrs.data(), length);
+  channel_convolution_avx2(&desc, input_ptrs.data(), output_test_ptrs.data(), length);
+
+  // Verify
+  for (int i = 0; i < nb_rx; i++) {
+    for (uint32_t j = 0; j < length; j++) {
+      EXPECT_NEAR(output_ref_data[i][j].r, output_test_data[i][j].r, 1) << "Mismatch at rx " << i << " sample " << j;
+      EXPECT_NEAR(output_ref_data[i][j].i, output_test_data[i][j].i, 1) << "Mismatch at rx " << i << " sample " << j;
+    }
+  }
+}
+
+TEST(convolve, compare_avx512_against_reference)
+{
+  // Parameters
+  const int nb_tx = 2;
+  const int nb_rx = 2;
+  const uint32_t length = 1000;
+  const uint32_t channel_length = 5;
+  const uint64_t channel_offset = 10;
+  const double path_loss_dB = -3.0;
+
+  // Setup channel descriptor
+  channel_desc_t desc;
+  desc.nb_tx = nb_tx;
+  desc.nb_rx = nb_rx;
+  desc.channel_length = channel_length;
+  desc.channel_offset = channel_offset;
+  desc.path_loss_dB = path_loss_dB;
+
+  // Allocate channel taps
+  std::vector<struct complexd *> ch_ptrs(nb_tx * nb_rx);
+  std::vector<std::vector<struct complexd>> ch_data(nb_tx * nb_rx, std::vector<struct complexd>(channel_length));
+
+  desc.ch = ch_ptrs.data();
+  for (int i = 0; i < nb_tx * nb_rx; i++) {
+    desc.ch[i] = ch_data[i].data();
+    for (uint32_t l = 0; l < channel_length; l++) {
+      desc.ch[i][l].r = (double)rand() / RAND_MAX;
+      desc.ch[i][l].i = (double)rand() / RAND_MAX;
+    }
+  }
+
+  // Input - allocate extra padding for SIMD reads
+  std::vector<c16_t *> input_ptrs(nb_tx);
+  std::vector<std::vector<c16_t>> input_data(nb_tx, std::vector<c16_t>(length + 8));
+  for (int i = 0; i < nb_tx; i++) {
+    input_ptrs[i] = input_data[i].data();
+    for (uint32_t j = 0; j < length; j++) {
+      input_data[i][j].r = rand() % 1000;
+      input_data[i][j].i = rand() % 1000;
+    }
+  }
+
+  // Output
+  std::vector<c16_t *> output_ref_ptrs(nb_rx);
+  std::vector<std::vector<c16_t>> output_ref_data(nb_rx, std::vector<c16_t>(length, {0, 0}));
+
+  std::vector<c16_t *> output_test_ptrs(nb_rx);
+  std::vector<std::vector<c16_t>> output_test_data(nb_rx, std::vector<c16_t>(length, {0, 0}));
+
+  for (int i = 0; i < nb_rx; i++) {
+    output_ref_ptrs[i] = output_ref_data[i].data();
+    output_test_ptrs[i] = output_test_data[i].data();
+  }
+
+  // Run
+  channel_convolution_ref(&desc, input_ptrs.data(), output_ref_ptrs.data(), length);
+  channel_convolution_avx512(&desc, input_ptrs.data(), output_test_ptrs.data(), length);
+
+  // Verify
+  for (int i = 0; i < nb_rx; i++) {
+    for (uint32_t j = 0; j < length; j++) {
+      EXPECT_NEAR(output_ref_data[i][j].r, output_test_data[i][j].r, 1) << "Mismatch at rx " << i << " sample " << j;
+      EXPECT_NEAR(output_ref_data[i][j].i, output_test_data[i][j].i, 1) << "Mismatch at rx " << i << " sample " << j;
+    }
+  }
+}
+
+int main(int argc, char **argv)
+{
+  logInit();
+  testing::InitGoogleTest(&argc, argv);
+  return RUN_ALL_TESTS();
+}

openair1/SIMULATION/TOOLS/multipath_channel.c

@@ -22,6 +22,8 @@
 #include <math.h>
 #include "PHY/TOOLS/tools_defs.h"
 #include "sim.h"
+#include <simde/x86/avx2.h>
+#include <simde/x86/avx512.h>
 
 //#define DEBUG_CH
 //#define DOPPLER_DEBUG
@@ -390,3 +392,270 @@ void multipath_channel_float(channel_desc_t *desc,
 }
 
 #endif
+
+void channel_convolution(channel_desc_t *desc, c16_t **input, c16_t **output, uint32_t length)
+{
+  float pathloss = powf(10.0f, (float)desc->path_loss_dB / 20.0f);
+  uint L = desc->channel_length;
+  uint64_t offset = desc->channel_offset;
+
+  // Reverse channel and bake in pathloss
+  cf_t rev_channel[desc->nb_rx][desc->nb_tx][L];
+  for (int rx = 0; rx < desc->nb_rx; rx++) {
+    for (int tx = 0; tx < desc->nb_tx; tx++) {
+      struct complexd *chan = desc->ch[rx + (tx * desc->nb_rx)];
+      for (uint i = 0U; i < L; i++) {
+        rev_channel[rx][tx][i].r = (float)chan[L - 1 - i].r * pathloss;
+        rev_channel[rx][tx][i].i = (float)chan[L - 1 - i].i * pathloss;
+      }
+    }
+  }
+
+  for (int rx = 0; rx < desc->nb_rx; rx++) {
+    for (uint i = 0U; i < length - offset; i++) {
+      cf_t rx_sum = {0, 0};
+
+      for (int tx = 0; tx < desc->nb_tx; tx++) {
+        // We look at tx_sig[i - offset] and go 'backwards' for L_eff samples.
+        // To make it a forward-crawl for speed, we start at the oldest
+        // sample in the window and move forward.
+
+        int oldest_idx = i - (L - 1);
+        int current_L = L;
+        int start_c = 0;
+
+        // Handle the "Past" boundary (nothing before index 0 exists)
+        if (oldest_idx < 0) {
+          start_c = -oldest_idx; // Skip the taps that point to negative time
+          current_L -= start_c;
+          oldest_idx = 0;
+        }
+
+        if (current_L > 0) {
+          c16_t *tx_ptr = &input[tx][oldest_idx];
+          cf_t *c_ptr = &rev_channel[rx][tx][start_c];
+
+          // Steady State inner loop: no branches, purely linear
+          for (int l = 0; l < current_L; l++) {
+            rx_sum.r += (tx_ptr->r * c_ptr->r) - (tx_ptr->i * c_ptr->i);
+            rx_sum.i += (tx_ptr->i * c_ptr->r) + (tx_ptr->r * c_ptr->i);
+            tx_ptr++;
+            c_ptr++;
+          }
+        }
+      }
+      output[rx][i + offset].r = rx_sum.r;
+      output[rx][i + offset].i = rx_sum.i;
+    }
+  }
+}
+
+void channel_convolution_avx2(channel_desc_t *desc, c16_t **input, c16_t **output, uint32_t length)
+{
+  float pathloss = powf(10.0f, (float)desc->path_loss_dB / 20.0f);
+  uint L = desc->channel_length;
+  uint64_t offset = desc->channel_offset;
+
+  cf_t rev_channel[desc->nb_rx][desc->nb_tx][L];
+
+  for (int rx = 0; rx < desc->nb_rx; rx++) {
+    for (int tx = 0; tx < desc->nb_tx; tx++) {
+      struct complexd *chan = desc->ch[rx + (tx * desc->nb_rx)];
+      for (uint i = 0U; i < L; i++) {
+        rev_channel[rx][tx][i].r = (float)chan[L - 1 - i].r * pathloss;
+        rev_channel[rx][tx][i].i = (float)chan[L - 1 - i].i * pathloss;
+      }
+    }
+  }
+
+  int end_idx = length - offset;
+
+  for (int rx = 0; rx < desc->nb_rx; rx++) {
+    int i = 0;
+
+    // 1. Transient part (scalar)
+    for (; i < (int)L - 1 && i < end_idx; i++) {
+      cf_t rx_sum = {0, 0};
+      for (int tx = 0; tx < desc->nb_tx; tx++) {
+        int oldest_idx = i - (L - 1);
+        int current_L = L;
+        int start_c = 0;
+        if (oldest_idx < 0) {
+          start_c = -oldest_idx;
+          current_L -= start_c;
+          oldest_idx = 0;
+        }
+        if (current_L > 0) {
+          c16_t *tx_ptr = &input[tx][oldest_idx];
+          cf_t *c_ptr = &rev_channel[rx][tx][start_c];
+          for (int l = 0; l < current_L; l++) {
+            rx_sum.r += (tx_ptr->r * c_ptr->r) - (tx_ptr->i * c_ptr->i);
+            rx_sum.i += (tx_ptr->i * c_ptr->r) + (tx_ptr->r * c_ptr->i);
+            tx_ptr++;
+            c_ptr++;
+          }
+        }
+      }
+      output[rx][i + offset].r = (int16_t)rx_sum.r;
+      output[rx][i + offset].i = (int16_t)rx_sum.i;
+    }
+
+    // 2. Steady state (SIMD)
+    for (; i <= end_idx - 4; i += 4) {
+      simde__m256 rx_sum_vec = simde_mm256_setzero_ps();
+
+      for (int tx = 0; tx < desc->nb_tx; tx++) {
+        int start_idx = i - (L - 1);
+        c16_t *tx_ptr = &input[tx][start_idx];
+        cf_t *c_ptr = &rev_channel[rx][tx][0];
+
+        for (int l = 0; l < (int)L; l++) {
+          simde__m128i in_128 = simde_mm_loadu_si128((simde__m128i const *)(tx_ptr + l)); // 4 complex samples
+          simde__m256i in_32 = simde_mm256_cvtepi16_epi32(in_128);
+          simde__m256 in_ps = simde_mm256_cvtepi32_ps(in_32);
+
+          float tr = c_ptr[l].r;
+          float ti = c_ptr[l].i;
+          simde__m256 tap = simde_mm256_set_ps(ti, tr, ti, tr, ti, tr, ti, tr);
+
+          simde__m256 in_re = simde_mm256_shuffle_ps(in_ps, in_ps, SIMDE_MM_SHUFFLE(2, 2, 0, 0));
+          simde__m256 in_im = simde_mm256_shuffle_ps(in_ps, in_ps, SIMDE_MM_SHUFFLE(3, 3, 1, 1));
+
+          simde__m256 term1 = simde_mm256_mul_ps(in_re, tap);
+          simde__m256 tap_swap = simde_mm256_set_ps(tr, -ti, tr, -ti, tr, -ti, tr, -ti);
+          simde__m256 term2 = simde_mm256_mul_ps(in_im, tap_swap);
+
+          rx_sum_vec = simde_mm256_add_ps(rx_sum_vec, simde_mm256_add_ps(term1, term2));
+        }
+      }
+      simde__m256i out_32 = simde_mm256_cvtps_epi32(rx_sum_vec);
+      simde__m128i out_16_low = simde_mm256_castsi256_si128(out_32);
+      simde__m128i out_16_high = simde_mm256_extracti128_si256(out_32, 1);
+      simde__m128i out_16 = simde_mm_packs_epi32(out_16_low, out_16_high);
+      simde_mm_storeu_si128((simde__m128i *)&output[rx][i + offset], out_16);
+    }
+
+    // 3. Remainder (scalar)
+    for (; i < end_idx; i++) {
+      cf_t rx_sum = {0, 0};
+      for (int tx = 0; tx < desc->nb_tx; tx++) {
+        int oldest_idx = i - (L - 1);
+        c16_t *tx_ptr = &input[tx][oldest_idx];
+        cf_t *c_ptr = &rev_channel[rx][tx][0];
+        for (int l = 0; l < (int)L; l++) {
+          rx_sum.r += (tx_ptr->r * c_ptr->r) - (tx_ptr->i * c_ptr->i);
+          rx_sum.i += (tx_ptr->i * c_ptr->r) + (tx_ptr->r * c_ptr->i);
+          tx_ptr++;
+          c_ptr++;
+        }
+      }
+      output[rx][i + offset].r = (int16_t)rx_sum.r;
+      output[rx][i + offset].i = (int16_t)rx_sum.i;
+    }
+  }
+}
+
+void channel_convolution_avx512(channel_desc_t *desc, c16_t **input, c16_t **output, uint32_t length)
+{
+  float pathloss = powf(10.0f, (float)desc->path_loss_dB / 20.0f);
+  uint L = desc->channel_length;
+  uint64_t offset = desc->channel_offset;
+
+  cf_t rev_channel[desc->nb_rx][desc->nb_tx][L];
+
+  for (int rx = 0; rx < desc->nb_rx; rx++) {
+    for (int tx = 0; tx < desc->nb_tx; tx++) {
+      struct complexd *chan = desc->ch[rx + (tx * desc->nb_rx)];
+      for (uint i = 0U; i < L; i++) {
+        rev_channel[rx][tx][i].r = (float)chan[L - 1 - i].r * pathloss;
+        rev_channel[rx][tx][i].i = (float)chan[L - 1 - i].i * pathloss;
+      }
+    }
+  }
+
+  int end_idx = length - offset;
+
+  for (int rx = 0; rx < desc->nb_rx; rx++) {
+    int i = 0;
+
+    // 1. Transient part (scalar)
+    for (; i < (int)L - 1 && i < end_idx; i++) {
+      cf_t rx_sum = {0, 0};
+      for (int tx = 0; tx < desc->nb_tx; tx++) {
+        int oldest_idx = i - (L - 1);
+        int current_L = L;
+        int start_c = 0;
+        if (oldest_idx < 0) {
+          start_c = -oldest_idx;
+          current_L -= start_c;
+          oldest_idx = 0;
+        }
+        if (current_L > 0) {
+          c16_t *tx_ptr = &input[tx][oldest_idx];
+          cf_t *c_ptr = &rev_channel[rx][tx][start_c];
+          for (int l = 0; l < current_L; l++) {
+            rx_sum.r += (tx_ptr->r * c_ptr->r) - (tx_ptr->i * c_ptr->i);
+            rx_sum.i += (tx_ptr->i * c_ptr->r) + (tx_ptr->r * c_ptr->i);
+            tx_ptr++;
+            c_ptr++;
+          }
+        }
+      }
+      output[rx][i + offset].r = (int16_t)rx_sum.r;
+      output[rx][i + offset].i = (int16_t)rx_sum.i;
+    }
+
+    // 2. Steady state (AVX512)
+    for (; i <= end_idx - 8; i += 8) {
+      simde__m512 rx_sum_vec = simde_mm512_setzero_ps();
+
+      for (int tx = 0; tx < desc->nb_tx; tx++) {
+        int start_idx = i - (L - 1);
+        c16_t *tx_ptr = &input[tx][start_idx];
+        cf_t *c_ptr = &rev_channel[rx][tx][0];
+
+        for (int l = 0; l < (int)L; l++) {
+          simde__m256i in_256i = simde_mm256_loadu_si256((simde__m256i const *)(tx_ptr + l));
+          simde__m512i in_512i = simde_mm512_cvtepi16_epi32(in_256i);
+          simde__m512 in_ps = simde_mm512_cvtepi32_ps(in_512i);
+
+          float tr = c_ptr[l].r;
+          float ti = c_ptr[l].i;
+          simde__m512 v_tr = simde_mm512_set1_ps(tr);
+          simde__m512 v_ti = simde_mm512_set1_ps(ti);
+          simde__m512 v_nti = simde_mm512_set1_ps(-ti);
+
+          simde__m512 tap = simde_mm512_mask_blend_ps(0xAAAA, v_tr, v_ti);
+          simde__m512 tap_swap = simde_mm512_mask_blend_ps(0xAAAA, v_nti, v_tr);
+
+          simde__m512 in_re = simde_mm512_shuffle_ps(in_ps, in_ps, 0xA0);
+          simde__m512 in_im = simde_mm512_shuffle_ps(in_ps, in_ps, 0xF5);
+
+          rx_sum_vec = simde_mm512_fmadd_ps(in_re, tap, rx_sum_vec);
+          rx_sum_vec = simde_mm512_fmadd_ps(in_im, tap_swap, rx_sum_vec);
+        }
+      }
+      simde__m512i out_32 = simde_mm512_cvtps_epi32(rx_sum_vec);
+      simde__m256i out_16 = simde_mm512_cvtsepi32_epi16(out_32);
+      simde_mm256_storeu_si256((simde__m256i *)&output[rx][i + offset], out_16);
+    }
+
+    // 3. Remainder (scalar)
+    for (; i < end_idx; i++) {
+      cf_t rx_sum = {0, 0};
+      for (int tx = 0; tx < desc->nb_tx; tx++) {
+        int oldest_idx = i - (L - 1);
+        c16_t *tx_ptr = &input[tx][oldest_idx];
+        cf_t *c_ptr = &rev_channel[rx][tx][0];
+        for (int l = 0; l < (int)L; l++) {
+          rx_sum.r += (tx_ptr->r * c_ptr->r) - (tx_ptr->i * c_ptr->i);
+          rx_sum.i += (tx_ptr->i * c_ptr->r) + (tx_ptr->r * c_ptr->i);
+          tx_ptr++;
+          c_ptr++;
+        }
+      }
+      output[rx][i + offset].r = (int16_t)rx_sum.r;
+      output[rx][i + offset].i = (int16_t)rx_sum.i;
+    }
+  }
+}

openair1/SIMULATION/TOOLS/sim.h

@@ -625,5 +625,9 @@ void do_DL_sig(sim_t *sim,
 
 void do_UL_sig(sim_t *sim, uint16_t subframe, uint8_t abstraction_flag, LTE_DL_FRAME_PARMS *frame_parms, uint32_t frame, int ru_id, uint8_t CC_id, int NB_UEs);
 int get_noise_power_dBFS(void);
+void channel_convolution(channel_desc_t *desc, c16_t **input, c16_t **output, uint32_t length);
+void channel_convolution_avx2(channel_desc_t *desc, c16_t **input, c16_t **output, uint32_t length);
+void channel_convolution_avx512(channel_desc_t *desc, c16_t **input, c16_t **output, uint32_t length);
+
 
 #endif