use the first 32 comb pilots to signal the mode

This commit is contained in:
Ahmet Inan 2025-06-30 12:36:35 +02:00
commit 4b99f82aac
2 changed files with 99 additions and 215 deletions

169
decode.cc
View file

@ -24,17 +24,12 @@ 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 "qam.hh"
#include "polar_tables.hh"
#include "polar_list_decoder.hh"
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];
}
#include "hadamard_encoder.hh"
#include "hadamard_decoder.hh"
template <typename value, typename cmplx, int rate>
struct Decoder
@ -61,8 +56,6 @@ struct Decoder
static const int mls0_poly = 0b1100110001;
static const int mls0_seed = 214;
static const int mls0_off = - mls0_len / 2;
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 = 5 * extended_len;
static const int search_pos = extended_len;
@ -73,15 +66,15 @@ struct Decoder
DSP::BipBuffer<cmplx, buffer_len> input_hist;
DSP::TheilSenEstimator<value, cols_max> tse;
SchmidlCox<value, cmplx, search_pos, symbol_len, guard_len> correlator;
CODE::CRC<uint16_t> crc0;
CODE::CRC<uint32_t> crc1;
CODE::OrderedStatisticsDecoder<255, 71, 4> osddec;
CODE::HadamardEncoder<6> hadamardenc;
CODE::HadamardDecoder<6> hadamarddec;
CODE::PolarListDecoder<mesg_type, code_max> polardec;
uint8_t output_data[data_max];
int8_t genmat[255*71];
mesg_type mesg[bits_max];
code_type code[bits_max], perm[bits_max];
cmplx cons[cons_max], prev[cols_max];
int8_t mode[32];
cmplx cons[cons_max], chan[cols_max];
cmplx fdom[symbol_len], tdom[symbol_len];
value index[cols_max], phase[cols_max];
value cfo_rad, sfo_rad;
@ -196,20 +189,18 @@ struct Decoder
return input_hist(tmp);
}
Decoder(DSP::ReadPCM<value> *pcm, const char *const *output_names, int output_count) :
pcm(pcm), correlator(mls0_seq()), crc0(0xA8F4), crc1(0x8F6E37A0)
pcm(pcm), correlator(mls0_seq()), crc1(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(filter_len);
DSP::Phasor<cmplx> osc;
const cmplx *buf;
int output_index = 0;
int comb_cols = 32;
int code_cols = 256;
int comb_dist = 9;
int comb_off = 4;
int cons_cols = code_cols + comb_cols;
int code_off = - cons_cols / 2;
while (output_index < output_count) {
do {
if (!pcm->good())
@ -223,108 +214,87 @@ struct Decoder
std::cerr << "coarse cfo: " << cfo_rad * (rate / Const::TwoPi()) << " Hz " << std::endl;
osc.omega(-cfo_rad);
for (int i = 0; i < symbol_len; ++i)
tdom[i] = buf[i+symbol_pos+symbol_len] * osc();
for (int i = 0; i < guard_len; ++i)
osc();
fwd(fdom, tdom);
CODE::MLS seq0(mls0_poly, mls0_seed);
for (int i = 0; i < mls0_len; ++i)
fdom[bin(i+mls0_off)] *= nrz(seq0());
for (int i = 0; i < cons_cols; ++i)
chan[i] = fdom[bin(i+code_off)];
for (int i = 0; i < symbol_len; ++i)
tdom[i] = buf[i+symbol_pos+symbol_len+extended_len] * osc();
for (int i = 0; i < guard_len; ++i)
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<value>(std::max<value>(
std::nearbyint(127 * demod_or_erase(
fdom[bin(i+mls1_off)], fdom[bin(i-1+mls1_off)]).real()),
-127), 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 < 23 || oper_mode > 30)) {
auto clamp = [](int v){ return v < -127 ? -127 : v > 127 ? 127 : v; };
for (int i = 0; i < comb_cols; ++i)
mode[i] = clamp(std::nearbyint(127 * demod_or_erase(fdom[bin(i*comb_dist+comb_off+code_off)], chan[i]).real() * nrz(seq1())));
oper_mode = hadamarddec(mode);
if (oper_mode < 0 || oper_mode > 8) {
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;
if (!oper_mode)
continue;
int data_bits = 0;
int cons_rows = 0;
int comb_cols = 32;
int code_cols = 256;
int comb_dist = 9;
int comb_off = 4;
switch (oper_mode) {
case 23:
case 1:
mod_bits = 2;
cons_rows = 8;
code_order = 12;
data_bits = 2048;
frozen_bits = frozen_4096_2080;
break;
case 24:
case 2:
mod_bits = 2;
cons_rows = 16;
code_order = 13;
data_bits = 4096;
frozen_bits = frozen_8192_4128;
break;
case 25:
case 3:
mod_bits = 2;
cons_rows = 32;
code_order = 14;
data_bits = 8192;
frozen_bits = frozen_16384_8224;
break;
case 26:
case 4:
mod_bits = 4;
cons_rows = 4;
code_order = 12;
data_bits = 2048;
frozen_bits = frozen_4096_2080;
break;
case 27:
case 5:
mod_bits = 4;
cons_rows = 8;
code_order = 13;
data_bits = 4096;
frozen_bits = frozen_8192_4128;
break;
case 28:
case 6:
mod_bits = 4;
cons_rows = 16;
code_order = 14;
data_bits = 8192;
frozen_bits = frozen_16384_8224;
break;
case 29:
case 7:
mod_bits = 6;
cons_rows = 6;
code_order = 13;
data_bits = 4096;
frozen_bits = frozen_8192_4128;
break;
case 30:
case 8:
mod_bits = 6;
cons_rows = 11;
code_order = 14;
@ -335,32 +305,27 @@ struct Decoder
return;
}
int data_bytes = data_bits / 8;
int cons_cols = code_cols + comb_cols;
int code_off = - cons_cols / 2;
for (int i = 0; i < symbol_pos+symbol_len+extended_len; ++i)
correlator(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 ";
CODE::MLS seq0(mls0_poly);
for (int j = 0; j < cons_rows; ++j) {
for (int i = 0; i < extended_len; ++i)
correlator(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)
cons[cons_cols*j+i] = demod_or_erase(fdom[bin(i+code_off)], prev[i]);
if (j) {
for (int i = 0; i < extended_len; ++i)
correlator(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);
} else {
for (int i = 0; i < symbol_pos+symbol_len+extended_len; ++i)
correlator(buf = next_sample());
seq1.reset();
hadamardenc(mode, oper_mode);
}
for (int i = 0; i < comb_cols; ++i)
cons[cons_cols*j+comb_dist*i+comb_off] *= nrz(seq0());
fdom[bin(comb_dist*i+comb_off+code_off)] *= nrz(seq1()) * mode[i];
for (int i = 0; i < cons_cols; ++i)
cons[cons_cols*j+i] = demod_or_erase(fdom[bin(i+code_off)], chan[i]);
for (int i = 0; i < comb_cols; ++i) {
index[i] = code_off + comb_dist * i + comb_off;
phase[i] = arg(cons[cons_cols*j+comb_dist*i+comb_off]);
@ -372,27 +337,9 @@ struct Decoder
cons[cons_cols*j+i] *= DSP::polar<value>(1, -tse(i+code_off));
for (int i = 0; i < cons_cols; ++i)
if (i % comb_dist == comb_off)
prev[i] = fdom[bin(i+code_off)];
chan[i] = fdom[bin(i+code_off)];
else
prev[i] *= DSP::polar<value>(1, tse(i+code_off));
for (int i = 0; i < cons_cols; ++i) {
index[i] = code_off + i;
if (i % comb_dist == comb_off) {
phase[i] = arg(cons[cons_cols*j+i]);
} else {
code_type tmp[mod_max];
mod_hard(tmp, cons[cons_cols*j+i]);
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;
for (int i = 0; i < cons_cols; ++i)
cons[cons_cols*j+i] *= DSP::polar<value>(1, -tse(i+code_off));
for (int i = 0; i < cons_cols; ++i)
if (i % comb_dist != comb_off)
prev[i] *= DSP::polar<value>(1, tse(i+code_off));
chan[i] *= DSP::polar<value>(1, tse(i+code_off));
std::cerr << ".";
}
std::cerr << " done" << std::endl;

145
encode.cc
View file

@ -21,7 +21,7 @@ Copyright 2021 Ahmet Inan <inan@aicodix.de>
#include "qam.hh"
#include "polar_tables.hh"
#include "polar_encoder.hh"
#include "bose_chaudhuri_hocquenghem_encoder.hh"
#include "hadamard_encoder.hh"
template <typename value, typename cmplx, int rate>
struct Encoder
@ -35,24 +35,22 @@ struct Encoder
static const int mls0_len = 320;
static const int mls0_poly = 0b1100110001;
static const int mls0_seed = 214;
static const int mls1_len = 255;
static const int mls1_poly = 0b100101011;
static const int mls2_poly = 0b100101010001;
DSP::WritePCM<value> *pcm;
DSP::FastFourierTransform<symbol_len, cmplx, -1> fwd;
DSP::FastFourierTransform<symbol_len, cmplx, 1> bwd;
CODE::CRC<uint16_t> crc0;
CODE::CRC<uint32_t> crc1;
CODE::BoseChaudhuriHocquenghemEncoder<255, 71> bchenc;
CODE::HadamardEncoder<6> hadamardenc;
CODE::PolarEncoder<code_type> polarenc;
uint8_t input_data[data_max];
code_type code[bits_max], perm[bits_max], mesg[bits_max];
int8_t mode[32];
cmplx fdom[symbol_len];
cmplx tdom[symbol_len];
cmplx temp[symbol_len];
cmplx kern[symbol_len];
cmplx guard[guard_len];
cmplx prev[cols_max];
value weight[guard_len];
value papr_min, papr_max;
int mod_bits;
@ -62,7 +60,6 @@ struct Encoder
int cons_cols;
int cons_rows;
int mls0_off;
int mls1_off;
static int bin(int carrier)
{
@ -113,7 +110,7 @@ struct Encoder
void symbol(bool papr_reduction = true, bool guard_interval = true)
{
bwd(tdom, fdom);
value scale = 2;
value scale = 2 / std::sqrt(value(symbol_len) / value(cons_cols));
for (int i = 0; i < symbol_len; ++i)
tdom[i] /= scale * std::sqrt(value(symbol_len));
clipping_and_filtering(scale, oper_mode > 0 && papr_reduction);
@ -145,56 +142,22 @@ struct Encoder
void pilot_block()
{
CODE::MLS seq2(mls2_poly);
value code_fac = std::sqrt(value(symbol_len) / value(cons_cols));
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());
fdom[bin(i)] = nrz(seq2());
symbol();
}
void schmidl_cox()
{
CODE::MLS seq0(mls0_poly, mls0_seed);
value mls0_fac = std::sqrt(value(symbol_len) / value(mls0_len));
for (int i = 0; i < symbol_len; ++i)
fdom[i] = 0;
for (int i = 0; i < mls0_len; ++i)
fdom[bin(i+mls0_off)] = mls0_fac * nrz(seq0());
fdom[bin(i+mls0_off)] = nrz(seq0());
symbol(false);
symbol(false, false);
}
void meta_data(uint64_t md)
{
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 seq1(mls1_poly);
value cons_fac = std::sqrt(value(symbol_len) / value(cons_cols));
for (int i = 0; i < symbol_len; ++i)
fdom[i] = 0;
fdom[bin(mls1_off-1)] = cons_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(seq1());
if (oper_mode > 0) {
for (int i = code_off; i < code_off + cons_cols; ++i) {
if (i == mls1_off-1)
i += mls1_len + 1;
fdom[bin(i)] = cons_fac * nrz(seq1());
}
}
symbol();
}
cmplx mod_map(code_type *b)
{
switch (mod_bits) {
@ -253,15 +216,8 @@ struct Encoder
dest[i] = src[seq()];
}
}
Encoder(DSP::WritePCM<value> *pcm, const char *const *input_names, int input_count, 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}),
oper_mode(oper_mode)
Encoder(DSP::WritePCM<value> *pcm, const char *const *input_names, int input_count, int freq_off, int oper_mode) :
pcm(pcm), crc1(0x8F6E37A0), oper_mode(oper_mode)
{
const uint32_t *frozen_bits = nullptr;
int code_cols = 256;
@ -272,59 +228,59 @@ struct Encoder
int reserved_tones = 32;
switch (oper_mode) {
case 0:
comb_cols = 0;
cons_rows = 1;
reserved_tones = 0;
break;
case 23:
case 1:
mod_bits = 2;
cons_rows = 8;
code_order = 12;
data_bits = 2048;
frozen_bits = frozen_4096_2080;
break;
case 24:
case 2:
mod_bits = 2;
cons_rows = 16;
code_order = 13;
data_bits = 4096;
frozen_bits = frozen_8192_4128;
break;
case 25:
case 3:
mod_bits = 2;
cons_rows = 32;
code_order = 14;
data_bits = 8192;
frozen_bits = frozen_16384_8224;
break;
case 26:
case 4:
mod_bits = 4;
cons_rows = 4;
code_order = 12;
data_bits = 2048;
frozen_bits = frozen_4096_2080;
break;
case 27:
case 5:
mod_bits = 4;
cons_rows = 8;
code_order = 13;
data_bits = 4096;
frozen_bits = frozen_8192_4128;
break;
case 28:
case 6:
mod_bits = 4;
cons_rows = 16;
code_order = 14;
data_bits = 8192;
frozen_bits = frozen_16384_8224;
break;
case 29:
case 7:
mod_bits = 6;
cons_rows = 6;
code_order = 13;
data_bits = 4096;
frozen_bits = frozen_8192_4128;
break;
case 30:
case 8:
mod_bits = 6;
cons_rows = 11;
code_order = 14;
@ -337,15 +293,14 @@ struct Encoder
int data_bytes = data_bits / 8;
int offset = (freq_off * symbol_len) / rate;
mls0_off = offset - mls0_len / 2;
mls1_off = offset - mls1_len / 2;
cons_cols = code_cols + comb_cols;
code_off = offset - cons_cols / 2;
if (oper_mode > 0) {
value kern_fac = 1 / value(10 * reserved_tones);
value mag = 1 / value(10 * reserved_tones);
for (int i = 0, j = code_off - reserved_tones / 2; i < reserved_tones; ++i, ++j) {
if (j == code_off)
j += cons_cols;
fdom[bin(j)] = kern_fac;
fdom[bin(j)] = mag;
}
bwd(kern, fdom);
}
@ -359,9 +314,18 @@ struct Encoder
weight[i] = 1;
papr_min = 1000, papr_max = -1000;
pilot_block();
hadamardenc(mode, oper_mode);
if (!oper_mode) {
schmidl_cox();
meta_data(call_sign << 8);
CODE::MLS seq1(mls1_poly);
for (int i = 0, m = 0; i < cons_cols; ++i) {
if (i % comb_dist == comb_off) {
fdom[bin(i+code_off)] = nrz(seq1()) * mode[m++];
} else {
fdom[bin(i+code_off)] = 0;
}
}
symbol(false);
}
for (int input_index = 0; input_index < input_count; ++input_index) {
const char *input_name = input_names[input_index];
@ -378,7 +342,6 @@ struct Encoder
for (int i = 0; i < data_bytes; ++i)
input_data[i] ^= scrambler();
schmidl_cox();
meta_data((call_sign << 8) | oper_mode);
for (int i = 0; i < data_bits; ++i)
mesg[i] = nrz(CODE::get_le_bit(input_data, i));
crc1.reset();
@ -388,16 +351,13 @@ struct Encoder
mesg[i+data_bits] = nrz((crc1()>>i)&1);
polarenc(code, mesg, frozen_bits, code_order);
shuffle(perm, code);
for (int i = 0; i < cons_cols; ++i)
prev[i] = fdom[bin(i+code_off)];
CODE::MLS seq0(mls0_poly);
CODE::MLS seq1(mls1_poly);
for (int j = 0, k = 0; j < cons_rows; ++j) {
for (int i = 0; i < cons_cols; ++i) {
for (int i = 0, m = 0; i < cons_cols; ++i) {
if (i % comb_dist == comb_off) {
prev[i] *= nrz(seq0());
fdom[bin(i+code_off)] = prev[i];
fdom[bin(i+code_off)] = nrz(seq1()) * mode[m++];
} else {
fdom[bin(i+code_off)] = prev[i] * mod_map(perm+k);
fdom[bin(i+code_off)] = mod_map(perm+k);
k += mod_bits;
}
}
@ -411,27 +371,10 @@ struct Encoder
}
};
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 < 8) {
std::cerr << "usage: " << argv[0] << " OUTPUT RATE BITS CHANNELS OFFSET MODE CALLSIGN INPUT.." << std::endl;
if (argc < 7) {
std::cerr << "usage: " << argv[0] << " OUTPUT RATE BITS CHANNELS OFFSET MODE INPUT.." << std::endl;
return 1;
}
@ -447,19 +390,13 @@ int main(int argc, char **argv)
std::cerr << "Frequency offset must be divisible by 50." << std::endl;
return 1;
}
int input_count = argc - 8;
int input_count = argc - 7;
int oper_mode = std::atoi(argv[6]);
if (!oper_mode != !input_count) {
std::cerr << "Using operation mode " << oper_mode << " but " << input_count << " input file" << (input_count == 1 ? "" : "s") << " provided." << 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;
}
if (oper_mode && (oper_mode < 23 || oper_mode > 30)) {
if (oper_mode < 0 || oper_mode > 8) {
std::cerr << "Unsupported operation mode." << std::endl;
return 1;
}
@ -474,16 +411,16 @@ int main(int argc, char **argv)
output_file.silence(output_rate);
switch (output_rate) {
case 8000:
delete new Encoder<value, cmplx, 8000>(&output_file, argv+8, input_count, freq_off, call_sign, oper_mode);
delete new Encoder<value, cmplx, 8000>(&output_file, argv+7, input_count, freq_off, oper_mode);
break;
case 16000:
delete new Encoder<value, cmplx, 16000>(&output_file, argv+8, input_count, freq_off, call_sign, oper_mode);
delete new Encoder<value, cmplx, 16000>(&output_file, argv+7, input_count, freq_off, oper_mode);
break;
case 44100:
delete new Encoder<value, cmplx, 44100>(&output_file, argv+8, input_count, freq_off, call_sign, oper_mode);
delete new Encoder<value, cmplx, 44100>(&output_file, argv+7, input_count, freq_off, oper_mode);
break;
case 48000:
delete new Encoder<value, cmplx, 48000>(&output_file, argv+8, input_count, freq_off, call_sign, oper_mode);
delete new Encoder<value, cmplx, 48000>(&output_file, argv+7, input_count, freq_off, oper_mode);
break;
default:
std::cerr << "Unsupported sample rate." << std::endl;