diff --git a/lib/include/srsran/phy/resampling/resampler.h b/lib/include/srsran/phy/resampling/resampler.h index 7a79cbb80..4d699d77d 100644 --- a/lib/include/srsran/phy/resampling/resampler.h +++ b/lib/include/srsran/phy/resampling/resampler.h @@ -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; /** diff --git a/lib/src/phy/resampling/resampler.c b/lib/src/phy/resampling/resampler.c index 6c398ae37..c4eba482b 100644 --- a/lib/src/phy/resampling/resampler.c +++ b/lib/src/phy/resampling/resampler.c @@ -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->mode = mode; + q->ratio = ratio; 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; -} \ No newline at end of file + return q->delay; +} diff --git a/lib/src/phy/resampling/test/resampler_test.c b/lib/src/phy/resampling/test/resampler_test.c index 9478f3ffc..df3061ed0 100644 --- a/lib/src/phy/resampling/test/resampler_test.c +++ b/lib/src/phy/resampling/test/resampler_test.c @@ -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);