Further minor updates and cleanup

pull/12/merge
Karlis Goba 2021-08-13 12:39:37 +03:00
rodzic d46686bccc
commit f7ca67be9f
4 zmienionych plików z 147 dodań i 145 usunięć

Wyświetl plik

@ -17,8 +17,8 @@
#define LOG_LEVEL LOG_INFO
const int kMin_score = 20; // Minimum sync score threshold for candidates
const int kMax_candidates = 160;
const int kMin_score = 10; // Minimum sync score threshold for candidates
const int kMax_candidates = 150;
const int kLDPC_iterations = 25;
const int kMax_decoded_messages = 50;
@ -30,7 +30,7 @@ const float kFSK_dev = 6.25f; // tone deviation in Hz and symbol rate
void usage()
{
fprintf(stderr, "Decode a 15-second WAV file.\n");
fprintf(stderr, "Decode a 15-second (or slighly shorter) WAV file.\n");
}
float hann_i(int i, int N)
@ -68,16 +68,16 @@ static float max2(float a, float b)
}
// Compute FFT magnitudes (log power) for each timeslot in the signal
void extract_power(const float signal[], waterfall_t *power)
void extract_power(const float signal[], waterfall_t *power, int block_size)
{
const int block_size = 2 * power->num_bins; // Average over 2 bins per FSK tone
const int subblock_size = block_size / power->time_osr;
const int nfft = block_size * power->freq_osr; // We take FFT of two blocks, advancing by one
const int nfft = block_size * power->freq_osr;
const float fft_norm = 2.0f / nfft;
float window[nfft];
for (int i = 0; i < nfft; ++i)
{
// window[i] = 1;
window[i] = hann_i(i, nfft);
// window[i] = (i < block_size) ? hamming_i(i, block_size) : 0;
// window[i] = blackman_i(i, nfft);
@ -200,8 +200,8 @@ int main(int argc, char **argv)
normalize_signal(signal, num_samples);
// Compute DSP parameters that depend on the sample rate
const int num_bins = (int)(sample_rate / (2 * kFSK_dev));
const int block_size = 2 * num_bins;
const int num_bins = (int)(sample_rate / (2 * kFSK_dev)); // number bins of FSK tone width that the spectrum can be divided into
const int block_size = (int)(sample_rate / kFSK_dev); // samples corresponding to one FSK symbol
const int subblock_size = block_size / kTime_osr;
const int nfft = block_size * kFreq_osr;
const int num_blocks = (num_samples - nfft + subblock_size) / block_size;
@ -216,7 +216,7 @@ int main(int argc, char **argv)
.time_osr = kTime_osr,
.freq_osr = kFreq_osr,
.mag = mag_power};
extract_power(signal, &power);
extract_power(signal, &power, block_size);
// Find top candidates by Costas sync score and localize them in time and frequency
candidate_t candidate_list[kMax_candidates];

Wyświetl plik

@ -29,7 +29,7 @@ static int get_index(const waterfall_t *power, int block, int time_sub, int freq
int find_sync(const waterfall_t *power, int num_candidates, candidate_t heap[], int min_score)
{
int heap_size = 0;
int num_alt = power->time_osr * power->freq_osr;
int sym_stride = power->time_osr * power->freq_osr * power->num_bins;
// Here we allow time offsets that exceed signal boundaries, as long as we still have all data bits.
// I.e. we can afford to skip the first 7 or the last 7 Costas symbols, as long as we track how many
@ -62,11 +62,10 @@ int find_sync(const waterfall_t *power, int num_candidates, candidate_t heap[],
// Weighted difference between the expected and all other symbols
// Does not work as well as the alternative score below
score += 8 * p8[kFT8_Costas_pattern[k]] -
p8[0] - p8[1] - p8[2] - p8[3] -
p8[4] - p8[5] - p8[6] - p8[7];
++num_symbols;
// score += 8 * p8[kFT8_Costas_pattern[k]] -
// p8[0] - p8[1] - p8[2] - p8[3] -
// p8[4] - p8[5] - p8[6] - p8[7];
// ++num_symbols;
// Check only the neighbors of the expected symbol frequency- and time-wise
int sm = kFT8_Costas_pattern[k]; // Index of the expected bin
@ -85,13 +84,13 @@ int find_sync(const waterfall_t *power, int num_candidates, candidate_t heap[],
if (k > 0)
{
// look one symbol back in time
score += p8[sm] - p8[sm - num_alt * power->num_bins];
score += p8[sm] - p8[sm - sym_stride];
++num_symbols;
}
if (k < 6)
{
// look one symbol forward in time
score += p8[sm] - p8[sm + num_alt * power->num_bins];
score += p8[sm] - p8[sm + sym_stride];
++num_symbols;
}
}
@ -133,7 +132,7 @@ int find_sync(const waterfall_t *power, int num_candidates, candidate_t heap[],
void extract_likelihood(const waterfall_t *power, const candidate_t *cand, float *log174)
{
int num_alt = power->time_osr * power->freq_osr;
int sym_stride = power->time_osr * power->freq_osr * power->num_bins;
int offset = get_index(power, cand->time_offset, cand->time_sub, cand->freq_sub, cand->freq_offset);
// Go over FSK tones and skip Costas sync symbols
@ -147,7 +146,7 @@ void extract_likelihood(const waterfall_t *power, const candidate_t *cand, float
int bit_idx = 3 * k;
// Pointer to 8 bins of the current symbol
const uint8_t *ps = power->mag + (offset + sym_idx * num_alt * power->num_bins);
const uint8_t *ps = power->mag + offset + (sym_idx * sym_stride);
decode_symbol(ps, bit_idx, log174);
}
@ -178,7 +177,7 @@ bool decode(const waterfall_t *power, const candidate_t *cand, message_t *messag
uint8_t plain174[FT8_N]; // message bits (0/1)
bp_decode(log174, max_iterations, plain174, &status->ldpc_errors);
// ldpc_decode(log174, kLDPC_iterations, plain174, &n_errors);
// ldpc_decode(log174, max_iterations, plain174, &status->ldpc_errors);
if (status->ldpc_errors > 0)
{
@ -292,10 +291,6 @@ static void decode_symbol(const uint8_t *power, int bit_idx, float *log174)
// Compute unnormalized log likelihood log(p(1) / p(0)) of bits corresponding to several FSK symbols at once
static void decode_multi_symbols(const uint8_t *power, int num_bins, int n_syms, int bit_idx, float *log174)
{
// The following section implements what seems to be multiple-symbol decode at one go,
// corresponding to WSJT-X's ft8b.f90. Experimentally found not to be any better than
// 1-symbol decode.
const int n_bits = 3 * n_syms;
const int n_tones = (1 << n_bits);
@ -321,10 +316,6 @@ static void decode_multi_symbols(const uint8_t *power, int num_bins, int n_syms,
s2[j] += (float)power[kFT8_Gray_map[j2] + 4 * num_bins];
s2[j] += (float)power[kFT8_Gray_map[j1] + 8 * num_bins];
}
// No need to go back to linear scale any more. Works better in dB.
// for (int j = 0; j < n_tones; ++j) {
// s2[j] = powf(10.0f, 0.1f * s2[j]);
// }
// Extract bit significance (and convert them to float)
// 8 FSK tones = 3 bits

Wyświetl plik

@ -83,7 +83,7 @@ void ldpc_decode(float codeword[], int max_iters, uint8_t plain[], int *ok)
a *= fast_tanh(-m[j][i2] / 2.0f);
}
}
e[j][i1] = logf((1 - a) / (1 + a));
e[j][i1] = -2.0f * fast_atanh(a);
}
}
@ -139,12 +139,12 @@ static int ldpc_check(uint8_t codeword[])
{
int errors = 0;
for (int j = 0; j < FT8_M; ++j)
for (int m = 0; m < FT8_M; ++m)
{
uint8_t x = 0;
for (int i = 0; i < kFT8_LDPC_num_rows[j]; ++i)
for (int i = 0; i < kFT8_LDPC_num_rows[m]; ++i)
{
x ^= codeword[kFT8_LDPC_Nm[j][i] - 1];
x ^= codeword[kFT8_LDPC_Nm[m][i] - 1];
}
if (x != 0)
{
@ -161,29 +161,20 @@ void bp_decode(float codeword[], int max_iters, uint8_t plain[], int *ok)
int min_errors = FT8_M;
int nclast = 0;
int ncnt = 0;
// initialize message data
for (int i = 0; i < FT8_N; ++i)
for (int n = 0; n < FT8_N; ++n)
{
for (int j = 0; j < 3; ++j)
{
tov[i][j] = 0;
}
tov[n][0] = tov[n][1] = tov[n][2] = 0;
}
for (int iter = 0; iter < max_iters; ++iter)
{
float zn[FT8_N];
// Do a hard decision guess (tov=0 in iter 0)
int plain_sum = 0;
// Update bit log likelihood ratios (tov=0 in iter 0)
for (int i = 0; i < FT8_N; ++i)
for (int n = 0; n < FT8_N; ++n)
{
zn[i] = codeword[i] + tov[i][0] + tov[i][1] + tov[i][2];
plain[i] = (zn[i] > 0) ? 1 : 0;
plain_sum += plain[i];
plain[n] = ((codeword[n] + tov[n][0] + tov[n][1] + tov[n][2]) > 0) ? 1 : 0;
plain_sum += plain[n];
}
if (plain_sum == 0)
@ -207,40 +198,40 @@ void bp_decode(float codeword[], int max_iters, uint8_t plain[], int *ok)
}
// Send messages from bits to check nodes
for (int i = 0; i < FT8_M; ++i)
for (int m = 0; m < FT8_M; ++m)
{
for (int j = 0; j < kFT8_LDPC_num_rows[i]; ++j)
for (int n_idx = 0; n_idx < kFT8_LDPC_num_rows[m]; ++n_idx)
{
int ibj = kFT8_LDPC_Nm[i][j] - 1;
float Tnm = zn[ibj];
for (int kk = 0; kk < 3; ++kk)
int n = kFT8_LDPC_Nm[m][n_idx] - 1;
// for each (n, m)
float Tnm = codeword[n];
for (int m_idx = 0; m_idx < 3; ++m_idx)
{
// subtract off what the bit had received from the check
if ((kFT8_LDPC_Mn[ibj][kk] - 1) == i)
if ((kFT8_LDPC_Mn[n][m_idx] - 1) != m)
{
Tnm -= tov[ibj][kk];
break;
Tnm += tov[n][m_idx];
}
}
toc[i][j] = fast_tanh(-Tnm / 2);
toc[m][n_idx] = fast_tanh(-Tnm / 2); // == (exp(-Tnm)-1) / (exp(-Tnm)+1)
}
}
// send messages from check nodes to variable nodes
for (int i = 0; i < FT8_N; ++i)
for (int n = 0; n < FT8_N; ++n)
{
for (int j = 0; j < 3; ++j)
for (int m_idx = 0; m_idx < 3; ++m_idx)
{
int ichk = kFT8_LDPC_Mn[i][j] - 1;
int m = kFT8_LDPC_Mn[n][m_idx] - 1;
// for each (n, m)
float Tmn = 1.0f;
for (int k = 0; k < kFT8_LDPC_num_rows[ichk]; ++k)
for (int n_idx = 0; n_idx < kFT8_LDPC_num_rows[m]; ++n_idx)
{
if ((kFT8_LDPC_Nm[ichk][k] - 1) != i)
if ((kFT8_LDPC_Nm[m][n_idx] - 1) != n)
{
Tmn *= toc[ichk][k];
Tmn *= toc[m][n_idx]; // tanh(q(n', m) / 2)
}
}
tov[i][j] = 2 * fast_atanh(-Tmn);
tov[n][m_idx] = 2 * fast_atanh(-Tmn); // == log( (1-Tmn) / (1+Tmn) )
}
}
}
@ -248,12 +239,12 @@ void bp_decode(float codeword[], int max_iters, uint8_t plain[], int *ok)
*ok = min_errors;
}
// https://varietyofsound.wordpress.com/2011/02/14/efficient-tanh-computation-using-lamberts-continued-fraction/
// http://functions.wolfram.com/ElementaryFunctions/ArcTanh/10/0001/
// https://mathr.co.uk/blog/2017-09-06_approximating_hyperbolic_tangent.html
// Ideas for approximating tanh/atanh:
// * https://varietyofsound.wordpress.com/2011/02/14/efficient-tanh-computation-using-lamberts-continued-fraction/
// * http://functions.wolfram.com/ElementaryFunctions/ArcTanh/10/0001/
// * https://mathr.co.uk/blog/2017-09-06_approximating_hyperbolic_tangent.html
// * https://math.stackexchange.com/a/446411
// thank you Douglas Bagnall
// https://math.stackexchange.com/a/446411
static float fast_tanh(float x)
{
if (x < -4.97f)

Wyświetl plik

@ -9,7 +9,7 @@
// n28 is a 28-bit integer, e.g. n28a or n28b, containing all the
// call sign bits from a packed message.
int unpack28(uint32_t n28, uint8_t ip, uint8_t i3, char *result)
int unpack_callsign(uint32_t n28, uint8_t ip, uint8_t i3, char *result)
{
// Check for special tokens DE, QRZ, CQ, CQ_nnn, CQ_aaaa
if (n28 < NTOKENS)
@ -105,7 +105,7 @@ int unpack28(uint32_t n28, uint8_t ip, uint8_t i3, char *result)
return 0; // Success
}
int unpack_type1(const uint8_t *a77, uint8_t i3, char *field1, char *field2, char *field3)
int unpack_type1(const uint8_t *a77, uint8_t i3, char *call_to, char *call_de, char *extra)
{
uint32_t n28a, n28b;
uint16_t igrid4;
@ -127,35 +127,36 @@ int unpack_type1(const uint8_t *a77, uint8_t i3, char *field1, char *field2, cha
igrid4 |= (a77[9] >> 6);
// Unpack both callsigns
if (unpack28(n28a >> 1, n28a & 0x01, i3, field1) < 0)
if (unpack_callsign(n28a >> 1, n28a & 0x01, i3, call_to) < 0)
{
return -1;
}
if (unpack28(n28b >> 1, n28b & 0x01, i3, field2) < 0)
if (unpack_callsign(n28b >> 1, n28b & 0x01, i3, call_de) < 0)
{
return -2;
}
// Fix "CQ_" to "CQ " -> already done in unpack28()
// Fix "CQ_" to "CQ " -> already done in unpack_callsign()
// TODO: add to recent calls
// if (field1[0] != '<' && strlen(field1) >= 4) {
// save_hash_call(field1)
// if (call_to[0] != '<' && strlen(call_to) >= 4) {
// save_hash_call(call_to)
// }
// if (field2[0] != '<' && strlen(field2) >= 4) {
// save_hash_call(field2)
// if (call_de[0] != '<' && strlen(call_de) >= 4) {
// save_hash_call(call_de)
// }
char *dst = extra;
if (igrid4 <= MAXGRID4)
{
// Extract 4 symbol grid locator
char *dst = field3;
uint16_t n = igrid4;
if (ir > 0)
{
// In case of ir=1 add an "R" before grid
dst = stpcpy(dst, "R ");
}
uint16_t n = igrid4;
dst[4] = '\0';
dst[3] = '0' + (n % 10);
n /= 10;
@ -164,33 +165,38 @@ int unpack_type1(const uint8_t *a77, uint8_t i3, char *field1, char *field2, cha
dst[1] = 'A' + (n % 18);
n /= 18;
dst[0] = 'A' + (n % 18);
// if (ir > 0 && strncmp(field1, "CQ", 2) == 0) return -1;
// if (ir > 0 && strncmp(call_to, "CQ", 2) == 0) return -1;
}
else
{
// Extract report
int irpt = igrid4 - MAXGRID4;
// Check special cases first
if (irpt == 1)
field3[0] = '\0';
else if (irpt == 2)
strcpy(field3, "RRR");
else if (irpt == 3)
strcpy(field3, "RR73");
else if (irpt == 4)
strcpy(field3, "73");
else if (irpt >= 5)
// Check special cases first (irpt > 0 always)
switch (irpt)
{
char *dst = field3;
case 1:
extra[0] = '\0';
break;
case 2:
strcpy(dst, "RRR");
break;
case 3:
strcpy(dst, "RR73");
break;
case 4:
strcpy(dst, "73");
break;
default:
// Extract signal report as a two digit number with a + or - sign
if (ir > 0)
{
*dst++ = 'R'; // Add "R" before report
}
int_to_dd(dst, irpt - 35, 2, true);
break;
}
// if (irpt >= 2 && strncmp(field1, "CQ", 2) == 0) return -1;
// if (irpt >= 2 && strncmp(call_to, "CQ", 2) == 0) return -1;
}
return 0; // Success
@ -201,6 +207,7 @@ int unpack_text(const uint8_t *a71, char *text)
// TODO: test
uint8_t b71[9];
// Shift 71 bits right by 1 bit, so that it's right-aligned in the byte array
uint8_t carry = 0;
for (int i = 0; i < 9; ++i)
{
@ -231,7 +238,7 @@ int unpack_telemetry(const uint8_t *a71, char *telemetry)
{
uint8_t b71[9];
// Shift bits in a71 right by 1
// Shift bits in a71 right by 1 bit
uint8_t carry = 0;
for (int i = 0; i < 9; ++i)
{
@ -256,7 +263,7 @@ int unpack_telemetry(const uint8_t *a71, char *telemetry)
//none standard for wsjt-x 2.0
//by KD8CEC
int unpack_nonstandard(const uint8_t *a77, char *field1, char *field2, char *field3)
int unpack_nonstandard(const uint8_t *a77, char *call_to, char *call_de, char *extra)
{
uint32_t n12, iflip, nrpt, icq;
uint64_t n58;
@ -290,11 +297,11 @@ int unpack_nonstandard(const uint8_t *a77, char *field1, char *field2, char *fie
char call_3[15];
// should replace with hash12(n12, call_3);
// strcpy(call_3, "<...>");
call_3[0] = '<';
int_to_dd(call_3 + 1, n12, 4, false);
call_3[5] = '>';
call_3[6] = '\0';
strcpy(call_3, "<...>");
// call_3[0] = '<';
// int_to_dd(call_3 + 1, n12, 4, false);
// call_3[5] = '>';
// call_3[6] = '\0';
char *call_1 = (iflip) ? c11 : call_3;
char *call_2 = (iflip) ? call_3 : c11;
@ -302,62 +309,65 @@ int unpack_nonstandard(const uint8_t *a77, char *field1, char *field2, char *fie
if (icq == 0)
{
strcpy(field1, trim(call_1));
strcpy(call_to, trim(call_1));
if (nrpt == 1)
strcpy(field3, "RRR");
strcpy(extra, "RRR");
else if (nrpt == 2)
strcpy(field3, "RR73");
strcpy(extra, "RR73");
else if (nrpt == 3)
strcpy(field3, "73");
strcpy(extra, "73");
else
{
field3[0] = '\0';
extra[0] = '\0';
}
}
else
{
strcpy(field1, "CQ");
field3[0] = '\0';
strcpy(call_to, "CQ");
extra[0] = '\0';
}
strcpy(field2, trim(call_2));
strcpy(call_de, trim(call_2));
return 0;
}
int unpack77_fields(const uint8_t *a77, char *field1, char *field2, char *field3)
int unpack77_fields(const uint8_t *a77, char *call_to, char *call_de, char *extra)
{
uint8_t n3, i3;
call_to[0] = call_de[0] = extra[0] = '\0';
// Extract n3 (bits 71..73) and i3 (bits 74..76)
n3 = ((a77[8] << 2) & 0x04) | ((a77[9] >> 6) & 0x03);
i3 = (a77[9] >> 3) & 0x07;
// Extract i3 (bits 74..76)
uint8_t i3 = (a77[9] >> 3) & 0x07;
field1[0] = field2[0] = field3[0] = '\0';
if (i3 == 0 && n3 == 0)
if (i3 == 0)
{
// 0.0 Free text
return unpack_text(a77, field1);
}
// else if (i3 == 0 && n3 == 1) {
// // 0.1 K1ABC RR73; W9XYZ <KH1/KH7Z> -11 28 28 10 5 71 DXpedition Mode
// }
// else if (i3 == 0 && n3 == 2) {
// // 0.2 PA3XYZ/P R 590003 IO91NP 28 1 1 3 12 25 70 EU VHF contest
// }
// else if (i3 == 0 && (n3 == 3 || n3 == 4)) {
// // 0.3 WA9XYZ KA1ABC R 16A EMA 28 28 1 4 3 7 71 ARRL Field Day
// // 0.4 WA9XYZ KA1ABC R 32A EMA 28 28 1 4 3 7 71 ARRL Field Day
// }
else if (i3 == 0 && n3 == 5)
{
// 0.5 0123456789abcdef01 71 71 Telemetry (18 hex)
return unpack_telemetry(a77, field1);
// Extract n3 (bits 71..73)
uint8_t n3 = ((a77[8] << 2) & 0x04) | ((a77[9] >> 6) & 0x03);
if (n3 == 0)
{
// 0.0 Free text
return unpack_text(a77, extra);
}
// else if (i3 == 0 && n3 == 1) {
// // 0.1 K1ABC RR73; W9XYZ <KH1/KH7Z> -11 28 28 10 5 71 DXpedition Mode
// }
// else if (i3 == 0 && n3 == 2) {
// // 0.2 PA3XYZ/P R 590003 IO91NP 28 1 1 3 12 25 70 EU VHF contest
// }
// else if (i3 == 0 && (n3 == 3 || n3 == 4)) {
// // 0.3 WA9XYZ KA1ABC R 16A EMA 28 28 1 4 3 7 71 ARRL Field Day
// // 0.4 WA9XYZ KA1ABC R 32A EMA 28 28 1 4 3 7 71 ARRL Field Day
// }
else if (n3 == 5)
{
// 0.5 0123456789abcdef01 71 71 Telemetry (18 hex)
return unpack_telemetry(a77, extra);
}
}
else if (i3 == 1 || i3 == 2)
{
// Type 1 (standard message) or Type 2 ("/P" form for EU VHF contest)
return unpack_type1(a77, i3, field1, field2, field3);
return unpack_type1(a77, i3, call_to, call_de, extra);
}
// else if (i3 == 3) {
// // Type 3: ARRL RTTY Contest
@ -367,7 +377,7 @@ int unpack77_fields(const uint8_t *a77, char *field1, char *field2, char *field3
// // Type 4: Nonstandard calls, e.g. <WA9XYZ> PJ4/KA1ABC RR73
// // One hashed call or "CQ"; one compound or nonstandard call with up
// // to 11 characters; and (if not "CQ") an optional RRR, RR73, or 73.
return unpack_nonstandard(a77, field1, field2, field3);
return unpack_nonstandard(a77, call_to, call_de, extra);
}
// else if (i3 == 5) {
// // Type 5: TU; W9XYZ K1ABC R-09 FN 1 28 28 1 7 9 74 WWROF contest
@ -379,22 +389,32 @@ int unpack77_fields(const uint8_t *a77, char *field1, char *field2, char *field3
int unpack77(const uint8_t *a77, char *message)
{
char field1[14];
char field2[14];
char field3[7];
char call_to[14];
char call_de[14];
char extra[7];
int rc = unpack77_fields(a77, field1, field2, field3);
int rc = unpack77_fields(a77, call_to, call_de, extra);
if (rc < 0)
return rc;
// int msg_sz = strlen(call_to) + strlen(call_de) + strlen(extra) + 2;
char *dst = message;
// int msg_sz = strlen(field1) + strlen(field2) + strlen(field3) + 2;
dst = stpcpy(dst, field1);
*dst++ = ' ';
dst = stpcpy(dst, field2);
*dst++ = ' ';
dst = stpcpy(dst, field3);
dst[0] = '\0';
if (call_to[0] != '\0')
{
dst = stpcpy(dst, call_to);
*dst++ = ' ';
}
if (call_de[0] != '\0')
{
dst = stpcpy(dst, call_de);
*dst++ = ' ';
}
dst = stpcpy(dst, extra);
*dst = '\0';
return 0;