From 81c02e5207b0e94ccd9f5519750d867c3df5a079 Mon Sep 17 00:00:00 2001 From: Ahmet Inan Date: Thu, 4 Jan 2024 07:23:08 +0100 Subject: [PATCH 1/7] include cstdint --- decode.cc | 1 + encode.cc | 1 + 2 files changed, 2 insertions(+) diff --git a/decode.cc b/decode.cc index c69a97f..5697c07 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 5ed87fa..264ca41 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 7d104033e170ac3c4e7bdfe38cbc71509a33a155 Mon Sep 17 00:00:00 2001 From: Ahmet Inan Date: Sat, 20 Jan 2024 10:50:42 +0100 Subject: [PATCH 2/7] added support for multiple input files --- Makefile | 2 +- README.md | 4 +- encode.cc | 110 +++++++++++++++++++++++++++--------------------------- 3 files changed, 57 insertions(+), 59 deletions(-) diff --git a/Makefile b/Makefile index afe04ed..67da107 100644 --- a/Makefile +++ b/Makefile @@ -11,7 +11,7 @@ CXX = clang++ -stdlib=libc++ -march=native all: encode decode test: encode decode - $(QEMU) ./encode encoded.wav 8000 8 1 /dev/urandom + $(QEMU) ./encode encoded.wav 8000 8 1 2000 6 ANONYMOUS /dev/urandom $(QEMU) ./decode /dev/null encoded.wav encode: encode.cc diff --git a/README.md b/README.md index c513c31..094e587 100644 --- a/README.md +++ b/README.md @@ -12,7 +12,7 @@ dd if=/dev/urandom of=uncoded.dat bs=1 count=5380 Encode file ```uncoded.dat``` to ```encoded.wav``` [WAV](https://en.wikipedia.org/wiki/WAV) audio file with 8000 Hz sample rate, 16 bits and only 1 (real) channel: ``` -./encode encoded.wav 8000 16 1 uncoded.dat +./encode encoded.wav 8000 16 1 2000 6 CALLSIGN uncoded.dat ``` Start recording to ```recorded.wav``` audio file and stop after 20 seconds: @@ -46,7 +46,7 @@ Prerequisite: [disorders](https://github.com/aicodix/disorders) Encode ```uncoded.dat``` to [analytic](https://en.wikipedia.org/wiki/Analytic_signal) audio signal, add [multipath](https://en.wikipedia.org/wiki/Multipath_propagation), [CFO, SFO](https://en.wikipedia.org/wiki/Carrier_frequency_offset), [AWGN](https://en.wikipedia.org/wiki/Additive_white_Gaussian_noise), decode and compare: ``` -./encode - 8000 16 2 uncoded.dat 2000 | ../disorders/multipath - - ../disorders/multipath.txt 10 | ../disorders/cfo - - 234.567 | ../disorders/sfo - - 147 | ../disorders/awgn - - -30 | ./decode - - | diff -s uncoded.dat - +./encode - 8000 16 2 2000 6 CALLSIGN uncoded.dat | ../disorders/multipath - - ../disorders/multipath.txt 10 | ../disorders/cfo - - 234.567 | ../disorders/sfo - - 147 | ../disorders/awgn - - -30 | ./decode - - | diff -s uncoded.dat - ``` ### Reading diff --git a/encode.cc b/encode.cc index 264ca41..07a889d 100644 --- a/encode.cc +++ b/encode.cc @@ -31,6 +31,7 @@ struct Encoder static const int symbol_len = (1280 * rate) / 8000; static const int guard_len = symbol_len / 8; static const int data_bits = 43040; + static const int data_bytes = data_bits / 8; static const int crc_bits = data_bits + 32; static const int mls0_len = 127; static const int mls0_poly = 0b10001001; @@ -267,7 +268,7 @@ struct Encoder cons_rows = cons_cnt / cons_cols; return true; } - Encoder(DSP::WritePCM *pcm, const uint8_t *inp, int freq_off, uint64_t call_sign, int oper_mode) : + Encoder(DSP::WritePCM *pcm, const uint8_t *inp, int count, int freq_off, uint64_t call_sign, int oper_mode) : pcm(pcm), crc0(0xA8F4), crc1(0xD419CC15), bchenc({ 0b100011101, 0b101110111, 0b111110011, 0b101101001, 0b110111101, 0b111100111, 0b100101011, 0b111010111, @@ -285,25 +286,27 @@ struct Encoder mls1_off = offset - mls1_len / 2; papr_min = cmplx(1000, 1000), papr_max = cmplx(-1000, -1000); pilot_block(); - schmidl_cox(); - meta_data((call_sign << 8) | oper_mode); - pilot_block(); - for (int i = 0; i < data_bits; ++i) - mesg[i] = nrz(CODE::get_le_bit(inp, i)); - crc1.reset(); - for (int i = 0; i < data_bits / 8; ++i) - crc1(inp[i]); - for (int i = 0; i < 32; ++i) - mesg[i+data_bits] = nrz((crc1()>>i)&1); - for (int i = crc_bits; i < mesg_bits; ++i) - mesg[i] = 1; - polarenc(code, mesg, frozen_bits, code_order); - shorten(); - for (int j = 0; j < cons_rows; ++j) { - for (int i = 0; i < cons_cols; ++i) - fdom[bin(i+code_off)] *= - mod_map(code+mod_bits*(cons_cols*j+i)); - symbol(); + for (int k = 0; k < count; ++k) { + schmidl_cox(); + meta_data((call_sign << 8) | oper_mode); + pilot_block(); + for (int i = 0; i < data_bits; ++i) + mesg[i] = nrz(CODE::get_le_bit(inp+k*data_bytes, i)); + crc1.reset(); + for (int i = 0; i < data_bytes; ++i) + crc1(inp[i+k*data_bytes]); + for (int i = 0; i < 32; ++i) + mesg[i+data_bits] = nrz((crc1()>>i)&1); + for (int i = crc_bits; i < mesg_bits; ++i) + mesg[i] = 1; + polarenc(code, mesg, frozen_bits, code_order); + shorten(); + for (int j = 0; j < cons_rows; ++j) { + for (int i = 0; i < cons_cols; ++i) + fdom[bin(i+code_off)] *= + mod_map(code+mod_bits*(cons_cols*j+i)); + symbol(); + } } for (int i = 0; i < symbol_len; ++i) fdom[i] = 0; @@ -333,8 +336,8 @@ long long int base37_encoder(const char *str) int main(int argc, char **argv) { - if (argc < 6 || argc > 9) { - std::cerr << "usage: " << argv[0] << " OUTPUT RATE BITS CHANNELS INPUT [OFFSET] [CALLSIGN] [MODE]" << std::endl; + if (argc < 9) { + std::cerr << "usage: " << argv[0] << " OUTPUT RATE BITS CHANNELS OFFSET CALLSIGN MODE INPUT.." << std::endl; return 1; } @@ -344,30 +347,18 @@ int main(int argc, char **argv) int output_rate = std::atoi(argv[2]); int output_bits = std::atoi(argv[3]); int output_chan = std::atoi(argv[4]); - const char *input_name = argv[5]; - if (input_name[0] == '-' && input_name[1] == 0) - input_name = "/dev/stdin"; - int freq_off = output_chan == 1 ? 2000 : 0; - if (argc >= 7) - freq_off = std::atoi(argv[6]); - - long long int call_sign = base37_encoder("ANONYMOUS"); - if (argc >= 8) - call_sign = base37_encoder(argv[7]); - - if (call_sign <= 0 || call_sign >= 129961739795077L) { - std::cerr << "Unsupported call sign." << std::endl; - return 1; - } - - int oper_mode = 6; - if (argc >= 9) - oper_mode = std::atoi(argv[8]); + int freq_off = std::atoi(argv[5]); + int oper_mode = std::atoi(argv[6]); if (oper_mode < 6 || oper_mode > 13) { std::cerr << "Unsupported operation mode." << std::endl; return 1; } + long long int call_sign = base37_encoder(argv[7]); + if (call_sign <= 0 || call_sign >= 129961739795077L) { + std::cerr << "Unsupported call sign." << std::endl; + return 1; + } int band_width; switch (oper_mode) { @@ -408,33 +399,40 @@ int main(int argc, char **argv) typedef float value; typedef DSP::Complex cmplx; - std::ifstream input_file(input_name, std::ios::binary); - if (input_file.bad()) { - std::cerr << "Couldn't open file \"" << input_name << "\" for reading." << std::endl; - return 1; - } + + const int input_count = argc - 8; const int data_len = 43040 / 8; - uint8_t *input_data = new uint8_t[data_len]; - for (int i = 0; i < data_len; ++i) - input_data[i] = input_file.get(); - CODE::Xorshift32 scrambler; - for (int i = 0; i < data_len; ++i) - input_data[i] ^= scrambler(); + uint8_t *input_data = new uint8_t[data_len*input_count]; + for (int j = 0; j < input_count; ++j) { + const char *input_name = argv[j+8]; + if (argc == 9 && input_name[0] == '-' && input_name[1] == 0) + input_name = "/dev/stdin"; + std::ifstream input_file(input_name, std::ios::binary); + if (input_file.bad()) { + std::cerr << "Couldn't open file \"" << input_name << "\" for reading." << std::endl; + return 1; + } + for (int i = 0; i < data_len; ++i) + input_data[j*data_len+i] = input_file.get(); + CODE::Xorshift32 scrambler; + for (int i = 0; i < data_len; ++i) + input_data[j*data_len+i] ^= scrambler(); + } DSP::WriteWAV output_file(output_name, output_rate, output_bits, output_chan); output_file.silence(output_rate); switch (output_rate) { case 8000: - delete new Encoder(&output_file, input_data, freq_off, call_sign, oper_mode); + delete new Encoder(&output_file, input_data, input_count, freq_off, call_sign, oper_mode); break; case 16000: - delete new Encoder(&output_file, input_data, freq_off, call_sign, oper_mode); + delete new Encoder(&output_file, input_data, input_count, freq_off, call_sign, oper_mode); break; case 44100: - delete new Encoder(&output_file, input_data, freq_off, call_sign, oper_mode); + delete new Encoder(&output_file, input_data, input_count, freq_off, call_sign, oper_mode); break; case 48000: - delete new Encoder(&output_file, input_data, freq_off, call_sign, oper_mode); + delete new Encoder(&output_file, input_data, input_count, freq_off, call_sign, oper_mode); break; default: std::cerr << "Unsupported sample rate." << std::endl; From fd43c54e567657771a7cc56542450063d393d774 Mon Sep 17 00:00:00 2001 From: Ahmet Inan Date: Sun, 28 Jan 2024 12:41:35 +0100 Subject: [PATCH 3/7] operation mode comes before call sign --- encode.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/encode.cc b/encode.cc index 07a889d..f09bd7a 100644 --- a/encode.cc +++ b/encode.cc @@ -337,7 +337,7 @@ long long int base37_encoder(const char *str) int main(int argc, char **argv) { if (argc < 9) { - std::cerr << "usage: " << argv[0] << " OUTPUT RATE BITS CHANNELS OFFSET CALLSIGN MODE INPUT.." << std::endl; + std::cerr << "usage: " << argv[0] << " OUTPUT RATE BITS CHANNELS OFFSET MODE CALLSIGN INPUT.." << std::endl; return 1; } From fdc11fcec75fe384fb05a6ba1fbcbf7c020f211d Mon Sep 17 00:00:00 2001 From: Ahmet Inan Date: Fri, 9 Feb 2024 13:48:32 +0100 Subject: [PATCH 4/7] candidates are sorted now --- decode.cc | 11 +++-------- 1 file changed, 3 insertions(+), 8 deletions(-) diff --git a/decode.cc b/decode.cc index 5697c07..3f593c4 100644 --- a/decode.cc +++ b/decode.cc @@ -523,20 +523,15 @@ struct Decoder for (int i = 0; i < cons_cnt; ++i) mod_soft(code+mod_bits*i, cons[i], precision); lengthen(); - 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 2039af90a2fac44b6593d0a380a2d3f32a713964 Mon Sep 17 00:00:00 2001 From: Ahmet Inan Date: Wed, 18 May 2022 12:05:47 +0200 Subject: [PATCH 5/7] simplified precision estimation --- decode.cc | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/decode.cc b/decode.cc index 3f593c4..c1a4337 100644 --- a/decode.cc +++ b/decode.cc @@ -513,12 +513,9 @@ struct Decoder sp += norm(hard); np += norm(error); } - value snr = DSP::decibel(sp / np); + precision = sp / np; + value snr = DSP::decibel(precision); std::cerr << "init Es/N0: " << snr << " dB" << std::endl; - // $LLR=log(\frac{p(x=+1|y)}{p(x=-1|y)})$ - // $p(x|\mu,\sigma)=\frac{1}{\sqrt{2\pi}\sigma}}e^{-\frac{(x-\mu)^2}{2\sigma^2}}$ - value sigma = std::sqrt(np / (2 * sp)); - precision = 1 / (sigma * sigma); } for (int i = 0; i < cons_cnt; ++i) mod_soft(code+mod_bits*i, cons[i], precision); From 567319214c748c26fa236291bd227d92b6886862 Mon Sep 17 00:00:00 2001 From: Ahmet Inan Date: Wed, 18 May 2022 12:24:32 +0200 Subject: [PATCH 6/7] use precision estimate from each OFDM symbol --- decode.cc | 33 ++++++++++++++++++++------------- 1 file changed, 20 insertions(+), 13 deletions(-) diff --git a/decode.cc b/decode.cc index c1a4337..3a9478e 100644 --- a/decode.cc +++ b/decode.cc @@ -502,23 +502,30 @@ struct Decoder std::cerr << "coarse sfo: " << 1000000 * sfo_rad / Const::TwoPi() << " ppm" << std::endl; std::cerr << "finer cfo: " << cfo_rad * (rate / Const::TwoPi()) << " Hz " << std::endl; } - value precision = 16; if (1) { + std::cerr << "Es/N0 (dB):"; value sp = 0, np = 0; - for (int i = 0; i < cons_cnt; ++i) { - code_type tmp[mod_max]; - mod_hard(tmp, cons[i]); - cmplx hard = mod_map(tmp); - cmplx error = cons[i] - hard; - sp += norm(hard); - np += norm(error); + for (int j = 0; j < cons_rows; ++j) { + for (int i = 0; i < cons_cols; ++i) { + code_type tmp[mod_max]; + mod_hard(tmp, cons[cons_cols*j+i]); + cmplx hard = mod_map(tmp); + cmplx error = cons[cons_cols*j+i] - hard; + sp += norm(hard); + np += norm(error); + } + value precision = sp / np; + value snr = DSP::decibel(precision); + std::cerr << " " << snr; + for (int i = 0; i < cons_cols; ++i) + mod_soft(code+mod_bits*(cons_cols*j+i), cons[cons_cols*j+i], precision); } - precision = sp / np; - value snr = DSP::decibel(precision); - std::cerr << "init Es/N0: " << snr << " dB" << std::endl; + std::cerr << std::endl; + } else { + value precision = 8; + for (int i = 0; i < cons_cnt; ++i) + mod_soft(code+mod_bits*i, cons[i], precision); } - for (int i = 0; i < cons_cnt; ++i) - mod_soft(code+mod_bits*i, cons[i], precision); lengthen(); polardec(nullptr, mesg, code, frozen_bits, code_order); systematic(); From 4d43a58b19ed96cf63cde6d922fef3e6e8899b37 Mon Sep 17 00:00:00 2001 From: Ahmet Inan Date: Wed, 18 May 2022 16:38:41 +0200 Subject: [PATCH 7/7] show number of bit flips --- decode.cc | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/decode.cc b/decode.cc index 3a9478e..cce016b 100644 --- a/decode.cc +++ b/decode.cc @@ -543,8 +543,16 @@ struct Decoder std::cerr << "payload decoding error." << std::endl; return; } - for (int i = 0; i < data_bits; ++i) - CODE::set_le_bit(out, i, mesg[i].v[best] < 0); + int flips = 0; + for (int i = 0, j = 0; i < data_bits; ++i, ++j) { + while ((frozen_bits[j / 32] >> (j % 32)) & 1) + ++j; + bool received = code[j] < 0; + bool decoded = mesg[i].v[best] < 0; + flips += received != decoded; + CODE::set_le_bit(out, i, decoded); + } + std::cerr << "bit flips: " << flips << std::endl; } };