Use new frequency+phase estimation method for preamble

impl/src/correlator.c

@@ -3,6 +3,8 @@
 #include <malloc.h>
 #include <math.h>
 
+#include "freq_est.h"
+
 #include "correlator.h"
 
 bool correlator_init(correlator_ctx_t *ctx, const float complex *search_seq, size_t nsymbols)
@@ -114,13 +116,13 @@ const float complex* correlator_get_conj_search_sequence(correlator_ctx_t *ctx)
 }
 
 
-float correlator_get_mean_phase_deviation(correlator_ctx_t *ctx)
+float correlator_get_mean_phase_deviation(correlator_ctx_t *ctx, size_t L)
 {
 	float complex mean_unrotated_symbol = 0;
 
 	size_t n = (ctx->history_ptr - ctx->search_sequence_len + 1) & ctx->buffer_size_mask;
 
-	for(size_t m = 0; m < ctx->search_sequence_len; m++) {
+	for(size_t m = ctx->search_sequence_len - L; m < ctx->search_sequence_len; m++) {
 		size_t nm = (n + m) & ctx->buffer_size_mask;
 		mean_unrotated_symbol += ctx->search_sequence[m] * ctx->input_history[nm];
 	}
@@ -130,30 +132,32 @@ float correlator_get_mean_phase_deviation(correlator_ctx_t *ctx)
 }
 
 
-/* Louise & Regiannini frequency estimator */
+float correlator_get_mean_frequency_deviation(correlator_ctx_t *ctx, size_t L, float *phase_offset)
-float correlator_get_mean_frequency_deviation(correlator_ctx_t *ctx)
 {
-	float complex z[ctx->search_sequence_len];
+	float complex z[L];
 
 	size_t n = (ctx->history_ptr - ctx->search_sequence_len + 1) & ctx->buffer_size_mask;
 
 	// remove the influence of the data from the received symbols
-	for(size_t m = 0; m < ctx->search_sequence_len; m++) {
+	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] = ctx->search_sequence[m] * ctx->input_history[nm];
+		z[m - m0] = ctx->search_sequence[m] * ctx->input_history[nm];
+		printf("%+.6f%+.6fj  -  %+.6f%+.6fj\n", crealf(ctx->search_sequence[m]), cimagf(ctx->search_sequence[m]), crealf(ctx->input_history[nm]), cimagf(ctx->input_history[nm]));
 	}
 
-
+	/* Louise & Regiannini frequency estimator */
+	/*
 	// Calculate averaged phase increments for <N> sub-sequences
-	size_t N = ctx->search_sequence_len/2;
+	size_t N = L/2;
 	float complex R_kappa[N];
 
 	for(size_t kappa = 0; kappa < N; kappa++) {
-		for(size_t k = 0; k < ctx->search_sequence_len - kappa; k++) {
+		for(size_t k = 0; k < L - kappa; k++) {
 			R_kappa[kappa] += z[k + kappa] * conj(z[k]);
 		}
 
-		R_kappa[kappa] /= ctx->search_sequence_len - kappa;
+		R_kappa[kappa] /= L - kappa;
 	}
 
 	// Calculate phase estimate (in radians/sample)
@@ -165,5 +169,8 @@ float correlator_get_mean_frequency_deviation(correlator_ctx_t *ctx)
 
 	float arg = cargf(sum_R_kappa);
 
-	return arg / ((float)M_PI * (1 + N));
+	return arg / (1 + N);
+	*/
+
+	return freq_est_data_free(z, L, phase_offset);
 }

impl/src/correlator.h

@@ -120,9 +120,10 @@ const float complex* correlator_get_conj_search_sequence(correlator_ctx_t *ctx);
  * found.
  *
  * \param ctx          The correlator context to use.
+ * \param L            Number of symbols to use (at the end of the sequence)
  * \returns            The mean phase deviation in radians.
  */
-float correlator_get_mean_phase_deviation(correlator_ctx_t *ctx);
+float correlator_get_mean_phase_deviation(correlator_ctx_t *ctx, size_t L);
 
 /*!
  * \brief Calculate mean frequency deviation between current input and search sequence.
@@ -135,9 +136,12 @@ float correlator_get_mean_phase_deviation(correlator_ctx_t *ctx);
  * symbol sequences.
  *
  * \param ctx          The correlator context to use.
+ * \param L            Number of symbols to use (at the end of the sequence)
+ * \param phase_offset The phase offset of the symbol *following* the preamble.
+ *                     May be NULL if not needed.
  * \returns            The mean frequency deviation in radians/symbol.
  */
-float correlator_get_mean_frequency_deviation(correlator_ctx_t *ctx);
+float correlator_get_mean_frequency_deviation(correlator_ctx_t *ctx, size_t L, float *phase_offset);
 
 
 #endif // CORRELATOR_H

impl/src/freq_est.c

