Implemented high performance AWGN generator

This commit is contained in:
Xavier Arteaga 2020-03-22 23:42:34 +01:00 committed by Xavier Arteaga
parent b5be0b94b8
commit c107b04f5a
4 changed files with 431 additions and 2 deletions

View File

@ -27,18 +27,77 @@
* Reference:
*********************************************************************************************/
#include <complex.h>
#include <stdint.h>
#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

View File

@ -25,7 +25,172 @@
#include "gauss.h"
#include <srslte/phy/channel/ch_awgn.h>
#include <srslte/phy/utils/simd.h>
#include <srslte/phy/utils/vector.h>
#include <srslte/srslte.h>
#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)
{

View File

@ -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)

View File

@ -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 <complex.h>
#include <srslte/phy/channel/ch_awgn.h>
#include <srslte/phy/dft/dft.h>
#include <srslte/phy/utils/debug.h>
#include <unistd.h>
#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;
}