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,5 +2,5 @@
 
 mkdir -p build
 cd build
-cmake ..
+cmake -DCMAKE_EXPORT_COMPILE_COMMANDS=true ..
 make

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.

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

impl/src/layer1/freq_est.c

@@ -1,9 +1,11 @@
 #include <liquid/liquid.h>
+#include <math.h>
 
 #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;

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.
  *

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];
+		if(freq_est_holdoff_samples == 0) { //freq_est_history_write_idx == FREQ_EST_L) {
-		for(size_t i = 0; i < FREQ_EST_L; i++) {
+			float complex linearized_history[FREQ_EST_L];
-			linearized_history[i] = freq_est_history[(i + freq_est_history_write_idx) % 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;
 }

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

impl/src/results.h

@@ -14,7 +14,9 @@ typedef enum {
 #ifdef DEBUG_LIQUID
 #include <stdio.h>
 
-#   define ERR_CHECK_LIQUID(call) \
+#include <liquid/liquid.h>
+
+#define ERR_CHECK_LIQUID(call) \
 	do { \
 		liquid_error_code lq_result; \
 		if((lq_result = (call)) != LIQUID_OK) { \

impl/test/layer1/test_loopback.c

@@ -3,6 +3,7 @@
 #include <string.h>
 #include <math.h>
 #include <ctype.h>
+#include <time.h>
 
 #include <liquid/liquid.h>
 
@@ -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_awgn(ch, -30, 10);
-	channel_cccf_add_carrier_offset(ch, 0.01f, 1.23f);
+	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);