diff --git a/Makefile b/Makefile index afe04ed..7b1379a 100644 --- a/Makefile +++ b/Makefile @@ -6,12 +6,15 @@ 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 test: encode decode - $(QEMU) ./encode encoded.wav 8000 8 1 /dev/urandom + $(QEMU) ./encode encoded.wav /dev/urandom $(QEMU) ./decode /dev/null encoded.wav encode: encode.cc diff --git a/README.md b/README.md index 4de45f3..db6778e 100644 --- a/README.md +++ b/README.md @@ -3,10 +3,10 @@ Quick start: -Create file ```uncoded.dat``` with ```1360``` bits of random data: +Create file ```uncoded.dat``` with ```2048``` bits of random data: ``` -dd if=/dev/urandom of=uncoded.dat bs=1 count=170 +dd if=/dev/urandom of=uncoded.dat bs=1 count=256 ``` 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: @@ -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 -q -s uncoded.dat - +./encode - uncoded.dat | ../disorders/multipath - - ../disorders/multipath.txt 10 | ../disorders/cfo - - 234.567 | ../disorders/sfo - - 147 | ../disorders/awgn - - -30 | ./decode - - | diff -q -s uncoded.dat - ``` ### Reading diff --git a/decode.cc b/decode.cc index d7e4aee..4b19b9a 100644 --- a/decode.cc +++ b/decode.cc @@ -1,7 +1,7 @@ /* OFDM modem decoder -Copyright 2021 Ahmet Inan +Copyright 2023 Ahmet Inan */ #include @@ -10,10 +10,10 @@ Copyright 2021 Ahmet Inan #include namespace DSP { using std::abs; using std::min; using std::cos; using std::sin; } #include "bip_buffer.hh" -#include "theil_sen.hh" #include "xorshift.hh" #include "trigger.hh" #include "complex.hh" +#include "permute.hh" #include "decibel.hh" #include "blockdc.hh" #include "hilbert.hh" @@ -26,8 +26,8 @@ namespace DSP { using std::abs; using std::min; using std::cos; using std::sin; #include "fft.hh" #include "mls.hh" #include "crc.hh" -#include "osd.hh" #include "psk.hh" +#include "simplex_decoder.hh" #include "polar_tables.hh" #include "polar_helper.hh" #include "polar_encoder.hh" @@ -42,13 +42,13 @@ struct SchmidlCox DSP::FastFourierTransform fwd; DSP::FastFourierTransform bwd; DSP::SMA4 cor; - DSP::SMA4 pwr; + DSP::SMA4 pwr; DSP::SMA4 match; - DSP::Delay delay; + DSP::Delay align; DSP::SchmittTrigger threshold; DSP::FallingEdgeTrigger falling; - cmplx tmp0[symbol_len], tmp1[symbol_len], tmp2[symbol_len]; - cmplx seq[symbol_len], kern[symbol_len]; + cmplx tmp0[symbol_len], tmp1[symbol_len]; + cmplx kern[symbol_len]; cmplx cmplx_shift = 0; value timing_max = 0; value phase_max = 0; @@ -72,26 +72,27 @@ public: value cfo_rad = 0; value frac_cfo = 0; - SchmidlCox(const cmplx *sequence) : threshold(value(0.17*match_len), value(0.19*match_len)) + SchmidlCox(const cmplx *sequence) : threshold(value(0.2*match_len), value(0.3*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; + cmplx P = cor(samples[search_pos] * conj(samples[search_pos+symbol_len])); + value R = value(0.5) * pwr(norm(samples[search_pos]) + norm(samples[search_pos+symbol_len])); + value min_R = 0.00001 * symbol_len; R = std::max(R, min_R); value timing = match(norm(P) / (R * R)); - value phase = delay(arg(P)); + value phase = align(arg(P)); bool collect = threshold(timing); bool process = falling(collect); + // std::cout << timing << " " << process << " " << index_max << std::endl; + // plot "< arecord -r 8000 -c 1 -f S16_LE | ./decode /dev/null - 1 | tee /tmp/data.txt" u ($1/33) w l, "/tmp/data.txt" u ($0-$3-1+16):2 w i, 0.2, 0.3 + if (!collect && !process) return false; @@ -101,6 +102,10 @@ public: index_max = match_del; } else if (index_max < symbol_len + guard_len + match_del) { ++index_max; + } else if (process) { + index_max = 0; + timing_max = 0; + return false; } if (!process) @@ -114,20 +119,20 @@ public: index_max = 0; timing_max = 0; for (int i = 0; i < symbol_len; ++i) - tmp1[i] = samples[i+symbol_pos+symbol_len] * osc(); + tmp1[i] = samples[i+symbol_pos] * 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); + bwd(tmp1, tmp0); int shift = 0; value peak = 0; value next = 0; for (int i = 0; i < symbol_len; ++i) { - value power = norm(tmp2[i]); + value power = norm(tmp1[i]); if (power > peak) { next = peak; peak = power; @@ -139,7 +144,7 @@ public: if (peak <= next * 4) return false; - int pos_err = std::nearbyint(arg(tmp2[shift]) * symbol_len / Const::TwoPi()); + 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; @@ -151,15 +156,10 @@ public: } }; -void base37_decoder(char *str, long long int val, int len) -{ - for (int i = len-1; i >= 0; --i, val /= 37) - str[i] = " 0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ"[val%37]; -} - -template struct Decoder { + typedef float value; + typedef DSP::Complex cmplx; typedef int8_t code_type; #ifdef __AVX2__ typedef SIMD mesg_type; @@ -167,52 +167,45 @@ struct Decoder typedef SIMD mesg_type; #endif typedef DSP::Const Const; - static const int code_order = 11; - static const int mod_bits = 2; + typedef PhaseShiftKeying<2, cmplx, code_type> bpsk; + typedef PhaseShiftKeying<4, cmplx, code_type> qpsk; + static const int sample_rate = 8000; + static const int code_order = 12; static const int code_len = 1 << code_order; - static const int symbol_len = (1280 * rate) / 8000; - static const int filter_len = (((21 * rate) / 8000) & ~3) | 1; + static const int meta_len = 63; + static const int symbol_len = 256; + static const int filter_len = 33; static const int guard_len = symbol_len / 8; - static const int max_bits = 1360 + 32; - static const int cons_cols = 256; - static const int cons_rows = 4; - static const int cons_total = cons_rows * cons_cols; - static const int code_off = - cons_cols / 2; - static const int mls0_len = 127; - static const int mls0_off = - mls0_len + 1; - static const int mls0_poly = 0b10001001; - static const int mls1_len = 255; - static const int mls1_off = - mls1_len / 2; - static const int mls1_poly = 0b100101011; - static const int buffer_len = 6 * (symbol_len + guard_len); - static const int search_pos = buffer_len - 4 * (symbol_len + guard_len); + static const int data_bits = 2048; + static const int mesg_bits = data_bits + 32; + static const int subcarrier_count = 64; + static const int payload_symbols = 32; + static const int first_subcarrier = -subcarrier_count / 2; + static const int cons_total = payload_symbols * subcarrier_count; + static const int buffer_len = 2 * symbol_len + guard_len; + static const int search_pos = symbol_len; DSP::ReadPCM *pcm; DSP::FastFourierTransform fwd; DSP::BlockDC blockdc; DSP::Hilbert hilbert; DSP::BipBuffer input_hist; - DSP::TheilSenEstimator tse; - SchmidlCox correlator; - CODE::CRC crc0; - CODE::CRC crc1; - CODE::OrderedStatisticsDecoder<255, 71, 4> osddec; + SchmidlCox correlator; + CODE::CRC crc; + CODE::SimplexDecoder<6> simplex; CODE::PolarEncoder polarenc; CODE::PolarListDecoder polardec; - int8_t genmat[255*71]; - mesg_type mesg[max_bits], mess[code_len]; + CODE::ReverseFisherYatesShuffle shuffle; + mesg_type mesg[mesg_bits], mess[code_len]; code_type code[code_len]; - cmplx cons[cons_total], prev[cons_cols]; + cmplx cons[cons_total], prev[subcarrier_count]; cmplx fdom[symbol_len], tdom[symbol_len]; - value index[cons_cols], phase[cons_cols]; value cfo_rad, sfo_rad; const uint32_t *frozen_bits; int symbol_pos; - int oper_mode; - int crc_bits; static int bin(int carrier) { - return (carrier + symbol_len) % symbol_len; + return (carrier + first_subcarrier + symbol_len) % symbol_len; } static int nrz(bool bit) { @@ -227,55 +220,34 @@ struct Decoder return 0; return cons; } - const cmplx *mls0_seq() + const cmplx *sync_seq() { - CODE::MLS seq0(mls0_poly); - for (int i = 0; i < symbol_len/2; ++i) + CODE::MLS seq(0b1100111); + for (int i = 0; i < symbol_len; ++i) fdom[i] = 0; - for (int i = 0; i < mls0_len; ++i) - fdom[(i+mls0_off/2+symbol_len/2)%(symbol_len/2)] = nrz(seq0()); + for (int i = 1; i < subcarrier_count; ++i) + fdom[bin(i)] = nrz(seq()); return fdom; } void systematic() { polarenc(mess, mesg, frozen_bits, code_order); int code_bits = 1 << code_order; - for (int i = 0, j = 0; i < code_bits && j < crc_bits; ++i) + for (int i = 0, j = 0; i < code_bits && j < mesg_bits; ++i) if (!((frozen_bits[i/32] >> (i%32)) & 1)) mesg[j++] = mess[i]; } - cmplx mod_map(code_type *b) - { - return PhaseShiftKeying<4, cmplx, code_type>::map(b); - } - void mod_hard(code_type *b, cmplx c) - { - PhaseShiftKeying<4, cmplx, code_type>::hard(b, c); - } - void mod_soft(code_type *b, cmplx c, value precision) - { - PhaseShiftKeying<4, cmplx, code_type>::soft(b, c, precision); - } const cmplx *next_sample() { - cmplx tmp; - pcm->read(reinterpret_cast(&tmp), 1); - if (pcm->channels() == 1) - tmp = hilbert(blockdc(tmp.real())); - return input_hist(tmp); + value real; + pcm->read(&real, 1); + return input_hist(hilbert(blockdc(real))); } Decoder(uint8_t *out, int *len, DSP::ReadPCM *pcm, int skip_count) : - pcm(pcm), correlator(mls0_seq()), crc0(0xA8F4), crc1(0x8F6E37A0) + pcm(pcm), correlator(sync_seq()), crc(0x8F6E37A0) { - CODE::BoseChaudhuriHocquenghemGenerator<255, 71>::matrix(genmat, true, { - 0b100011101, 0b101110111, 0b111110011, 0b101101001, - 0b110111101, 0b111100111, 0b100101011, 0b111010111, - 0b000010011, 0b101100101, 0b110001011, 0b101100011, - 0b100011011, 0b100111111, 0b110001101, 0b100101101, - 0b101011111, 0b111111001, 0b111000011, 0b100111001, - 0b110101001, 0b000011111, 0b110000111, 0b110110001}); - - blockdc.samples(2*(symbol_len+guard_len)); + frozen_bits = frozen_4096_2080; + blockdc.samples(filter_len); DSP::Phasor osc; const cmplx *buf; bool okay; @@ -290,70 +262,11 @@ struct Decoder symbol_pos = correlator.symbol_pos; cfo_rad = correlator.cfo_rad; std::cerr << "symbol pos: " << symbol_pos << std::endl; - std::cerr << "coarse cfo: " << cfo_rad * (rate / Const::TwoPi()) << " Hz " << std::endl; + std::cerr << "coarse cfo: " << cfo_rad * (sample_rate / Const::TwoPi()) << " Hz " << std::endl; + for (int i = 0; i < symbol_pos; ++i) + next_sample(); osc.omega(-cfo_rad); - for (int i = 0; i < symbol_len; ++i) - tdom[i] = buf[i+symbol_pos+(symbol_len+guard_len)] * osc(); - fwd(fdom, tdom); - CODE::MLS seq1(mls1_poly); - for (int i = 0; i < mls1_len; ++i) - fdom[bin(i+mls1_off)] *= nrz(seq1()); - int8_t soft[mls1_len]; - uint8_t data[(mls1_len+7)/8]; - for (int i = 0; i < mls1_len; ++i) - soft[i] = std::min(std::max( - std::nearbyint(127 * demod_or_erase( - fdom[bin(i+mls1_off)], fdom[bin(i-1+mls1_off)]).real()), - -128), 127); - bool unique = osddec(data, soft, genmat); - if (!unique) { - std::cerr << "OSD error." << std::endl; - continue; - } - uint64_t md = 0; - for (int i = 0; i < 55; ++i) - md |= (uint64_t)CODE::get_be_bit(data, i) << i; - uint16_t cs = 0; - for (int i = 0; i < 16; ++i) - cs |= (uint16_t)CODE::get_be_bit(data, i+55) << i; - crc0.reset(); - if (crc0(md<<9) != cs) { - std::cerr << "header CRC error." << std::endl; - continue; - } - oper_mode = md & 255; - if (oper_mode && (oper_mode < 14 || oper_mode > 16)) { - std::cerr << "operation mode " << oper_mode << " unsupported." << std::endl; - continue; - } - std::cerr << "oper mode: " << oper_mode << std::endl; - if ((md>>8) == 0 || (md>>8) >= 129961739795077L) { - std::cerr << "call sign unsupported." << std::endl; - continue; - } - char call_sign[10]; - base37_decoder(call_sign, md>>8, 9); - call_sign[9] = 0; - std::cerr << "call sign: " << call_sign << std::endl; - okay = true; - } while (skip_count--); - - *len = 0; - if (!okay || !oper_mode) - return; - - for (int i = 0; i < symbol_pos+(symbol_len+guard_len); ++i) - buf = next_sample(); - for (int i = 0; i < symbol_len; ++i) - tdom[i] = buf[i] * osc(); - for (int i = 0; i < guard_len; ++i) - osc(); - fwd(fdom, tdom); - for (int i = 0; i < cons_cols; ++i) - prev[i] = fdom[bin(i+code_off)]; - std::cerr << "demod "; - for (int j = 0; j < cons_rows; ++j) { for (int i = 0; i < symbol_len+guard_len; ++i) buf = next_sample(); for (int i = 0; i < symbol_len; ++i) @@ -361,83 +274,69 @@ struct Decoder for (int i = 0; i < guard_len; ++i) osc(); fwd(fdom, tdom); - for (int i = 0; i < cons_cols; ++i) - cons[cons_cols*j+i] = demod_or_erase(fdom[bin(i+code_off)], prev[i]); - for (int i = 0; i < cons_cols; ++i) - prev[i] = fdom[bin(i+code_off)]; + for (int i = 0; i < meta_len; ++i) + cons[i] = demod_or_erase(fdom[bin(i+1)], fdom[bin(i)]); + for (int i = 0; i < subcarrier_count; ++i) + prev[i] = fdom[bin(i)]; + for (int i = 0; i < meta_len; ++i) + bpsk::soft(code+i, cons[i], 8); + CODE::MLS seq(0b1000011); + for (int i = 0; i < meta_len; ++i) + code[i] *= nrz(seq()); + int oper_mode = simplex(code); + if (oper_mode != 1) { + std::cerr << "operation mode " << oper_mode << " unsupported." << std::endl; + continue; + } + std::cerr << "oper mode: " << oper_mode << std::endl; + okay = true; + } while (skip_count--); + + if (!okay) + return; + + std::cerr << "demod "; + for (int j = 0; j < payload_symbols; ++j) { + for (int i = 0; i < symbol_len+guard_len; ++i) + buf = next_sample(); + for (int i = 0; i < symbol_len; ++i) + tdom[i] = buf[i] * osc(); + for (int i = 0; i < guard_len; ++i) + osc(); + fwd(fdom, tdom); + for (int i = 0; i < subcarrier_count; ++i) + cons[subcarrier_count*j+i] = demod_or_erase(fdom[bin(i)], prev[i]); + for (int i = 0; i < subcarrier_count; ++i) + prev[i] = fdom[bin(i)]; std::cerr << "."; } std::cerr << " done" << std::endl; - if (1) { - value sum_slope = 0, sum_yint = 0; - for (int j = 0; j < cons_rows; ++j) { - for (int i = 0; i < cons_cols; ++i) { - code_type tmp[mod_bits]; - mod_hard(tmp, cons[cons_cols*j+i]); - index[i] = i + code_off; - phase[i] = arg(cons[cons_cols*j+i] * conj(mod_map(tmp))); - } - tse.compute(index, phase, cons_cols); - //std::cerr << "Theil-Sen slope = " << tse.slope() << std::endl; - //std::cerr << "Theil-Sen yint = " << tse.yint() << std::endl; - sum_slope += tse.slope(); - sum_yint += tse.yint(); - for (int i = 0; i < cons_cols; ++i) - cons[cons_cols*j+i] *= DSP::polar(1, -tse(i+code_off)); - } - value avg_slope = sum_slope / cons_rows; - value avg_yint = sum_yint / cons_rows; - //for (int i = 0; i < cons_cnt; ++i) - // cons[i] *= DSP::polar(1, -(avg_yint+avg_slope*((i%cons_cols)+code_off))); - sfo_rad -= avg_slope * symbol_len / value(symbol_len+guard_len); - cfo_rad += avg_yint / (symbol_len+guard_len); - std::cerr << "coarse sfo: " << 1000000 * sfo_rad / Const::TwoPi() << " ppm" << std::endl; - std::cerr << "finer cfo: " << cfo_rad * (rate / Const::TwoPi()) << " Hz " << std::endl; - } if (1) { std::cerr << "Es/N0 (dB):"; value sp = 0, np = 0; - for (int j = 0; j < cons_rows; ++j) { - for (int i = 0; i < cons_cols; ++i) { - code_type tmp[mod_bits]; - mod_hard(tmp, cons[cons_cols*j+i]); - cmplx hard = mod_map(tmp); - cmplx error = cons[cons_cols*j+i] - hard; + for (int j = 0; j < payload_symbols; ++j) { + for (int i = 0; i < subcarrier_count; ++i) { + code_type tmp[2]; + qpsk::hard(tmp, cons[subcarrier_count*j+i]); + cmplx hard = qpsk::map(tmp); + cmplx error = cons[subcarrier_count*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+2*(cons_cols*j+i), cons[cons_cols*j+i], precision); + for (int i = 0; i < subcarrier_count; ++i) + qpsk::soft(code+2*(subcarrier_count*j+i), cons[subcarrier_count*j+i], precision); } std::cerr << std::endl; } else { value precision = 8; for (int i = 0; i < cons_total; ++i) - mod_soft(code+mod_bits*i, cons[i], precision); + qpsk::soft(code+2*i, cons[i], precision); } - int data_bits = 0; - switch (oper_mode) { - case 14: - data_bits = 1360; - frozen_bits = frozen_2048_1392; - break; - case 15: - data_bits = 1024; - frozen_bits = frozen_2048_1056; - break; - case 16: - data_bits = 680; - frozen_bits = frozen_2048_712; - break; - default: - return; - } - *len = data_bits / 8; - crc_bits = data_bits + 32; CODE::PolarHelper::PATH metric[mesg_type::SIZE]; + shuffle(code); polardec(metric, mesg, code, frozen_bits, code_order); systematic(); int order[mesg_type::SIZE]; @@ -446,19 +345,19 @@ struct Decoder 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); - if (crc1() == 0) { + crc.reset(); + for (int i = 0; i < mesg_bits; ++i) + crc(mesg[i].v[order[k]] < 0); + if (crc() == 0) { best = order[k]; break; } } if (best < 0) { std::cerr << "payload decoding error." << std::endl; - *len = 0; return; } + *len = data_bits / 8; int flips = 0; for (int i = 0, j = 0; i < data_bits; ++i, ++j) { while ((frozen_bits[j / 32] >> (j % 32)) & 1) @@ -479,9 +378,6 @@ int main(int argc, char **argv) return 1; } - typedef float value; - typedef DSP::Complex cmplx; - const char *output_name = argv[1]; if (output_name[0] == '-' && output_name[1] == 0) output_name = "/dev/stdout"; @@ -489,10 +385,10 @@ int main(int argc, char **argv) if (input_name[0] == '-' && input_name[1] == 0) input_name = "/dev/stdin"; - DSP::ReadWAV input_file(input_name); + DSP::ReadWAV input_file(input_name); - if (input_file.channels() < 1 || input_file.channels() > 2) { - std::cerr << "Only real or analytic signal (one or two channels) supported." << std::endl; + if (input_file.channels() != 1) { + std::cerr << "Only real signal (one channel) supported." << std::endl; return 1; } @@ -500,28 +396,17 @@ int main(int argc, char **argv) if (argc > 3) skip_count = std::atoi(argv[3]); - const int data_max = 1360 / 8; + const int data_max = 2048 / 8; uint8_t *output_data = new uint8_t[data_max]; int data_len = 0; - switch (input_file.rate()) { - case 8000: - delete new Decoder(output_data, &data_len, &input_file, skip_count); - break; - case 16000: - delete new Decoder(output_data, &data_len, &input_file, skip_count); - break; - case 44100: - delete new Decoder(output_data, &data_len, &input_file, skip_count); - break; - case 48000: - delete new Decoder(output_data, &data_len, &input_file, skip_count); - break; - default: + if (input_file.rate() != 8000) { std::cerr << "Unsupported sample rate." << std::endl; return 1; } + delete new Decoder(output_data, &data_len, &input_file, skip_count); + std::ofstream output_file(output_name, std::ios::binary | std::ios::trunc); if (output_file.bad()) { std::cerr << "Couldn't open file \"" << output_name << "\" for writing." << std::endl; diff --git a/encode.cc b/encode.cc index 046be60..fed1426 100644 --- a/encode.cc +++ b/encode.cc @@ -1,7 +1,7 @@ /* OFDM modem encoder -Copyright 2021 Ahmet Inan +Copyright 2023 Ahmet Inan */ #include @@ -9,95 +9,53 @@ Copyright 2021 Ahmet Inan #include #include "xorshift.hh" #include "complex.hh" -#include "utils.hh" +#include "permute.hh" +#include "phasor.hh" #include "bitman.hh" -#include "decibel.hh" +#include "utils.hh" #include "fft.hh" #include "wav.hh" #include "pcm.hh" #include "mls.hh" #include "crc.hh" #include "psk.hh" +#include "simplex_encoder.hh" #include "polar_tables.hh" #include "polar_helper.hh" #include "polar_encoder.hh" -#include "bose_chaudhuri_hocquenghem_encoder.hh" -template struct Encoder { + typedef float value; + typedef DSP::Complex cmplx; typedef int8_t code_type; - static const int mod_bits = 2; - static const int code_order = 11; + typedef PhaseShiftKeying<4, cmplx, code_type> qpsk; + static const int code_order = 12; static const int code_len = 1 << code_order; - static const int symbol_len = (1280 * rate) / 8000; + static const int meta_len = 63; + static const int symbol_len = 256; static const int guard_len = symbol_len / 8; - static const int max_bits = 1360 + 32; - static const int cons_cols = 256; - static const int cons_rows = 4; - static const int mls0_len = 127; - static const int mls0_poly = 0b10001001; - static const int mls1_len = 255; - static const int mls1_poly = 0b100101011; - static const int mls2_poly = 0b100101010001; + static const int data_bits = 2048; + static const int mesg_bits = data_bits + 32; + static const int subcarrier_count = 64; + static const int payload_symbols = 32; + static const int first_subcarrier = 16; DSP::WritePCM *pcm; DSP::FastFourierTransform bwd; - DSP::FastFourierTransform<4*symbol_len, cmplx, -1> fwd4; - DSP::FastFourierTransform<4*symbol_len, cmplx, 1> bwd4; - CODE::CRC crc0; - CODE::CRC crc1; - CODE::BoseChaudhuriHocquenghemEncoder<255, 71> bchenc; + CODE::CRC crc; + CODE::SimplexEncoder<6> simplex; CODE::PolarSysEnc polarenc; - code_type code[code_len], mesg[max_bits]; - cmplx fdom[symbol_len], fdom4[4*symbol_len]; - cmplx tdom[symbol_len], tdom4[4*symbol_len]; - cmplx temp[symbol_len]; - cmplx guard[guard_len]; - cmplx papr_min, papr_max; - int code_off; - int mls0_off; - int mls1_off; + CODE::FisherYatesShuffle shuffle; + code_type code[code_len], mesg[mesg_bits]; + cmplx fdom[symbol_len], tdom[symbol_len], guard[guard_len]; - static int bin(int carrier) - { - return (carrier + symbol_len) % symbol_len; - } - static int bin4(int carrier) - { - return (carrier + 4*symbol_len) % (4*symbol_len); - } static int nrz(bool bit) { return 1 - 2 * bit; } - void improve_papr() + void symbol(bool output_guard = true) { - for (int i = 0; i < 4*symbol_len; ++i) - fdom4[i] = 0; - for (int i = -symbol_len/2; i < symbol_len/2; ++i) - fdom4[bin4(i)] = fdom[bin(i)]; - bwd4(tdom4, fdom4); - for (int i = 0; i < 4*symbol_len; ++i) - tdom4[i] /= std::sqrt(value(4*symbol_len)); - for (int i = 0; i < 4*symbol_len; ++i) { - value amp = std::max(std::abs(tdom4[i].real()), std::abs(tdom4[i].imag())); - if (amp > value(1)) - tdom4[i] /= amp; - } - fwd4(fdom4, tdom4); - for (int i = -symbol_len/2; i < symbol_len/2; ++i) - if (norm(temp[bin(i)])) - temp[bin(i)] = fdom4[bin4(i)] / std::sqrt(value(4*symbol_len)); - else - temp[bin(i)] = 0; - } - void symbol(bool papr_reduction = true) - { - for (int i = 0; i < symbol_len; ++i) - temp[i] = fdom[i]; - if (papr_reduction) - improve_papr(); - bwd(tdom, temp); + bwd(tdom, fdom); for (int i = 0; i < symbol_len; ++i) tdom[i] /= std::sqrt(value(8*symbol_len)); for (int i = 0; i < guard_len; ++i) { @@ -107,239 +65,103 @@ struct Encoder x = value(0.5) * (value(1) - std::cos(DSP::Const::Pi() * x)); guard[i] = DSP::lerp(guard[i], tdom[i+symbol_len-guard_len], x); } - cmplx peak, mean; - for (int i = 0; i < symbol_len; ++i) { - cmplx power(tdom[i].real() * tdom[i].real(), tdom[i].imag() * tdom[i].imag()); - peak = cmplx(std::max(peak.real(), power.real()), std::max(peak.imag(), power.imag())); - mean += power; - } - if (mean.real() > 0 && mean.imag() > 0) { - cmplx papr(peak.real() / mean.real(), peak.imag() / mean.imag()); - papr *= value(symbol_len); - papr_min = cmplx(std::min(papr_min.real(), papr.real()), std::min(papr_min.imag(), papr.imag())); - papr_max = cmplx(std::max(papr_max.real(), papr.real()), std::max(papr_max.imag(), papr.imag())); - } - pcm->write(reinterpret_cast(guard), guard_len, 2); + if (output_guard) + pcm->write(reinterpret_cast(guard), guard_len, 2); pcm->write(reinterpret_cast(tdom), symbol_len, 2); for (int i = 0; i < guard_len; ++i) guard[i] = tdom[i]; } - void pilot_block() - { - CODE::MLS seq2(mls2_poly); - value code_fac = std::sqrt(value(symbol_len) / value(cons_cols)); + void leading_noise() { + CODE::MLS seq(0b100101010001); + value factor = std::sqrt(value(symbol_len) / value(subcarrier_count)); for (int i = 0; i < symbol_len; ++i) fdom[i] = 0; - for (int i = code_off; i < code_off + cons_cols; ++i) - fdom[bin(i)] = code_fac * nrz(seq2()); - symbol(); + for (int j = 0; j < 14; ++j) { + for (int i = 0; i < subcarrier_count; ++i) + fdom[first_subcarrier+i] = factor * cmplx(nrz(seq()), nrz(seq())); + symbol(); + } } void schmidl_cox() { - CODE::MLS seq0(mls0_poly); - value mls0_fac = std::sqrt(value(2 * symbol_len) / value(mls0_len)); + CODE::MLS seq(0b1100111); for (int i = 0; i < symbol_len; ++i) fdom[i] = 0; - fdom[bin(mls0_off-2)] = mls0_fac; - for (int i = 0; i < mls0_len; ++i) - fdom[bin(2*i+mls0_off)] = nrz(seq0()); - for (int i = 0; i < mls0_len; ++i) - fdom[bin(2*i+mls0_off)] *= fdom[bin(2*(i-1)+mls0_off)]; + fdom[first_subcarrier] = std::sqrt(value(2 * symbol_len) / value(subcarrier_count)); + for (int i = first_subcarrier + 1; i < first_subcarrier + subcarrier_count; ++i) + fdom[i] = fdom[i-1] * cmplx(nrz(seq())); + symbol(); symbol(false); } - void meta_data(uint64_t md) + void meta_data(int data) { - uint8_t data[9] = { 0 }, parity[23] = { 0 }; - for (int i = 0; i < 55; ++i) - CODE::set_be_bit(data, i, (md>>i)&1); - crc0.reset(); - uint16_t cs = crc0(md << 9); - for (int i = 0; i < 16; ++i) - CODE::set_be_bit(data, i+55, (cs>>i)&1); - bchenc(data, parity); - CODE::MLS seq4(mls1_poly); - value mls1_fac = std::sqrt(value(symbol_len) / value(mls1_len)); - for (int i = 0; i < symbol_len; ++i) - fdom[i] = 0; - fdom[bin(mls1_off-1)] = mls1_fac; - for (int i = 0; i < 71; ++i) - fdom[bin(i+mls1_off)] = nrz(CODE::get_be_bit(data, i)); - for (int i = 71; i < mls1_len; ++i) - fdom[bin(i+mls1_off)] = nrz(CODE::get_be_bit(parity, i-71)); - for (int i = 0; i < mls1_len; ++i) - fdom[bin(i+mls1_off)] *= fdom[bin(i-1+mls1_off)]; - for (int i = 0; i < mls1_len; ++i) - fdom[bin(i+mls1_off)] *= nrz(seq4()); + simplex(code, data); + CODE::MLS seq(0b1000011); + fdom[first_subcarrier] = std::sqrt(value(symbol_len) / value(subcarrier_count)); + for (int i = 0; i < meta_len; ++i) + fdom[first_subcarrier+1+i] = fdom[first_subcarrier+i] * cmplx(code[i] * nrz(seq())); symbol(); } - cmplx mod_map(code_type *b) + Encoder(DSP::WritePCM *pcm, const uint8_t *inp) : pcm(pcm), crc(0x8F6E37A0) { - return PhaseShiftKeying<4, cmplx, code_type>::map(b); - } - Encoder(DSP::WritePCM *pcm, const uint8_t *inp, int freq_off, uint64_t call_sign, int oper_mode) : - pcm(pcm), crc0(0xA8F4), crc1(0x8F6E37A0), bchenc({ - 0b100011101, 0b101110111, 0b111110011, 0b101101001, - 0b110111101, 0b111100111, 0b100101011, 0b111010111, - 0b000010011, 0b101100101, 0b110001011, 0b101100011, - 0b100011011, 0b100111111, 0b110001101, 0b100101101, - 0b101011111, 0b111111001, 0b111000011, 0b100111001, - 0b110101001, 0b000011111, 0b110000111, 0b110110001}) - { - int offset = (freq_off * symbol_len) / rate; - code_off = offset - cons_cols / 2; - mls0_off = offset - mls0_len + 1; - mls1_off = offset - mls1_len / 2; - papr_min = cmplx(1000, 1000), papr_max = cmplx(-1000, -1000); - pilot_block(); + leading_noise(); schmidl_cox(); - meta_data((call_sign << 8) | oper_mode); - if (oper_mode > 0) { - const uint32_t *frozen_bits = nullptr; - int data_bits = 0; - switch (oper_mode) { - case 14: - data_bits = 1360; - frozen_bits = frozen_2048_1392; - break; - case 15: - data_bits = 1024; - frozen_bits = frozen_2048_1056; - break; - case 16: - data_bits = 680; - frozen_bits = frozen_2048_712; - break; - default: - return; - } - 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); - polarenc(code, mesg, frozen_bits, code_order); - 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(); - } + meta_data(1); + for (int i = 0; i < data_bits; ++i) + mesg[i] = nrz(CODE::get_le_bit(inp, i)); + crc.reset(); + for (int i = 0; i < data_bits / 8; ++i) + crc(inp[i]); + for (int i = 0; i < 32; ++i) + mesg[i+data_bits] = nrz((crc()>>i)&1); + polarenc(code, mesg, frozen_4096_2080, code_order); + shuffle(code); + for (int j = 0; j < payload_symbols; ++j) { + for (int i = 0; i < subcarrier_count; ++i) + fdom[first_subcarrier+i] *= + qpsk::map(code+2*(subcarrier_count*j+i)); + symbol(); } for (int i = 0; i < symbol_len; ++i) fdom[i] = 0; symbol(); - std::cerr << "real PAPR: " << DSP::decibel(papr_min.real()) << " .. " << DSP::decibel(papr_max.real()) << " dB" << std::endl; - if (pcm->channels() == 2) - std::cerr << "imag PAPR: " << DSP::decibel(papr_min.imag()) << " .. " << DSP::decibel(papr_max.imag()) << " dB" << std::endl; } }; -long long int base37_encoder(const char *str) -{ - long long int acc = 0; - for (char c = *str++; c; c = *str++) { - acc *= 37; - if (c >= '0' && c <= '9') - acc += c - '0' + 1; - else if (c >= 'a' && c <= 'z') - acc += c - 'a' + 11; - else if (c >= 'A' && c <= 'Z') - acc += c - 'A' + 11; - else if (c != ' ') - return -1; - } - return acc; -} - int main(int argc, char **argv) { - if (argc < 6 || argc > 8) { - std::cerr << "usage: " << argv[0] << " OUTPUT RATE BITS CHANNELS INPUT [OFFSET] [CALLSIGN]" << std::endl; + if (argc != 3) { + std::cerr << "usage: " << argv[0] << " OUTPUT INPUT" << std::endl; return 1; } const char *output_name = argv[1]; if (output_name[0] == '-' && output_name[1] == 0) output_name = "/dev/stdout"; - 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]; + const char *input_name = argv[2]; 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; - } - - const int band_width = 1600; - - if ((output_chan == 1 && freq_off < band_width / 2) || freq_off < band_width / 2 - output_rate / 2 || freq_off > output_rate / 2 - band_width / 2) { - std::cerr << "Unsupported frequency offset." << std::endl; - return 1; - } - - if (freq_off % 50) { - std::cerr << "Frequency offset must be divisible by 50." << std::endl; - return 1; - } - - 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 data_len = 1360 / 8; + const int data_len = 2048 / 8; uint8_t *input_data = new uint8_t[data_len]; for (int i = 0; i < data_len; ++i) input_data[i] = std::max(input_file.get(), 0); - int oper_mode = 0; - for (int i = 128; i < 170; ++i) - if (!oper_mode && input_data[i]) - oper_mode = 14; - for (int i = 85; i < 128; ++i) - if (!oper_mode && input_data[i]) - oper_mode = 15; - for (int i = 0; i < 85; ++i) - if (!oper_mode && input_data[i]) - oper_mode = 16; + CODE::Xorshift32 scrambler; for (int i = 0; i < data_len; ++i) input_data[i] ^= scrambler(); - DSP::WriteWAV output_file(output_name, output_rate, output_bits, output_chan); + const int output_rate = 8000; + const int output_bits = 16; + const int output_chan = 1; + 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); - break; - case 16000: - delete new Encoder(&output_file, input_data, freq_off, call_sign, oper_mode); - break; - case 44100: - delete new Encoder(&output_file, input_data, freq_off, call_sign, oper_mode); - break; - case 48000: - delete new Encoder(&output_file, input_data, freq_off, call_sign, oper_mode); - break; - default: - std::cerr << "Unsupported sample rate." << std::endl; - return 1; - } + delete new Encoder(&output_file, input_data); output_file.silence(output_rate); delete []input_data; diff --git a/freezer.cc b/freezer.cc index c60e62e..b398b18 100644 --- a/freezer.cc +++ b/freezer.cc @@ -33,8 +33,6 @@ void code(int N, int K) int main() { - code<11>(2048, 1360+32); - code<11>(2048, 1024+32); - code<11>(2048, 680+32); + code<12>(4096, 2048+32); return 0; } diff --git a/polar_tables.hh b/polar_tables.hh index 56c0b90..b477b47 100644 --- a/polar_tables.hh +++ b/polar_tables.hh @@ -1,3 +1 @@ -static const uint32_t frozen_2048_1392[64] = { 0xffffffff, 0xffffffff, 0xffffffff, 0xffffffff, 0xffffffff, 0xffffffff, 0x7fffffff, 0x11f7fff, 0xffffffff, 0x7fffffff, 0x17ffffff, 0x117177f, 0x177f7fff, 0x1037f, 0x1011f, 0x1, 0xffffffff, 0x177fffff, 0x77f7fff, 0x1011f, 0x1173fff, 0x10117, 0x10117, 0x0, 0x117177f, 0x17, 0x3, 0x0, 0x1, 0x0, 0x0, 0x0, 0x7fffffff, 0x11f7fff, 0x11717ff, 0x117, 0x17177f, 0x3, 0x1, 0x0, 0x1037f, 0x1, 0x1, 0x0, 0x1, 0x0, 0x0, 0x0, 0x1011f, 0x1, 0x1, 0x0, 0x1, 0x0, 0x0, 0x0, 0x1, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, }; -static const uint32_t frozen_2048_1056[64] = { 0xffffffff, 0xffffffff, 0xffffffff, 0xffffffff, 0xffffffff, 0xffffffff, 0xffffffff, 0x7fffffff, 0xffffffff, 0xffffffff, 0xffffffff, 0x7fffffff, 0xffffffff, 0x177fffff, 0x177f7fff, 0x1017f, 0xffffffff, 0xffffffff, 0xffffffff, 0x177f7fff, 0x7fffffff, 0x13f7fff, 0x1171fff, 0x117, 0x3fffffff, 0x11717ff, 0x7177f, 0x1, 0x1017f, 0x1, 0x1, 0x0, 0xffffffff, 0x7fffffff, 0x7fffffff, 0x1171fff, 0x17ffffff, 0x7177f, 0x1037f, 0x1, 0x77f7fff, 0x1013f, 0x10117, 0x1, 0x10117, 0x0, 0x0, 0x0, 0x1173fff, 0x10117, 0x117, 0x0, 0x7, 0x0, 0x0, 0x0, 0x1, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, }; -static const uint32_t frozen_2048_712[64] = { 0xffffffff, 0xffffffff, 0xffffffff, 0xffffffff, 0xffffffff, 0xffffffff, 0xffffffff, 0xffffffff, 0xffffffff, 0xffffffff, 0xffffffff, 0xffffffff, 0xffffffff, 0xffffffff, 0xffffffff, 0x177fffff, 0xffffffff, 0xffffffff, 0xffffffff, 0xffffffff, 0xffffffff, 0xffffffff, 0x7fffffff, 0x11f7fff, 0xffffffff, 0x7fffffff, 0x1fffffff, 0x17177f, 0x177fffff, 0x1037f, 0x1011f, 0x1, 0xffffffff, 0xffffffff, 0xffffffff, 0x7fffffff, 0xffffffff, 0x1fffffff, 0x177fffff, 0x1077f, 0xffffffff, 0x177f7fff, 0x13f7fff, 0x10117, 0x1171fff, 0x117, 0x7, 0x0, 0x7fffffff, 0x1173fff, 0x11717ff, 0x7, 0x3077f, 0x1, 0x1, 0x0, 0x1013f, 0x1, 0x1, 0x0, 0x1, 0x0, 0x0, 0x0, }; +static const uint32_t frozen_4096_2080[128] = { 0xffffffff, 0xffffffff, 0xffffffff, 0xffffffff, 0xffffffff, 0xffffffff, 0xffffffff, 0xffffffff, 0xffffffff, 0xffffffff, 0xffffffff, 0xffffffff, 0xffffffff, 0xffffffff, 0xffffffff, 0x7fffffff, 0xffffffff, 0xffffffff, 0xffffffff, 0xffffffff, 0xffffffff, 0xffffffff, 0xffffffff, 0x17ffffff, 0xffffffff, 0xffffffff, 0xffffffff, 0x77f7fff, 0x7fffffff, 0x1173fff, 0x11717ff, 0x117, 0xffffffff, 0xffffffff, 0xffffffff, 0xffffffff, 0xffffffff, 0xffffffff, 0x7fffffff, 0x11f7fff, 0xffffffff, 0x7fffffff, 0x1fffffff, 0x117177f, 0x177fffff, 0x1077f, 0x1013f, 0x1, 0xffffffff, 0x177fffff, 0x177f7fff, 0x1013f, 0x11f7fff, 0x10117, 0x10117, 0x1, 0x11717ff, 0x117, 0x7, 0x0, 0x1, 0x0, 0x0, 0x0, 0xffffffff, 0xffffffff, 0xffffffff, 0x7fffffff, 0xffffffff, 0x17ffffff, 0x177fffff, 0x1077f, 0xffffffff, 0x177f7fff, 0x11f7fff, 0x10117, 0x11717ff, 0x10117, 0x17, 0x0, 0x7fffffff, 0x1171fff, 0x117177f, 0x17, 0x3177f, 0x1, 0x1, 0x0, 0x1017f, 0x1, 0x1, 0x0, 0x1, 0x0, 0x0, 0x0, 0x177fffff, 0x7177f, 0x1037f, 0x1, 0x1011f, 0x1, 0x1, 0x0, 0x10117, 0x1, 0x1, 0x0, 0x0, 0x0, 0x0, 0x0, 0x10117, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, }; diff --git a/psk.hh b/psk.hh index a07e85f..f20e3bc 100644 --- a/psk.hh +++ b/psk.hh @@ -26,7 +26,7 @@ struct PhaseShiftKeying<2, TYPE, CODE> if (std::is_integral::value) value = std::nearbyint(value); if (std::is_same::value) - value = std::min(std::max(value, -128), 127); + value = std::min(std::max(value, -127), 127); return value; } @@ -65,7 +65,7 @@ struct PhaseShiftKeying<4, TYPE, CODE> if (std::is_integral::value) value = std::nearbyint(value); if (std::is_same::value) - value = std::min(std::max(value, -128), 127); + value = std::min(std::max(value, -127), 127); return value; } @@ -111,7 +111,7 @@ struct PhaseShiftKeying<8, TYPE, CODE> if (std::is_integral::value) value = std::nearbyint(value); if (std::is_same::value) - value = std::min(std::max(value, -128), 127); + value = std::min(std::max(value, -127), 127); return value; }