@@ -22,8 +22,7 @@ static void linear_regression(size_t nx, const float *y, float *m, float *t)
 	}
 
 	*m = numerator / denominator;
-	//*t = mean_y - (*m) * mean_x;
+	*t = mean_y - (*m) * mean_x;
-	*t = y[0]; // This should work because x always starts at 0.
 }
 
 
@@ -52,12 +51,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)
+{
+	float complex symbols[n];
+
+	// remove the data influence from the symbols
+	for(size_t i = 0; i < n; i++) {
+		symbols[i] = conjf(ref[i]) * recv[i];
+	}
+
+	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];
 
-	// square the symbols and calculate phase
 	for(size_t i = 0; i < n; i++) {
-		phases[i] = cargf(conjf(ref[i]) * recv[i]);
+		phases[i] = cargf(symbols[i]);
 	}
 
 	liquid_unwrap_phase(phases, n);

impl/src/freq_est.h

@@ -23,8 +23,8 @@ float freq_est_unknown_bpsk(const float complex *recv, size_t n, float *final_ph
  * \brief Estimate the frequency of known symbols.
  *
  * This function removes the data influence by calculating
- * recv[k]*conj(ref[k]), calculates the phase values and then performs linear
+ * recv[k]*conj(ref[k]), and then calls \ref freq_est_data_free() with the
- * regression on those to determine the frequency.
+ * result.
  *
  * \param recv        Pointer to the array of received symbols.
  * \param ref         Pointer to the array of reference symbols (e.g. a preamble).
@@ -35,4 +35,18 @@ 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 from data-free symbols.
+ *
+ * This function calculates the phase values from the given symbols and then
+ * performs linear regression on those to determine the frequency.
+ *
+ * \param symbols     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_data_free(const float complex *symbols, size_t n, float *final_phase);
+
 #endif // FREQ_EST_H

impl/src/main.c

@@ -141,15 +141,15 @@ int main(void)
 						// Preamble found!
 						printf("Preamble found at symbol %u: %.3f > %.3f\n", i/RRC_SPS, cabsf(corr_out), 0.5f * preamble_get_symbol_count());
 
-						float mean_phase_error = correlator_get_mean_phase_deviation(&preamble_correlator);
+						float phase_offset;
-						float mean_frequency_error = correlator_get_mean_frequency_deviation(&preamble_correlator);
+						float mean_frequency_error = correlator_get_mean_frequency_deviation(&preamble_correlator, FREQ_EST_L, &phase_offset);
 
-						printf("Preamble phase deviation:     %.6f rad\n", mean_phase_error);
+						printf("Phase offset:                 %.6f rad\n", phase_offset);
 						printf("Preamble frequency deviation: %.6f rad/symbol\n", mean_frequency_error);
 
 						// adjust the frequency and phase of the NCO with the estimations from the preamble
-						nco_crcf_adjust_frequency(carrier_nco, -mean_frequency_error / RRC_SPS);
+						nco_crcf_adjust_frequency(carrier_nco, mean_frequency_error / RRC_SPS * 0.5f);
-						nco_crcf_adjust_phase(carrier_nco, mean_phase_error);
+						nco_crcf_adjust_phase(carrier_nco, phase_offset);
 
 						printf("New estimated carrier frequency: %.6f\n", nco_crcf_get_frequency(carrier_nco));
 
@@ -160,7 +160,7 @@ int main(void)
 						printf("import matplotlib.pyplot as pp\n");
 						print_complex_array("pre", preamble_get_symbols(), preamble_get_symbol_count());
 						print_complex_array("recv", input_history, preamble_get_symbol_count());
-						printf("pp.plot(recv * pre.conj())\n");
+						printf("pp.plot(np.angle(recv * pre.conj()))\n");
 						printf("pp.show()\n");
 
 						// receive and decode the header
@@ -174,7 +174,7 @@ int main(void)
 						if(((i/RRC_SPS) % FREQ_EST_L) == 0) {
 							float freq_est = freq_est_unknown_bpsk(symsync_out + symsync_out_len - FREQ_EST_L, FREQ_EST_L, NULL);
 
-							float freq_adjustment = (freq_est / RRC_SPS) * 0.5f;
+							float freq_adjustment = (freq_est / RRC_SPS) * 0.3f;
 							nco_crcf_adjust_frequency(carrier_nco, freq_adjustment);
 
 							printf("Frequency adjustment: %.6f - carrier frequency: %.6f\n", freq_adjustment, nco_crcf_get_frequency(carrier_nco));

impl/test/CMakeLists.txt

@@ -2,10 +2,13 @@ add_executable(
 	test_correlator
 	../src/correlator.c
 	../src/correlator.h
+	../src/freq_est.c
+	../src/freq_est.h
 	test_correlator.c
 )
 
 target_link_libraries(
 	test_correlator
 	m
+	liquid
 )