moved Schmidl & Cox correlator to own header

Makefile

@@ -6,6 +6,9 @@ CXX = clang++ -stdlib=libc++ -march=native
 #CXX = armv7a-hardfloat-linux-gnueabi-g++ -static -mfpu=neon -march=armv7-a
 #QEMU = qemu-arm
 
+#CXX = aarch64-unknown-linux-gnu-g++ -static -march=armv8-a+crc+simd -mtune=cortex-a72
+#QEMU = qemu-aarch64
+
 .PHONY: all
 
 all: encode decode

decode.cc

@@ -4,15 +4,15 @@ OFDM modem decoder
 Copyright 2021 Ahmet Inan <inan@aicodix.de>
 */
 
-#include <algorithm>
 #include <iostream>
 #include <cassert>
+#include <cstdint>
 #include <cmath>
 namespace DSP { using std::abs; using std::min; using std::cos; using std::sin; }
+#include "schmidl_cox.hh"
 #include "bip_buffer.hh"
 #include "theil_sen.hh"
 #include "xorshift.hh"
-#include "trigger.hh"
 #include "complex.hh"
 #include "decibel.hh"
 #include "blockdc.hh"
@@ -20,7 +20,6 @@ namespace DSP { using std::abs; using std::min; using std::cos; using std::sin;
 #include "phasor.hh"
 #include "bitman.hh"
 #include "delay.hh"
-#include "sma.hh"
 #include "wav.hh"
 #include "pcm.hh"
 #include "fft.hh"
@@ -33,124 +32,6 @@ namespace DSP { using std::abs; using std::min; using std::cos; using std::sin;
 #include "polar_encoder.hh"
 #include "polar_list_decoder.hh"
 
-template <typename value, typename cmplx, int search_pos, int symbol_len, int guard_len>
-struct SchmidlCox
-{
-	typedef DSP::Const<value> Const;
-	static const int match_len = guard_len | 1;
-	static const int match_del = (match_len - 1) / 2;
-	DSP::FastFourierTransform<symbol_len, cmplx, -1> fwd;
-	DSP::FastFourierTransform<symbol_len, cmplx, 1> bwd;
-	DSP::SMA4<cmplx, value, symbol_len, false> cor;
-	DSP::SMA4<value, value, 2*symbol_len, false> pwr;
-	DSP::SMA4<value, value, match_len, false> match;
-	DSP::Delay<value, match_del> delay;
-	DSP::SchmittTrigger<value> threshold;
-	DSP::FallingEdgeTrigger falling;
-	cmplx tmp0[symbol_len], tmp1[symbol_len], tmp2[symbol_len];
-	cmplx seq[symbol_len], kern[symbol_len];
-	cmplx cmplx_shift = 0;
-	value timing_max = 0;
-	value phase_max = 0;
-	int index_max = 0;
-
-	static int bin(int carrier)
-	{
-		return (carrier + symbol_len) % symbol_len;
-	}
-	static cmplx demod_or_erase(cmplx curr, cmplx prev)
-	{
-		if (!(norm(prev) > 0))
-			return 0;
-		cmplx cons = curr / prev;
-		if (!(norm(cons) <= 4))
-			return 0;
-		return cons;
-	}
-public:
-	int symbol_pos = 0;
-	value cfo_rad = 0;
-	value frac_cfo = 0;
-
-	SchmidlCox(const cmplx *sequence) : threshold(value(0.17*match_len), value(0.19*match_len))
-	{
-		for (int i = 0; i < symbol_len; ++i)
-			seq[i] = sequence[i];
-		fwd(kern, sequence);
-		for (int i = 0; i < symbol_len; ++i)
-			kern[i] = conj(kern[i]) / value(symbol_len);
-	}
-	bool operator()(const cmplx *samples)
-	{
-		cmplx P = cor(samples[search_pos+symbol_len] * conj(samples[search_pos+2*symbol_len]));
-		value R = value(0.5) * pwr(norm(samples[search_pos+2*symbol_len]));
-		value min_R = 0.0001 * symbol_len;
-		R = std::max(R, min_R);
-		value timing = match(norm(P) / (R * R));
-		value phase = delay(arg(P));
-
-		bool collect = threshold(timing);
-		bool process = falling(collect);
-
-		if (!collect && !process)
-			return false;
-
-		if (timing_max < timing) {
-			timing_max = timing;
-			phase_max = phase;
-			index_max = match_del;
-		} else if (index_max < symbol_len + guard_len + match_del) {
-			++index_max;
-		}
-
-		if (!process)
-			return false;
-
-		frac_cfo = phase_max / value(symbol_len);
-
-		DSP::Phasor<cmplx> osc;
-		osc.omega(frac_cfo);
-		symbol_pos = search_pos - index_max;
-		index_max = 0;
-		timing_max = 0;
-		for (int i = 0; i < symbol_len; ++i)
-			tmp1[i] = samples[i+symbol_pos+symbol_len] * osc();
-		fwd(tmp0, tmp1);
-		for (int i = 0; i < symbol_len; ++i)
-			tmp1[i] = demod_or_erase(tmp0[i], tmp0[bin(i-1)]);
-		fwd(tmp0, tmp1);
-		for (int i = 0; i < symbol_len; ++i)
-			tmp0[i] *= kern[i];
-		bwd(tmp2, tmp0);
-
-		int shift = 0;
-		value peak = 0;
-		value next = 0;
-		for (int i = 0; i < symbol_len; ++i) {
-			value power = norm(tmp2[i]);
-			if (power > peak) {
-				next = peak;
-				peak = power;
-				shift = i;
-			} else if (power > next) {
-				next = power;
-			}
-		}
-		if (peak <= next * 4)
-			return false;
-
-		int pos_err = std::nearbyint(arg(tmp2[shift]) * symbol_len / Const::TwoPi());
-		if (abs(pos_err) > guard_len / 2)
-			return false;
-		symbol_pos -= pos_err;
-
-		cfo_rad = shift * (Const::TwoPi() / symbol_len) - frac_cfo;
-		if (cfo_rad >= Const::Pi())
-			cfo_rad -= Const::TwoPi();
-		return true;
-	}
-};
-
 void base37_decoder(char *str, long long int val, int len)
 {
 	for (int i = len-1; i >= 0; --i, val /= 37)
@@ -409,6 +290,8 @@ struct Decoder
 				value precision = sp / np;
 				value snr = DSP::decibel(precision);
 				std::cerr << " " << snr;
+				if (std::is_same<code_type, int8_t>::value && precision > 8)
+					precision = 8;
 				for (int i = 0; i < cons_cols; ++i)
 					mod_soft(code+2*(cons_cols*j+i), cons[cons_cols*j+i], precision);
 			}
@@ -437,20 +320,15 @@ struct Decoder
 		}
 		*len = data_bits / 8;
 		crc_bits = data_bits + 32;
-		CODE::PolarHelper<mesg_type>::PATH metric[mesg_type::SIZE];
+		polardec(nullptr, mesg, code, frozen_bits, code_order);
-		polardec(metric, mesg, code, frozen_bits, code_order);
 		systematic();
-		int order[mesg_type::SIZE];
-		for (int k = 0; k < mesg_type::SIZE; ++k)
-			order[k] = k;
-		std::sort(order, order+mesg_type::SIZE, [metric](int a, int b){ return metric[a] < metric[b]; });
 		int best = -1;
 		for (int k = 0; k < mesg_type::SIZE; ++k) {
 			crc1.reset();
 			for (int i = 0; i < crc_bits; ++i)
-				crc1(mesg[i].v[order[k]] < 0);
+				crc1(mesg[i].v[k] < 0);
 			if (crc1() == 0) {
-				best = order[k];
+				best = k;
 				break;
 			}
 		}

encode.cc

@@ -6,6 +6,7 @@ Copyright 2021 Ahmet Inan <inan@aicodix.de>
 
 #include <iostream>
 #include <cassert>
+#include <cstdint>
 #include <cmath>
 #include "xorshift.hh"
 #include "complex.hh"

schmidl_cox.hh

@@ -0,0 +1,129 @@
+/*
+Schmidl & Cox correlator
+
+Copyright 2021 Ahmet Inan <inan@aicodix.de>
+*/
+
+#pragma once
+
+#include "fft.hh"
+#include "sma.hh"
+#include "phasor.hh"
+#include "trigger.hh"
+
+template<typename value, typename cmplx, int search_pos, int symbol_len, int guard_len>
+class SchmidlCox {
+	typedef DSP::Const<value> Const;
+	static const int match_len = guard_len | 1;
+	static const int match_del = (match_len - 1) / 2;
+	DSP::FastFourierTransform<symbol_len, cmplx, -1> fwd;
+	DSP::FastFourierTransform<symbol_len, cmplx, 1> bwd;
+	DSP::SMA4<cmplx, value, symbol_len, false> cor;
+	DSP::SMA4<value, value, 2 * symbol_len, false> pwr;
+	DSP::SMA4<value, value, match_len, false> match;
+	DSP::Delay<value, match_del> align;
+	DSP::SchmittTrigger<value> threshold;
+	DSP::FallingEdgeTrigger falling;
+	cmplx tmp0[symbol_len], tmp1[symbol_len];
+	cmplx kern[symbol_len];
+	value timing_max = 0;
+	value phase_max = 0;
+	int index_max = 0;
+
+	static int bin(int carrier) {
+		return (carrier + symbol_len) % symbol_len;
+	}
+
+	static cmplx demod_or_erase(cmplx curr, cmplx prev, value pwr) {
+		if (norm(curr) > pwr && norm(prev) > pwr) {
+			cmplx cons = curr / prev;
+			if (norm(cons) < value(4))
+				return cons;
+		}
+		return 0;
+	}
+
+public:
+	int symbol_pos = 0;
+	value cfo_rad = 0;
+	value frac_cfo = 0;
+
+	SchmidlCox(const cmplx *sequence) : threshold(value(0.17 * match_len), value(0.19 * match_len)) {
+		fwd(kern, sequence);
+		for (int i = 0; i < symbol_len; ++i)
+			kern[i] = conj(kern[i]) / value(symbol_len);
+	}
+
+	bool operator()(const cmplx *samples) {
+		cmplx P = cor(samples[search_pos + symbol_len] * conj(samples[search_pos + 2 * symbol_len]));
+		value R = value(0.5) * pwr(norm(samples[search_pos + 2 * symbol_len]));
+		value min_R = 0.00001 * symbol_len;
+		R = std::max(R, min_R);
+		value timing = match(norm(P) / (R * R));
+		value phase = align(arg(P));
+
+		bool collect = threshold(timing);
+		bool process = falling(collect);
+
+		if (!collect && !process)
+			return false;
+
+		if (timing_max < timing) {
+			timing_max = timing;
+			phase_max = phase;
+			index_max = match_del;
+		} else if (index_max < symbol_len + guard_len + match_del) {
+			++index_max;
+		}
+
+		if (!process)
+			return false;
+
+		frac_cfo = phase_max / value(symbol_len);
+
+		DSP::Phasor<cmplx> osc;
+		osc.omega(frac_cfo);
+		symbol_pos = search_pos - index_max;
+		index_max = 0;
+		timing_max = 0;
+		for (int i = 0; i < symbol_len; ++i)
+			tmp1[i] = samples[i + symbol_pos + symbol_len] * osc();
+		fwd(tmp0, tmp1);
+		value min_pwr = 0;
+		for (int i = 0; i < symbol_len; ++i)
+			min_pwr += norm(tmp0[i]);
+		min_pwr /= symbol_len;
+		for (int i = 0; i < symbol_len; ++i)
+			tmp1[i] = demod_or_erase(tmp0[i], tmp0[bin(i - 1)], min_pwr);
+		fwd(tmp0, tmp1);
+		for (int i = 0; i < symbol_len; ++i)
+			tmp0[i] *= kern[i];
+		bwd(tmp1, tmp0);
+
+		int shift = 0;
+		value peak = 0;
+		value next = 0;
+		for (int i = 0; i < symbol_len; ++i) {
+			value power = norm(tmp1[i]);
+			if (power > peak) {
+				next = peak;
+				peak = power;
+				shift = i;
+			} else if (power > next) {
+				next = power;
+			}
+		}
+		if (peak <= next * 4)
+			return false;
+
+		int pos_err = std::nearbyint(arg(tmp1[shift]) * symbol_len / Const::TwoPi());
+		if (abs(pos_err) > guard_len / 2)
+			return false;
+		symbol_pos -= pos_err;
+
+		cfo_rad = shift * (Const::TwoPi() / symbol_len) - frac_cfo;
+		if (cfo_rad >= Const::Pi())
+			cfo_rad -= Const::TwoPi();
+		return true;
+	}
+};