From c107b04f5a99ee56af1cb3aba2b0f079def25e2c Mon Sep 17 00:00:00 2001 From: Xavier Arteaga Date: Sun, 22 Mar 2020 23:42:34 +0100 Subject: [PATCH] Implemented high performance AWGN generator --- lib/include/srslte/phy/channel/ch_awgn.h | 63 +++++- lib/src/phy/channel/ch_awgn.c | 165 ++++++++++++++++ lib/src/phy/channel/test/CMakeLists.txt | 7 + lib/src/phy/channel/test/awgn_channel_test.c | 198 +++++++++++++++++++ 4 files changed, 431 insertions(+), 2 deletions(-) create mode 100644 lib/src/phy/channel/test/awgn_channel_test.c diff --git a/lib/include/srslte/phy/channel/ch_awgn.h b/lib/include/srslte/phy/channel/ch_awgn.h index 4b1630247..32a56d8cf 100644 --- a/lib/include/srslte/phy/channel/ch_awgn.h +++ b/lib/include/srslte/phy/channel/ch_awgn.h @@ -27,18 +27,77 @@ * Reference: *********************************************************************************************/ -#include #include - #include "srslte/config.h" #ifndef SRSLTE_CH_AWGN_H #define SRSLTE_CH_AWGN_H +#ifdef __cplusplus +extern "C" { +#endif + +/** + * The srsLTE channel AWGN implements an efficient Box-Muller Method accelerated with SIMD. + */ +typedef struct { + cf_t* table_exp; + float* table_log; + uint32_t rand_state; + float std_dev; +} srslte_channel_awgn_t; + +/** + * Initialization function of the channel AWGN object + * + * @param q AWGN channel object + * @param seed random generator seed + */ +SRSLTE_API int srslte_channel_awgn_init(srslte_channel_awgn_t* q, uint32_t seed); + +/** + * Sets the noise level N0 in decibels full scale (dBfs) + * + * @param q AWGN channel object + * @param n0_dBfs noise level + */ +SRSLTE_API int srslte_channel_awgn_set_n0(srslte_channel_awgn_t* q, float n0_dBfs); + +/** + * Runs the complex AWGN channel + * + * @param q AWGN channel object + * @param in complex input array + * @param out complex output array + * @param length number of samples + */ +SRSLTE_API void srslte_channel_awgn_run_c(srslte_channel_awgn_t* q, const cf_t* in, cf_t* out, uint32_t length); + +/** + * Runs the real AWGN channel + * + * @param q AWGN channel object + * @param in real input array + * @param out real output array + * @param length number of samples + */ +SRSLTE_API void srslte_channel_awgn_run_f(srslte_channel_awgn_t* q, const float* in, float* out, uint32_t length); + +/** + * Free AWGN channel generator data + * + * @param q AWGN channel object + */ +SRSLTE_API void srslte_channel_awgn_free(srslte_channel_awgn_t* q); + SRSLTE_API void srslte_ch_awgn_c(const cf_t* input, cf_t* output, float variance, uint32_t len); SRSLTE_API void srslte_ch_awgn_f(const float* x, float* y, float variance, uint32_t len); SRSLTE_API float srslte_ch_awgn_get_variance(float ebno_db, float rate); +#ifdef __cplusplus +} +#endif + #endif // SRSLTE_CH_AWGN_H diff --git a/lib/src/phy/channel/ch_awgn.c b/lib/src/phy/channel/ch_awgn.c index 4b7d967e5..f43752c7e 100644 --- a/lib/src/phy/channel/ch_awgn.c +++ b/lib/src/phy/channel/ch_awgn.c @@ -25,7 +25,172 @@ #include "gauss.h" #include +#include #include +#include + +#define AWGN_TABLE_SIZE_POW 10U +#define AWGN_TABLE_SIZE (1U << AWGN_TABLE_SIZE_POW) +#define AWGN_TABLE_ALLOC_SIZE (AWGN_TABLE_SIZE + SRSLTE_SIMD_F_SIZE) + +// Linear congruential Generator parameters +#define AWGN_LCG_A 12345U +#define AWGN_LCG_C 1103515245U + +#undef ENABLE_GUI /* Disable GUI by default */ +#ifdef ENABLE_GUI +#include "srsgui/srsgui.h" +static bool enable_gui = false; +#endif /* ENABLE_GUI */ + +static inline int32_t channel_awgn_rand(srslte_channel_awgn_t* q) +{ + q->rand_state = (AWGN_LCG_A + AWGN_LCG_C * q->rand_state) & (AWGN_TABLE_SIZE - 1); + return q->rand_state; +} + +int srslte_channel_awgn_init(srslte_channel_awgn_t* q, uint32_t seed) +{ + if (!q) { + return SRSLTE_ERROR_INVALID_INPUTS; + } + + // Initialise random generator + q->rand_state = seed; + + // Allocate complex exponential and logarithmic tables + q->table_exp = srslte_vec_cf_malloc(AWGN_TABLE_ALLOC_SIZE); + q->table_log = srslte_vec_f_malloc(AWGN_TABLE_ALLOC_SIZE); + if (!q->table_exp || !q->table_log) { + ERROR("Malloc\n"); + } + + // Fill tables + for (uint32_t i = 0; i < AWGN_TABLE_SIZE; i++) { + float temp1 = (float)i / (float)AWGN_TABLE_SIZE; + float temp2 = (float)(i + 1) / (float)AWGN_TABLE_SIZE; + + q->table_exp[i] = cexpf(I * 2.0f * (float)M_PI * temp1); + q->table_log[i] = sqrtf(-2.0f * logf(temp2)); + } + + // Shuffle values in tables to break correlation in SIMD registers + for (uint32_t i = 0; i < AWGN_TABLE_SIZE; i++) { + int32_t idx; + + do { + idx = channel_awgn_rand(q) % AWGN_TABLE_SIZE; + } while (idx == i); + float temp_log = q->table_log[i]; + q->table_log[i] = q->table_log[idx]; + q->table_log[idx] = temp_log; + + do { + idx = channel_awgn_rand(q) % AWGN_TABLE_SIZE; + } while (idx == i); + cf_t temp_exp = q->table_exp[i]; + q->table_exp[i] = q->table_exp[idx]; + q->table_exp[idx] = temp_exp; + } + + // Copy head in tail for keeping continuity in SIMD registers + for (uint32_t i = 0; i < SRSLTE_SIMD_F_SIZE; i++) { + q->table_log[i + AWGN_TABLE_SIZE] = q->table_log[i]; + q->table_exp[i + AWGN_TABLE_SIZE] = q->table_exp[i]; + } + + return SRSLTE_SUCCESS; +} + +int srslte_channel_awgn_set_n0(srslte_channel_awgn_t* q, float n0_dBfs) +{ + if (!q) { + return SRSLTE_ERROR_INVALID_INPUTS; + } + + q->std_dev = powf(10.0f, n0_dBfs / 20.0f); + + return isnormal(q->std_dev) ? SRSLTE_SUCCESS : SRSLTE_ERROR; +} + +static inline void channel_awgn_run(srslte_channel_awgn_t* q, const float* in, float* out, uint32_t size, float std_dev) +{ + if (!q) { + return; + } + + int i = 0; + +#if SRSLTE_SIMD_F_SIZE + for (; i < (int)size - SRSLTE_SIMD_F_SIZE + 1; i += SRSLTE_SIMD_F_SIZE) { + // Load SIMD registers + simd_f_t t1 = srslte_simd_f_loadu((float*)&q->table_exp[channel_awgn_rand(q)]); + simd_f_t t2 = srslte_simd_f_loadu(&q->table_log[channel_awgn_rand(q)]); + simd_f_t in_ = srslte_simd_f_loadu(&in[i]); + simd_f_t out_ = srslte_simd_f_set1(std_dev); + + // Compute random number + out_ = srslte_simd_f_mul(t1, out_); + out_ = srslte_simd_f_mul(t2, out_); + out_ = srslte_simd_f_add(in_, out_); + + // Store random number + srslte_simd_f_storeu(&out[i], out_); + } +#endif /* SRSLTE_SIMD_F_SIZE */ + + int32_t idx1 = channel_awgn_rand(q); + int32_t idx2 = channel_awgn_rand(q); + for (; i < size; i++) { + + if (i % 8 == 0) { + idx1 = channel_awgn_rand(q); + idx2 = channel_awgn_rand(q); + } else { + idx1 = (idx1 + 1) % AWGN_TABLE_SIZE; + idx2 = (idx2 + 1) % AWGN_TABLE_SIZE; + } + + float n = std_dev; + n *= q->table_log[idx1]; + + cf_t t_exp = q->table_exp[idx2]; + float n1 = n * __real__ t_exp; + float n2 = n * __imag__ t_exp; + + out[i] = in[i] + n1; + i++; + + if (i < size) { + out[i] = in[i] + n2; + } + } +} + +void srslte_channel_awgn_run_c(srslte_channel_awgn_t* q, const cf_t* in, cf_t* out, uint32_t size) +{ + channel_awgn_run(q, (float*)in, (float*)out, 2 * size, q->std_dev * (float)M_SQRT1_2); +} + +void srslte_channel_awgn_run_f(srslte_channel_awgn_t* q, const float* in, float* out, uint32_t size) +{ + channel_awgn_run(q, in, out, size, q->std_dev); +} + +void srslte_channel_awgn_free(srslte_channel_awgn_t* q) +{ + if (!q) { + return; + } + + if (q->table_exp) { + free(q->table_exp); + } + + if (q->table_log) { + free(q->table_log); + } +} float srslte_ch_awgn_get_variance(float ebno_db, float rate) { diff --git a/lib/src/phy/channel/test/CMakeLists.txt b/lib/src/phy/channel/test/CMakeLists.txt index a6db348dd..e904cde55 100644 --- a/lib/src/phy/channel/test/CMakeLists.txt +++ b/lib/src/phy/channel/test/CMakeLists.txt @@ -37,3 +37,10 @@ add_executable(hst_channel_test hst_channel_test.c) target_link_libraries(hst_channel_test srslte_phy srslte_common srslte_phy ${SEC_LIBRARIES} ${CMAKE_THREAD_LIBS_INIT}) add_test(hst_channel_test hst_channel_test -f 750 -t 7.2 -i 0 -T 1 -s 1.92e6) +add_executable(awgn_channel_test awgn_channel_test.c) +if(SRSGUI_FOUND) + target_link_libraries(awgn_channel_test ${SRSGUI_LIBRARIES}) +endif(SRSGUI_FOUND) +target_link_libraries(awgn_channel_test srslte_phy srslte_common srslte_phy ${SEC_LIBRARIES} ${CMAKE_THREAD_LIBS_INIT}) +add_test(awgn_channel_test awgn_channel_test) + diff --git a/lib/src/phy/channel/test/awgn_channel_test.c b/lib/src/phy/channel/test/awgn_channel_test.c new file mode 100644 index 000000000..3cfde7af6 --- /dev/null +++ b/lib/src/phy/channel/test/awgn_channel_test.c @@ -0,0 +1,198 @@ +/* + * Copyright 2013-2020 Software Radio Systems Limited + * + * This file is part of srsLTE. + * + * srsLTE is free software: you can redistribute it and/or modify + * it under the terms of the GNU Affero General Public License as + * published by the Free Software Foundation, either version 3 of + * the License, or (at your option) any later version. + * + * srsLTE is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Affero General Public License for more details. + * + * A copy of the GNU Affero General Public License can be found in + * the LICENSE file in the top-level directory of this distribution + * and at http://www.gnu.org/licenses/. + * + */ + +#include "srslte/phy/utils/vector.h" +#include +#include +#include +#include +#include + +#undef ENABLE_GUI +#ifdef ENABLE_GUI +#include "srsgui/srsgui.h" +#endif /* ENABLE_GUI */ + +static srslte_channel_awgn_t awgn = {}; + +static uint32_t nof_samples = 1920 * 8; +static float n0_min = 0.0f; +static float n0_max = 6.0f; +static float n0_step = 0.25f; +static float tolerance = 0.05f; + +static void usage(char* prog) +{ + printf("Usage: %s [nmMst]\n", prog); + printf("\t-n number of samples to simulate [Default %d]\n", nof_samples); + printf("\t-m Minimum n0 (in dBfs) to simulate [Default %.3f]\n", n0_min); + printf("\t-M Maximum n0 (in dBfs) to simulate: [Default %.3f]\n", n0_max); + printf("\t-s n0 step size: [Default %.3f]\n", n0_step); + printf("\t-t tolerance: [Default %.3f]\n", tolerance); +} + +static void parse_args(int argc, char** argv) +{ + int opt; + while ((opt = getopt(argc, argv, "nmMst")) != -1) { + switch (opt) { + case 'm': + n0_min = strtof(argv[optind], NULL); + break; + case 'M': + n0_max = strtof(argv[optind], NULL); + break; + case 's': + n0_step = strtof(argv[optind], NULL); + break; + case 't': + tolerance = strtof(argv[optind], NULL); + break; + case 'n': + nof_samples = (uint32_t)strtol(argv[optind], NULL, 10); + break; + default: + usage(argv[0]); + exit(-1); + } + } +} + +int main(int argc, char** argv) +{ + int ret = SRSLTE_SUCCESS; + cf_t* input_buffer = NULL; + cf_t* output_buffer = NULL; + uint64_t count_samples = 0; + uint64_t count_us = 0; + + // Parse arguments + parse_args(argc, argv); + + // Initialise buffers + input_buffer = srslte_vec_cf_malloc(nof_samples); + output_buffer = srslte_vec_cf_malloc(nof_samples); + + if (!input_buffer || !output_buffer) { + ERROR("Error: Allocating memory\n"); + ret = SRSLTE_ERROR; + } + + // Initialise input + srslte_vec_cf_zero(input_buffer, nof_samples); + +#ifdef ENABLE_GUI + sdrgui_init(); + sdrgui_init_title("SRS Fading channel"); + plot_real_t plot_fft = NULL; + + plot_real_init(&plot_fft); + plot_real_setTitle(&plot_fft, "AWGN"); + plot_real_addToWindowGrid(&plot_fft, (char*)"Spectrum", 0, 0); + + cf_t* fft_out = srslte_vec_cf_malloc(nof_samples); + srslte_dft_plan_t fft = {}; + if (srslte_dft_plan_c(&fft, nof_samples, SRSLTE_DFT_FORWARD)) { + ERROR("Error: init DFT\n"); + ret = SRSLTE_ERROR; + } +#endif /* ENABLE_GUI */ + + // Initialise AWGN channel + if (ret == SRSLTE_SUCCESS) { + ret = srslte_channel_awgn_init(&awgn, 0); + } + + float n0 = n0_min; + while (!isnan(n0) && !isinf(n0) && n0 < n0_max) { + struct timeval t[3] = {}; + + srslte_channel_awgn_set_n0(&awgn, n0); + + // Run actual test + gettimeofday(&t[1], NULL); + srslte_channel_awgn_run_c(&awgn, input_buffer, output_buffer, nof_samples); + gettimeofday(&t[2], NULL); + get_time_interval(t); + + float power = srslte_vec_avg_power_cf(output_buffer, nof_samples); + float power_dB = srslte_convert_power_to_dB(power); + + if ((n0 + tolerance) < power_dB || (n0 - tolerance) > power_dB) { + printf("-- failed: %.3f<%.3f<%.3f\n", n0 - tolerance, power_dB, n0 + tolerance); + ret = SRSLTE_ERROR; + } + +#ifdef ENABLE_GUI + srslte_dft_run_c(&fft, output_buffer, fft_out); + + float min = +INFINITY; + float max = -INFINITY; + float* fft_mag = (float*)fft_out; + for (uint32_t i = 0; i < nof_samples; i++) { + float mag = srslte_convert_amplitude_to_dB(cabsf(fft_out[i])); + + min = SRSLTE_MIN(min, mag); + max = SRSLTE_MAX(max, mag); + + fft_mag[i] = srslte_convert_amplitude_to_dB(mag); + } + + plot_real_setYAxisScale(&plot_fft, min, max); + + plot_real_setNewData(&plot_fft, fft_mag, nof_samples); + sleep(5); +#endif /* ENABLE_GUI */ + + count_samples += nof_samples; + count_us += t->tv_usec + t->tv_sec * 1000000; + + n0 += n0_step; + } + + // Free + srslte_channel_awgn_free(&awgn); + + if (input_buffer) { + free(input_buffer); + } + + if (output_buffer) { + free(output_buffer); + } + +#ifdef ENABLE_GUI + if (fft_out) { + free(fft_out); + } + srslte_dft_plan_free(&fft); +#endif /* ENABLE_GUI */ + + // Print result and exit + printf("Test n0_min=%.3f; n0_max=%.3f; n0_step=%.3f; nof_samples=%d; %s ... %.1f MSps\n", + n0_min, + n0_max, + n0_step, + nof_samples, + (ret == SRSLTE_SUCCESS) ? "Passed" : "Failed", + (double)nof_samples / (double)count_us); + return ret; +}