From 186954fd81dd8b8098a1a01c4229c518b086e99a Mon Sep 17 00:00:00 2001 From: Ahmet Inan Date: Fri, 9 Feb 2024 13:45:07 +0100 Subject: [PATCH 1/7] candidates are sorted now --- decode.cc | 11 +++-------- 1 file changed, 3 insertions(+), 8 deletions(-) diff --git a/decode.cc b/decode.cc index d7e4aee..e9f3612 100644 --- a/decode.cc +++ b/decode.cc @@ -437,20 +437,15 @@ struct Decoder } *len = data_bits / 8; crc_bits = data_bits + 32; - CODE::PolarHelper::PATH metric[mesg_type::SIZE]; - polardec(metric, mesg, code, frozen_bits, code_order); + polardec(nullptr, 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; } } From 3112e07854ecd01b483dc20e79fb871c4a969591 Mon Sep 17 00:00:00 2001 From: Ahmet Inan Date: Wed, 21 Feb 2024 11:28:35 +0100 Subject: [PATCH 2/7] avoid oversaturation decoder errors at high SNR --- decode.cc | 2 ++ 1 file changed, 2 insertions(+) diff --git a/decode.cc b/decode.cc index e9f3612..56d98fd 100644 --- a/decode.cc +++ b/decode.cc @@ -409,6 +409,8 @@ struct Decoder value precision = sp / np; value snr = DSP::decibel(precision); std::cerr << " " << snr; + if (std::is_same::value && precision > 32) + precision = 32; for (int i = 0; i < cons_cols; ++i) mod_soft(code+2*(cons_cols*j+i), cons[cons_cols*j+i], precision); } From 459ef628c5d86ccd39d1518e363a77f2a0ade905 Mon Sep 17 00:00:00 2001 From: Ahmet Inan Date: Fri, 1 Mar 2024 09:06:00 +0100 Subject: [PATCH 3/7] made gcc happy --- decode.cc | 1 + encode.cc | 1 + 2 files changed, 2 insertions(+) diff --git a/decode.cc b/decode.cc index 56d98fd..63d38bc 100644 --- a/decode.cc +++ b/decode.cc @@ -7,6 +7,7 @@ Copyright 2021 Ahmet Inan #include #include #include +#include #include namespace DSP { using std::abs; using std::min; using std::cos; using std::sin; } #include "bip_buffer.hh" diff --git a/encode.cc b/encode.cc index 7ffc42f..98afa70 100644 --- a/encode.cc +++ b/encode.cc @@ -6,6 +6,7 @@ Copyright 2021 Ahmet Inan #include #include +#include #include #include "xorshift.hh" #include "complex.hh" From 4d64b30affd1100afab80b72f9245cc1ae643c7a Mon Sep 17 00:00:00 2001 From: Ahmet Inan Date: Fri, 1 Mar 2024 09:14:49 +0100 Subject: [PATCH 4/7] lower factor even further to improve corner cases --- decode.cc | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/decode.cc b/decode.cc index 63d38bc..c370a12 100644 --- a/decode.cc +++ b/decode.cc @@ -410,8 +410,8 @@ struct Decoder value precision = sp / np; value snr = DSP::decibel(precision); std::cerr << " " << snr; - if (std::is_same::value && precision > 32) - precision = 32; + if (std::is_same::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); } From 8fb7a2aeb896501a472c75a40f9a2e6f898e695c Mon Sep 17 00:00:00 2001 From: Ahmet Inan Date: Thu, 11 May 2023 09:06:40 +0200 Subject: [PATCH 5/7] added aarch64 example --- Makefile | 3 +++ 1 file changed, 3 insertions(+) diff --git a/Makefile b/Makefile index afe04ed..0bb0366 100644 --- a/Makefile +++ b/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 From 02de94951f73d922a6d3e249cec5fe89ec373373 Mon Sep 17 00:00:00 2001 From: Ahmet Inan Date: Sat, 9 Mar 2024 19:29:46 +0100 Subject: [PATCH 6/7] algorithm header not needed anymore --- decode.cc | 1 - 1 file changed, 1 deletion(-) diff --git a/decode.cc b/decode.cc index c370a12..41e729b 100644 --- a/decode.cc +++ b/decode.cc @@ -4,7 +4,6 @@ OFDM modem decoder Copyright 2021 Ahmet Inan */ -#include #include #include #include From 4f5511f78f70e1ccf7e060efba2addcab5dc56c4 Mon Sep 17 00:00:00 2001 From: Ahmet Inan Date: Sun, 17 Mar 2024 12:04:03 +0100 Subject: [PATCH 7/7] moved Schmidl & Cox correlator to own header --- decode.cc | 121 +--------------------------------------------- schmidl_cox.hh | 129 +++++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 130 insertions(+), 120 deletions(-) create mode 100644 schmidl_cox.hh diff --git a/decode.cc b/decode.cc index 41e729b..4582457 100644 --- a/decode.cc +++ b/decode.cc @@ -9,10 +9,10 @@ Copyright 2021 Ahmet Inan #include #include 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 -struct SchmidlCox -{ - typedef DSP::Const Const; - static const int match_len = guard_len | 1; - static const int match_del = (match_len - 1) / 2; - DSP::FastFourierTransform fwd; - DSP::FastFourierTransform bwd; - DSP::SMA4 cor; - DSP::SMA4 pwr; - DSP::SMA4 match; - DSP::Delay delay; - DSP::SchmittTrigger 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 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) diff --git a/schmidl_cox.hh b/schmidl_cox.hh new file mode 100644 index 0000000..16d07f1 --- /dev/null +++ b/schmidl_cox.hh @@ -0,0 +1,129 @@ +/* +Schmidl & Cox correlator + +Copyright 2021 Ahmet Inan +*/ + +#pragma once + +#include "fft.hh" +#include "sma.hh" +#include "phasor.hh" +#include "trigger.hh" + +template +class SchmidlCox { + typedef DSP::Const Const; + static const int match_len = guard_len | 1; + static const int match_del = (match_len - 1) / 2; + DSP::FastFourierTransform fwd; + DSP::FastFourierTransform bwd; + DSP::SMA4 cor; + DSP::SMA4 pwr; + DSP::SMA4 match; + DSP::Delay align; + DSP::SchmittTrigger 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 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; + } +};