PentHertz
/
srsLTE
mirror of https://github.com/PentHertz/srsLTE.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Fix resampler and improved unit test

This commit is contained in:
Xavier Arteaga 2021-05-11 11:44:17 +02:00 committed by Xavier Arteaga
parent 32228389a9
commit 9517b78c03
3 changed files with 129 additions and 41 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

23
lib/include/srsran/phy/resampling/resampler.h
View File

@ -40,19 +40,20 @@ typedef enum {
} srsran_resampler_mode_t;
/**
* Resampler internal buffers and subcomponents
* @brief Resampler internal buffers and subcomponents
*/
typedef struct {
srsran_resampler_mode_t mode;
uint32_t ratio;
uint32_t window_sz;
srsran_dft_plan_t fft;
srsran_dft_plan_t ifft;
uint32_t state_len;
cf_t* in_buffer;
cf_t* out_buffer;
cf_t* state;
cf_t* filter;
srsran_resampler_mode_t mode; ///< Interpolate or decimate mode
uint32_t ratio; ///< Decimation/Interpolation ratio
uint32_t window_sz; ///< Maximum number of processed samples
uint32_t delay; ///< Filter delay in samples
srsran_dft_plan_t fft; ///< Forward DFT
srsran_dft_plan_t ifft; ///< Backward DFT
uint32_t state_len; ///< Number of acccumulated samples in the internal state
cf_t* in_buffer; ///< DFT input buffer
cf_t* out_buffer; ///< DFT output buffer
cf_t* state; ///< Filter state
cf_t* filter; ///< Frequency domain filter
} srsran_resampler_fft_t;
/**

53
lib/src/phy/resampling/resampler.c
View File

@ -25,6 +25,11 @@
*/
#define RESAMPLER_BETA 0.45
/**
* Filter delay in multiples of ratio
*/
#define RESAMPLER_DELAY 7
/**
* The FFT size power is determined from the ratio logarithm in base 2 plus the following parameter
*/
@ -61,21 +66,26 @@ int srsran_resampler_fft_init(srsran_resampler_fft_t* q, srsran_resampler_mode_t
uint32_t output_fft_size = 0;
uint32_t high_size = base_size * ratio;
// Select FFT/IFFT sizes filter delay and window size. For best performance and avoid aliasing, the window size shall
// be as big as the input DFT subtracting the filter length at the input rate
switch (mode) {
case SRSRAN_RESAMPLER_MODE_INTERPOLATE:
input_fft_size = base_size;
output_fft_size = high_size;
q->delay = RESAMPLER_DELAY * ratio;
q->window_sz = input_fft_size - 2 * RESAMPLER_DELAY;
break;
case SRSRAN_RESAMPLER_MODE_DECIMATE:
default:
input_fft_size = high_size;
output_fft_size = base_size;
q->delay = RESAMPLER_DELAY * ratio;
q->window_sz = input_fft_size - 2 * q->delay;
break;
}
q->mode = mode;
q->ratio = ratio;
q->window_sz = input_fft_size / 4;
q->in_buffer = srsran_vec_cf_malloc(high_size);
if (q->in_buffer == NULL) {
@ -111,11 +121,20 @@ int srsran_resampler_fft_init(srsran_resampler_fft_t* q, srsran_resampler_mode_t
return SRSRAN_ERROR;
}
// Calculate absolute filter delay
double delay = (double)q->delay;
if (mode == SRSRAN_RESAMPLER_MODE_INTERPOLATE) {
delay = (double)(high_size - q->delay);
}
// Compute time domain filter coefficients, see raised cosine formula in section "1.2 Impulse Response" of
// https://dspguru.com/dsp/reference/raised-cosine-and-root-raised-cosine-formulas/
double T = (double)1.0;
for (int32_t i = 0; i < high_size; i++) {
double t = ((double)i - (double)high_size / 2.0) / (double)ratio;
// Convert to time
double t = ((double)i - delay) / (double)ratio;
// Compute coefficient
double h = 1.0 / T;
if (isnormal(t)) {
h = sin(M_PI * t / T);
@ -126,6 +145,11 @@ int srsran_resampler_fft_init(srsran_resampler_fft_t* q, srsran_resampler_mode_t
q->in_buffer[i] = (float)h;
}
if (srsran_verbose >= SRSRAN_VERBOSE_INFO && !handler_registered) {
printf("h_%s=", q->mode == SRSRAN_RESAMPLER_MODE_INTERPOLATE ? "interp" : "decimate");
srsran_vec_fprint_c(stdout, q->in_buffer, high_size);
}
// Compute frequency domain coefficients, since the filter is symmetrical, it does not matter whether FFT or iFFT
if (mode == SRSRAN_RESAMPLER_MODE_INTERPOLATE) {
srsran_dft_run_guru_c(&q->ifft);
@ -170,14 +194,11 @@ static void resampler_fft_interpolate(srsran_resampler_fft_t* q, const cf_t* inp
// Execute FFT
srsran_dft_run_guru_c(&q->fft);
// Replicate input spectrum
for (uint32_t i = 1; i < q->ratio; i++) {
srsran_vec_cf_copy(&q->out_buffer[q->fft.size * i], q->out_buffer, q->fft.size);
// Replicate input spectrum and filter at same time
for (uint32_t i = 0; i < q->ratio; i++) {
srsran_vec_prod_ccc(q->out_buffer, &q->filter[q->fft.size * i], &q->in_buffer[q->fft.size * i], q->fft.size);
}
// Apply filtering
srsran_vec_prod_ccc(q->out_buffer, q->filter, q->in_buffer, q->ifft.size);
// Execute iFFT
srsran_dft_run_guru_c(&q->ifft);
} else {
@ -223,12 +244,14 @@ static void resampler_fft_decimate(srsran_resampler_fft_t* q, const cf_t* input,
// Execute FFT
srsran_dft_run_guru_c(&q->fft);
// Apply filtering and cut
srsran_vec_prod_ccc(q->out_buffer, q->filter, q->in_buffer, q->ifft.size / 2);
srsran_vec_prod_ccc(&q->out_buffer[q->fft.size - q->ifft.size / 2],
&q->filter[q->fft.size - q->ifft.size / 2],
&q->in_buffer[q->ifft.size / 2],
q->ifft.size / 2);
// Apply filter
srsran_vec_prod_ccc(q->out_buffer, q->filter, q->out_buffer, q->fft.size);
// Decimate
srsran_vec_cf_copy(q->in_buffer, q->out_buffer, q->ifft.size);
for (uint32_t i = 1; i < q->ratio; i++) {
srsran_vec_sum_ccc(&q->out_buffer[q->ifft.size * i], q->in_buffer, q->in_buffer, q->ifft.size);
}
// Execute iFFT
srsran_dft_run_guru_c(&q->ifft);
@ -307,5 +330,5 @@ uint32_t srsran_resampler_fft_get_delay(srsran_resampler_fft_t* q)
return UINT32_MAX;
}
return q->ifft.size / 2;
return q->delay;
}

88
lib/src/phy/resampling/test/resampler_test.c
View File

@ -10,6 +10,7 @@
*
*/
#include "srsran/phy/channel/ch_awgn.h"
#include "srsran/phy/resampling/resampler.h"
#include "srsran/phy/utils/debug.h"
#include "srsran/phy/utils/vector.h"
@ -20,12 +21,19 @@
static uint32_t buffer_size = 1920;
static uint32_t factor = 2;
static uint32_t repetitions = 2;
static enum {
WAVE_SINE = 0,
WAVE_DELTA,
WAVE_STEP,
WAVE_GAUSS,
} wave = WAVE_SINE;
static void usage(char* prog)
{
printf("Usage: %s [sfr]\n", prog);
printf("\t-s Buffer size [Default %d]\n", buffer_size);
printf("\t-f Buffer size [Default %d]\n", factor);
printf("\t-f Interpolation/Decimation factor [Default %d]\n", factor);
printf("\t-w Wave type: sine, step, delta [Default sine]\n");
printf("\t-f r [Default %d]\n", repetitions);
}
@ -33,7 +41,7 @@ static void parse_args(int argc, char** argv)
{
int opt;
while ((opt = getopt(argc, argv, "sfr")) != -1) {
while ((opt = getopt(argc, argv, "sfrvw")) != -1) {
switch (opt) {
case 's':
buffer_size = (uint32_t)strtol(argv[optind], NULL, 10);
@ -44,6 +52,33 @@ static void parse_args(int argc, char** argv)
case 'r':
repetitions = (uint32_t)strtol(argv[optind], NULL, 10);
break;
case 'v':
srsran_verbose++;
break;
case 'w':
if (strcmp(argv[optind], "sine") == 0) {
wave = WAVE_SINE;
break;
}
if (strcmp(argv[optind], "delta") == 0) {
wave = WAVE_DELTA;
break;
}
if (strcmp(argv[optind], "step") == 0) {
wave = WAVE_STEP;
break;
}
if (strcmp(argv[optind], "gauss") == 0) {
wave = WAVE_GAUSS;
break;
}
printf("Invalid wave '%s'\n", argv[optind]);
usage(argv[0]);
break;
default:
usage(argv[0]);
exit(-1);
@ -56,6 +91,7 @@ int main(int argc, char** argv)
struct timeval t[3] = {};
srsran_resampler_fft_t interp = {};
srsran_resampler_fft_t decim = {};
srsran_channel_awgn_t awgn = {};
parse_args(argc, argv);
@ -72,7 +108,31 @@ int main(int argc, char** argv)
}
srsran_vec_cf_zero(src, buffer_size);
srsran_vec_gen_sine(1.0f, 0.01f, src, buffer_size / 10);
switch (wave) {
case WAVE_SINE:
srsran_vec_gen_sine(1.0f, 0.01f, src, buffer_size / 2);
break;
case WAVE_DELTA:
src[0] = 1.0f;
break;
case WAVE_STEP:
for (uint32_t i = 0; i < buffer_size; i++) {
src[i] = 1.0f;
}
break;
case WAVE_GAUSS:
srsran_channel_awgn_init(&awgn, 0);
srsran_channel_awgn_set_n0(&awgn, 0);
srsran_channel_awgn_run_c(&awgn, src, src, buffer_size);
srsran_channel_awgn_free(&awgn);
break;
}
if (srsran_verbose >= SRSRAN_VERBOSE_INFO && !handler_registered) {
printf("signal=");
srsran_vec_fprint_c(stdout, src, buffer_size);
}
gettimeofday(&t[1], NULL);
for (uint32_t r = 0; r < repetitions; r++) {
@ -82,22 +142,26 @@ int main(int argc, char** argv)
gettimeofday(&t[2], NULL);
get_time_interval(t);
uint64_t duration_us = (uint64_t)(t[0].tv_sec * 1000000UL + t[0].tv_usec);
printf("Done %.1f Msps\n", factor * buffer_size * repetitions / (double)duration_us);
// printf("interp=");
// srsran_vec_fprint_c(stdout, interpolated, buffer_size * factor);
if (srsran_verbose >= SRSRAN_VERBOSE_INFO && !handler_registered) {
printf("interp=");
srsran_vec_fprint_c(stdout, interpolated, buffer_size * factor);
printf("decim=");
srsran_vec_fprint_c(stdout, decimated, buffer_size);
}
// Check error
uint32_t delay = srsran_resampler_fft_get_delay(&decim) * 2;
uint32_t delay = (srsran_resampler_fft_get_delay(&decim) + srsran_resampler_fft_get_delay(&interp)) / factor;
uint32_t nsamples = buffer_size - delay;
srsran_vec_sub_ccc(src, &decimated[delay], interpolated, nsamples);
float mse = sqrtf(srsran_vec_avg_power_cf(interpolated, nsamples));
printf("MSE: %f\n", mse);
// printf("src=");
// srsran_vec_fprint_c(stdout, src, nsamples);
// printf("decim=");
// srsran_vec_fprint_c(stdout, &decimated[delay], nsamples);
if (srsran_verbose >= SRSRAN_VERBOSE_INFO && !handler_registered) {
printf("recovered=");
srsran_vec_fprint_c(stdout, &decimated[delay], nsamples);
}
printf("Done %.1f Msps; MSE: %.6f\n", factor * buffer_size * repetitions / (double)duration_us, mse);
srsran_resampler_fft_free(&interp);
srsran_resampler_fft_free(&decim);