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