From f5a367464f3b63033a3c39a0773abbe7d28f293e Mon Sep 17 00:00:00 2001 From: Thomas Kolb Date: Wed, 17 May 2023 22:28:18 +0200 Subject: [PATCH] Optimized one-shot frequency synchronization - Ramp-up length increased to 128 symbols (here is room for improvement!) - Try to detect the frequency once during ramp-up. To do so, every second symbol is inverted (to remove the +/-1 symbol toggling) and the phase difference between neigboring resulting symbols is checked. When it is low enough for all symbols, the frequency is estimated and corrected. When frequency estimation was done, it is not retried for a number of incoming symbols in order to allow the timing estimator to converge again. - This approach was verified in a simulated loopback test with frequency offset and AWGN. --- impl/make.sh | 2 +- impl/src/config.h | 4 +- impl/src/layer1/correlator.c | 30 ++--------- impl/src/layer1/freq_est.c | 54 ++++++++++++++++++-- impl/src/layer1/freq_est.h | 17 +++++++ impl/src/layer1/rx.c | 86 ++++++++++++++++++++++---------- impl/src/main.c | 10 ++++ impl/src/results.h | 4 +- impl/test/layer1/test_loopback.c | 8 +-- 9 files changed, 150 insertions(+), 65 deletions(-) diff --git a/impl/make.sh b/impl/make.sh index e80f080..336e931 100755 --- a/impl/make.sh +++ b/impl/make.sh @@ -2,5 +2,5 @@ mkdir -p build cd build -cmake .. +cmake -DCMAKE_EXPORT_COMPILE_COMMANDS=true .. make diff --git a/impl/src/config.h b/impl/src/config.h index 6a33984..f82a9bb 100644 --- a/impl/src/config.h +++ b/impl/src/config.h @@ -29,7 +29,7 @@ #define RRC_DELAY 7 // delay in symbols #define RRC_BETA 0.2f -#define TRANSMISSION_RAMP_UP_LEN 64 // symbols +#define TRANSMISSION_RAMP_UP_LEN 128 // symbols #define TRANSMISSION_RAMP_DOWN_LEN 32 // symbols /*** SDR CONFIG ***/ @@ -43,7 +43,7 @@ // actually transmitted or received signal frequencies, NOT the SDR center frequency. #define SDR_TX_FREQ 434.100e6f -#define SDR_RX_FREQ 434.100e6f +#define SDR_RX_FREQ 434.115e6f // shift applied in the baseband, to get rid of SDR DC peak. If the value here // is not 0, software mixing will be done on the received signal. diff --git a/impl/src/layer1/correlator.c b/impl/src/layer1/correlator.c index 5db3df2..d69ac00 100644 --- a/impl/src/layer1/correlator.c +++ b/impl/src/layer1/correlator.c @@ -5,6 +5,8 @@ #include "freq_est.h" +#include "utils.h" + #include "correlator.h" bool correlator_init(correlator_ctx_t *ctx, const float complex *search_seq, size_t nsymbols) @@ -142,34 +144,8 @@ float correlator_get_mean_frequency_deviation(correlator_ctx_t *ctx, size_t L, f size_t m0 = ctx->search_sequence_len - L; for(size_t m = m0; m < ctx->search_sequence_len; m++) { size_t nm = (n + m) & ctx->buffer_size_mask; - z[m - m0] = ctx->search_sequence[m] * ctx->input_history[nm]; + z[m - m0] = conjf(ctx->search_sequence[m]) * ctx->input_history[nm]; } - /* Louise & Regiannini frequency estimator */ - /* - // Calculate averaged phase increments for sub-sequences - size_t N = L/2; - float complex R_kappa[N]; - - for(size_t kappa = 0; kappa < N; kappa++) { - for(size_t k = 0; k < L - kappa; k++) { - R_kappa[kappa] += z[k + kappa] * conj(z[k]); - } - - R_kappa[kappa] /= L - kappa; - } - - // Calculate phase estimate (in radians/sample) - float complex sum_R_kappa = 0; - - for(size_t kappa = 0; kappa < N; kappa++) { - sum_R_kappa += R_kappa[kappa]; - } - - float arg = cargf(sum_R_kappa); - - return arg / (1 + N); - */ - return freq_est_data_free(z, L, phase_offset); } diff --git a/impl/src/layer1/freq_est.c b/impl/src/layer1/freq_est.c index 7011bd4..51f2041 100644 --- a/impl/src/layer1/freq_est.c +++ b/impl/src/layer1/freq_est.c @@ -1,9 +1,11 @@ #include +#include #include "freq_est.h" #include "utils.h" +#if 0 // Fits a line y' = m*x' + t to the given points (x, y), where x is the // zero-based index in the array. static void linear_regression(size_t nx, const float *y, float *m, float *t) @@ -27,7 +29,6 @@ static void linear_regression(size_t nx, const float *y, float *m, float *t) *t = mean_y - (*m) * mean_x; } -#if 0 // Calculates an estimation of the true peak frequency in the signal by // calculating the FFT und interpolating around the bin with the highest value. static float fft_true_peak(const float complex *recv, size_t n) @@ -100,7 +101,7 @@ float freq_est_lr(const float complex *recv, size_t l0) float arg = cargf(R_total); // convert to frequency offset in radians per sample - float freq_offset = arg / (N+1) / 2.0f; /* FIXME: not sure where the factor 2 comes from. */ + float freq_offset = 2.0f * arg / (N+1); /* FIXME: check where the factor 2.0 comes from… */ return freq_offset; } @@ -128,6 +129,42 @@ float freq_est_unknown_bpsk(const float complex *recv, size_t n, float *final_ph } +float freq_est_in_rampup(const float complex *recv, size_t n, float *final_phase) +{ + float complex rotated[n]; + + float prev_phase = 0.0f; + + // rotate every second symbol by 180°. As a plausibility check that we are + // in fact inside the ramp-up, we verify that there is no phase rotation by + // more than 90° between the symbols. + for(size_t i = 0; i < n; i++) { + if((i % 2) == 0) { + rotated[i] = recv[i]; + } else { + rotated[i] = -recv[i]; + } + + float phase = cargf(rotated[i]); + if(i > 0) { + float phase_delta = phase - prev_phase; + if(fabsf(phase_delta) > 3.14159f/4.0f) { + // abort because we either have the wrong signal or too much noise. + if(final_phase) { + *final_phase = 0.0f; + } + return 0.0f; + } + } + prev_phase = phase; + } + + // signal is clean (and all symbol information removed), so we run the estimator + + return freq_est_data_free(rotated, n, final_phase); +} + + float freq_est_known_symbols(const float complex *recv, const float complex *ref, size_t n, float *final_phase) { float complex symbols[n]; @@ -140,6 +177,7 @@ float freq_est_known_symbols(const float complex *recv, const float complex *ref return freq_est_data_free(symbols, n, final_phase); } + float freq_est_data_free(const float complex *symbols, size_t n, float *final_phase) { float phases[n]; @@ -151,10 +189,18 @@ float freq_est_data_free(const float complex *symbols, size_t n, float *final_ph liquid_unwrap_phase(phases, n); float freq, phi0; - linear_regression(n, phases, &freq, &phi0); + freq = freq_est_lr(symbols, n); + + // calculate the mean of unwrapped symbol phases and from that the + // de-noised phase at the end of the sequence. + phi0 = 0.0f; + for(size_t i = 0; i < n; i++) { + phi0 += phases[i]; + } + phi0 /= n; if(final_phase) { - *final_phase = phi0 + n * freq; + *final_phase = phi0 + (n+1)/2.0f * freq; } return freq; diff --git a/impl/src/layer1/freq_est.h b/impl/src/layer1/freq_est.h index 6ee29d0..5c44a7e 100644 --- a/impl/src/layer1/freq_est.h +++ b/impl/src/layer1/freq_est.h @@ -35,6 +35,23 @@ float freq_est_unknown_bpsk(const float complex *recv, size_t n, float *final_ph */ float freq_est_known_symbols(const float complex *recv, const float complex *ref, size_t n, float *final_phase); +/*! + * \brief Estimate the frequency during ramp-up. + * + * This is a special version for frequency estimation during signal ramp-up + * where alternating +1/-1 BPSK symbols are transmitted. Therefore, every + * second symbol is inverted (rotated by 180°). As a plausibility check, the + * phase between neighboring symbols is calculated after rotation and if it + * changes by more than 90°, frequency estimation is skipped. + * + * \param recv Pointer to the array of received symbols. + * \param n Number of symbols given. + * \param final_phase Pointer to a float where the phase correction value is + * stored. May be NULL if that value is not required. + * \returns The estimated frequency in radians per sample. + */ +float freq_est_in_rampup(const float complex *recv, size_t n, float *final_phase); + /*! * \brief Estimate the frequency from data-free symbols. * diff --git a/impl/src/layer1/rx.c b/impl/src/layer1/rx.c index 5c7e613..fdcfbf2 100644 --- a/impl/src/layer1/rx.c +++ b/impl/src/layer1/rx.c @@ -15,28 +15,32 @@ #include "utils.h" #define HEADER_SIZE_BYTES 4 -#define FREQ_EST_L 8 +#define FREQ_EST_L 24 #define AGC_BW_ACQUISITION 5e-2f #define AGC_BW_TRACKING 1e-4f #define MAX_COARSE_FREQ_OFFSET 0.20f -#define SHOW_DEBUG_LOG 0 +#define SHOW_DEBUG_LOG 1 #if SHOW_DEBUG_LOG -# define DEBUG_LOG(...) fprintf(stderr, __VA_ARGS__) +# define DEBUG_LOG(...) fprintf(stderr, "DBG: " __VA_ARGS__) #else # define DEBUG_LOG(...) {} #endif static void update_nco_pll(nco_crcf nco, float phase_error) { - static const float pll_alpha = 0.20f; // phase adjustment factor + static const float pll_alpha = 0.03f; // phase adjustment factor static const float pll_beta = (pll_alpha * pll_alpha) / 2.0f; // frequency adjustment factor nco_crcf_adjust_phase(nco, pll_alpha * phase_error); nco_crcf_adjust_frequency(nco, pll_beta * phase_error); + + DEBUG_LOG("NCO adjusted: f=%.6f, φ=%.6f\n", + nco_crcf_get_frequency(nco), + nco_crcf_get_phase(nco)); } @@ -45,6 +49,8 @@ static bool acquire_preamble(layer1_rx_t *rx, const float complex sample) static float complex freq_est_history[FREQ_EST_L]; static unsigned int freq_est_history_write_idx = 0; + static uint32_t freq_est_holdoff_samples = 0; + // preamble search float complex corr_out = correlator_step(&rx->preamble_correlator, sample); @@ -59,10 +65,12 @@ static bool acquire_preamble(layer1_rx_t *rx, const float complex sample) DEBUG_LOG("Preamble frequency deviation: %.6f rad/symbol\n", mean_frequency_error); // Set the fine carrier correction NCO to the values estimated from the preamble - nco_crcf_set_frequency(rx->carrier_fine_nco, mean_frequency_error / RRC_SPS * 0.5f); + nco_crcf_set_frequency(rx->carrier_fine_nco, mean_frequency_error); nco_crcf_set_phase(rx->carrier_fine_nco, phase_offset); - DEBUG_LOG("New estimated carrier frequency: %.6f + %.6f\n", nco_crcf_get_frequency(rx->carrier_coarse_nco), nco_crcf_get_frequency(rx->carrier_fine_nco)); + float fcoarse = nco_crcf_get_frequency(rx->carrier_coarse_nco); + float ffine = nco_crcf_get_frequency(rx->carrier_fine_nco); + DEBUG_LOG("New estimated carrier frequency: %.6f + %.6f/4 = %.6f rad/sample\n", fcoarse, ffine, fcoarse+ffine/4); //float complex input_history[preamble_get_symbol_count()]; //correlator_get_input_history(&rx->preamble_correlator, input_history); @@ -79,32 +87,44 @@ static bool acquire_preamble(layer1_rx_t *rx, const float complex sample) return true; // preamble found! } else { // preamble not found. + //DEBUG_LOG("Preamble not found: %.3f < %.3f\n", cabsf(corr_out), 0.5f * preamble_get_symbol_count()); // Run the frequency estimator for every incoming sample overy the last // FREQ_EST_L samples. This is an implementation that works with unknown // BPSK symbols and therefore can be used during ramp-up and preamble. - float complex linearized_history[FREQ_EST_L]; - for(size_t i = 0; i < FREQ_EST_L; i++) { - linearized_history[i] = freq_est_history[(i + freq_est_history_write_idx) % FREQ_EST_L]; + if(freq_est_holdoff_samples == 0) { //freq_est_history_write_idx == FREQ_EST_L) { + float complex linearized_history[FREQ_EST_L]; + for(size_t i = 0; i < FREQ_EST_L; i++) { + linearized_history[i] = freq_est_history[(i + freq_est_history_write_idx) % FREQ_EST_L]; + } + + float freq_est = freq_est_in_rampup(freq_est_history, FREQ_EST_L, NULL); + + float freq_adjustment = freq_est / RRC_SPS; + nco_crcf_adjust_frequency(rx->carrier_coarse_nco, freq_adjustment); + + float actual_freq = nco_crcf_get_frequency(rx->carrier_coarse_nco); + if(actual_freq > MAX_COARSE_FREQ_OFFSET) { + //fprintf(stderr, ">"); + nco_crcf_set_frequency(rx->carrier_coarse_nco, MAX_COARSE_FREQ_OFFSET); + } else if(actual_freq < -MAX_COARSE_FREQ_OFFSET) { + //fprintf(stderr, "<"); + nco_crcf_set_frequency(rx->carrier_coarse_nco, -MAX_COARSE_FREQ_OFFSET); + } + + //freq_est_history_write_idx = 0; + + if(freq_est != 0.0f) { + DEBUG_LOG("Freq. est (x%d): %0.6f - Adj (x1): %.6f - carrier frequency (x1): %.6f\n", + RRC_SPS, freq_est, freq_adjustment, nco_crcf_get_frequency(rx->carrier_coarse_nco)); + + freq_est_holdoff_samples = preamble_get_symbol_count() * 2; + } + } else { + freq_est_holdoff_samples--; } - float freq_est = freq_est_unknown_bpsk(linearized_history, FREQ_EST_L, NULL); - - float freq_adjustment = freq_est / 32.0f; - nco_crcf_adjust_frequency(rx->carrier_coarse_nco, freq_adjustment); - - float actual_freq = nco_crcf_get_frequency(rx->carrier_coarse_nco); - if(actual_freq > MAX_COARSE_FREQ_OFFSET) { - //fprintf(stderr, ">"); - nco_crcf_set_frequency(rx->carrier_coarse_nco, MAX_COARSE_FREQ_OFFSET); - } else if(actual_freq < -MAX_COARSE_FREQ_OFFSET) { - //fprintf(stderr, "<"); - nco_crcf_set_frequency(rx->carrier_coarse_nco, -MAX_COARSE_FREQ_OFFSET); - } - - DEBUG_LOG("Freq. est: %0.6f - Adj: %.6f - carrier frequency: %.6f\n", freq_est, freq_adjustment, nco_crcf_get_frequency(rx->carrier_coarse_nco)); - assert(freq_est_history_write_idx < FREQ_EST_L); freq_est_history[freq_est_history_write_idx] = sample; @@ -123,11 +143,15 @@ result_t layer1_rx_process(layer1_rx_t *rx, const float complex *samples, size_t float complex samples2dump[sample_count]; size_t nsamples2dump = 0; + float complex samples2dump2[sample_count]; + size_t nsamples2dump2 = 0; // filter the input float complex filtered_samples[sample_count]; firfilt_crcf_execute_block(rx->channel_filter, (float complex *)samples, sample_count, filtered_samples); + fprintf(stderr, "\nagc: %f\n", agc_crcf_get_gain(rx->agc)); + for(unsigned int i = 0; i < sample_count; i++) { rx_state_t last_state = rx->state; @@ -149,6 +173,7 @@ result_t layer1_rx_process(layer1_rx_t *rx, const float complex *samples, size_t symsync_crcf_execute(rx->symsync, &mixed_sample, 1, &symsync_out, &out_len); if(out_len != 0) { + switch(rx->state) { // Try to acquire packets by synchronizing the frequency // (symbol-independent search) and correlating the preamble. @@ -175,13 +200,16 @@ result_t layer1_rx_process(layer1_rx_t *rx, const float complex *samples, size_t break; case RX_STATE_HEADER: + nco_crcf_step(rx->carrier_fine_nco); nco_crcf_mix_down(rx->carrier_fine_nco, symsync_out, &mixed_sample); + if(symbol_counter < rx->hdr_len_symbols) { unsigned int sym_demod; modem_demodulate(rx->hdr_demod, mixed_sample, &sym_demod); float phase_error = modem_get_demodulator_phase_error(rx->hdr_demod); - DEBUG_LOG("Sym: %d; Phase error: %f\n", sym_demod, phase_error); + //DEBUG_LOG("Sym: %d; Phase error: %f\n", sym_demod, phase_error); + update_nco_pll(rx->carrier_fine_nco, phase_error); @@ -230,13 +258,16 @@ result_t layer1_rx_process(layer1_rx_t *rx, const float complex *samples, size_t break; case RX_STATE_DATA: + nco_crcf_step(rx->carrier_fine_nco); nco_crcf_mix_down(rx->carrier_fine_nco, symsync_out, &mixed_sample); + + samples2dump2[nsamples2dump2++] = mixed_sample; if(symbol_counter < rx->payload_len_symbols) { unsigned int sym_demod; modem_demodulate(rx->payload_demod, mixed_sample, &sym_demod); float phase_error = modem_get_demodulator_phase_error(rx->payload_demod); - DEBUG_LOG("Sym: %d; Phase error: %f\n", sym_demod, phase_error); + //DEBUG_LOG("Sym: %d; Phase error: %f\n", sym_demod, phase_error); update_nco_pll(rx->carrier_fine_nco, phase_error); @@ -290,6 +321,7 @@ result_t layer1_rx_process(layer1_rx_t *rx, const float complex *samples, size_t assert(nsamples2dump <= (sizeof(samples2dump) / sizeof(samples2dump[0]))); dump_array_cf(samples2dump, nsamples2dump, 1.0f/(RRC_SPS * SYMBOL_RATE), "/tmp/rx_dbg.cpx64"); + dump_array_cf(samples2dump2, nsamples2dump2, 1.0f/(SYMBOL_RATE), "/tmp/rx_dbg2.cpx64"); return OK; } diff --git a/impl/src/main.c b/impl/src/main.c index 101f1c7..033bb47 100644 --- a/impl/src/main.c +++ b/impl/src/main.c @@ -218,6 +218,9 @@ int main(void) unsigned rx_retries = 0; + double old = get_hires_time(); + size_t total_samples = 0; + while(m_running) { double now = get_hires_time(); @@ -290,6 +293,13 @@ int main(void) } } + double new = get_hires_time(); + + total_samples += n_rf_samples; + double rate = total_samples / (new - old); + + fprintf(stderr, "\nEstimated rate: %.3f MS/s\n", rate / 1e6); + rx_retries = 0; fprintf(stderr, "r"); diff --git a/impl/src/results.h b/impl/src/results.h index a59d1af..fb1f821 100644 --- a/impl/src/results.h +++ b/impl/src/results.h @@ -14,7 +14,9 @@ typedef enum { #ifdef DEBUG_LIQUID #include -# define ERR_CHECK_LIQUID(call) \ +#include + +#define ERR_CHECK_LIQUID(call) \ do { \ liquid_error_code lq_result; \ if((lq_result = (call)) != LIQUID_OK) { \ diff --git a/impl/test/layer1/test_loopback.c b/impl/test/layer1/test_loopback.c index 917ce2c..739d62e 100644 --- a/impl/test/layer1/test_loopback.c +++ b/impl/test/layer1/test_loopback.c @@ -3,6 +3,7 @@ #include #include #include +#include #include @@ -107,11 +108,12 @@ int main(void) // ** Initialize ** + srand(time(NULL)); + RESULT_CHECK(layer1_tx_init(&tx)); RESULT_CHECK(layer1_rx_init(&rx, cb_rx)); - // ** Generate packets ** RESULT_CHECK(layer1_tx_reset(&tx)); @@ -139,8 +141,8 @@ int main(void) // ** Simulate channel effects ** channel_cccf ch = channel_cccf_create(); - //channel_cccf_add_awgn(ch, -30, 15); - channel_cccf_add_carrier_offset(ch, 0.01f, 1.23f); + channel_cccf_add_awgn(ch, -30, 10); + channel_cccf_add_carrier_offset(ch, 0.05f, 1.23f); float complex whole_burst_ch[burst_len]; channel_cccf_execute_block(ch, (float complex*)whole_burst, burst_len, whole_burst_ch);