rpitx/ft8_lib
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity

Improved decode accuracy (better sync score), configurable time/freq OSR

Browse source
This commit is contained in:
Karlis Goba 2019-11-10 09:04:00 +02:00
parent c727f2b572
commit 43266baf76
4 changed files with 137 additions and 73 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

75
decode_ft8.cpp
View file

@ -19,7 +19,10 @@ const int kMax_candidates = 100;
const int kLDPC_iterations = 20;
const int kMax_decoded_messages = 50;
const int kMax_message_length = 20;
const int kMax_message_length = 25;
const int kFreq_osr = 2;
const int kTime_osr = 2;
void usage() {
@ -55,21 +58,27 @@ float blackman_i(int i, int N) {
return a0 - a1*x1 + a2*x2;
}
static float max2(float a, float b) {
return (a >= b) ? a : b;
}
// Compute FFT magnitudes (log power) for each timeslot in the signal
void extract_power(const float signal[], int num_blocks, int num_bins, uint8_t power[]) {
const int block_size = 2 * num_bins; // Average over 2 bins per FSK tone
const int nfft = 2 * block_size; // We take FFT of two blocks, advancing by one
void extract_power(const float signal[], ft8::MagArray * power) {
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 float fft_norm = 2.0f / nfft;
float window[nfft];
for (int i = 0; i < nfft; ++i) {
window[i] = blackman_i(i, nfft);
window[i] = hann_i(i, nfft);
}
size_t fft_work_size;
kiss_fftr_alloc(nfft, 0, 0, &fft_work_size);
LOG(LOG_INFO, "Block size = %d\n", block_size);
LOG(LOG_INFO, "Subblock size = %d\n", subblock_size);
LOG(LOG_INFO, "N_FFT = %d\n", nfft);
LOG(LOG_INFO, "FFT work area = %lu\n", fft_work_size);
@ -78,16 +87,16 @@ void extract_power(const float signal[], int num_blocks, int num_bins, uint8_t p
int offset = 0;
float max_mag = -100.0f;
for (int i = 0; i < num_blocks; ++i) {
for (int i = 0; i < power->num_blocks; ++i) {
// Loop over two possible time offsets (0 and block_size/2)
for (int time_sub = 0; time_sub <= block_size/2; time_sub += block_size/2) {
for (int time_sub = 0; time_sub < power->time_osr; ++time_sub) {
kiss_fft_scalar timedata[nfft];
kiss_fft_cpx freqdata[nfft/2 + 1];
float mag_db[nfft/2 + 1];
// Extract windowed signal block
for (int j = 0; j < nfft; ++j) {
timedata[j] = window[j] * signal[(i * block_size) + (j + time_sub)];
timedata[j] = window[j] * signal[(i * block_size) + (j + time_sub * subblock_size)];
}
kiss_fftr(fft_cfg, timedata, freqdata);
@ -99,15 +108,17 @@ void extract_power(const float signal[], int num_blocks, int num_bins, uint8_t p
}
// Loop over two possible frequency bin offsets (for averaging)
for (int freq_sub = 0; freq_sub < 2; ++freq_sub) {
for (int j = 0; j < num_bins; ++j) {
float db1 = mag_db[j * 2 + freq_sub];
float db2 = mag_db[j * 2 + freq_sub + 1];
float db = (db1 + db2) / 2;
for (int freq_sub = 0; freq_sub < power->freq_osr; ++freq_sub) {
for (int j = 0; j < power->num_bins; ++j) {
float db1 = mag_db[j * power->freq_osr + freq_sub];
//float db2 = mag_db[j * 2 + freq_sub + 1];
//float db = (db1 + db2) / 2;
float db = db1;
//float db = sqrtf(db1 * db2);
// Scale decibels to unsigned 8-bit range and clamp the value
int scaled = (int)(2 * (db + 120));
power[offset] = (scaled < 0) ? 0 : ((scaled > 255) ? 255 : scaled);
power->mag[offset] = (scaled < 0) ? 0 : ((scaled > 255) ? 255 : scaled);
++offset;
if (db > max_mag) max_mag = db;
@ -171,17 +182,26 @@ int main(int argc, char **argv) {
// Compute DSP parameters that depend on the sample rate
const int num_bins = (int)(sample_rate / (2 * fsk_dev));
const int block_size = 2 * num_bins;
const int num_blocks = (num_samples - (block_size/2) - block_size) / block_size;
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;
LOG(LOG_INFO, "%d blocks, %d bins\n", num_blocks, num_bins);
LOG(LOG_INFO, "Sample rate %d Hz, %d blocks, %d bins\n", sample_rate, num_blocks, num_bins);
// Compute FFT over the whole signal and store it
uint8_t power[num_blocks * 4 * num_bins];
extract_power(signal, num_blocks, num_bins, power);
uint8_t mag_power[num_blocks * kFreq_osr * kTime_osr * num_bins];
ft8::MagArray power = {
.num_blocks = num_blocks,
.num_bins = num_bins,
.time_osr = kTime_osr,
.freq_osr = kFreq_osr,
.mag = mag_power
};
extract_power(signal, &power);
// Find top candidates by Costas sync score and localize them in time and frequency
ft8::Candidate candidate_list[kMax_candidates];
int num_candidates = ft8::find_sync(power, num_blocks, num_bins, ft8::kCostas_map, kMax_candidates, candidate_list);
int num_candidates = ft8::find_sync(&power, ft8::kCostas_map, kMax_candidates, candidate_list);
// TODO: sort the candidates by strongest sync first?
@ -190,23 +210,32 @@ int main(int argc, char **argv) {
int num_decoded = 0;
for (int idx = 0; idx < num_candidates; ++idx) {
ft8::Candidate &cand = candidate_list[idx];
float freq_hz = (cand.freq_offset + cand.freq_sub / 2.0f) * fsk_dev;
float time_sec = (cand.time_offset + cand.time_sub / 2.0f) / fsk_dev;
float freq_hz = (cand.freq_offset + (float)cand.freq_sub / kFreq_osr) * fsk_dev;
float time_sec = (cand.time_offset + (float)cand.time_sub / kTime_osr) / fsk_dev;
float log174[ft8::N];
ft8::extract_likelihood(power, num_bins, cand, ft8::kGray_map, log174);
ft8::extract_likelihood(&power, cand, ft8::kGray_map, log174);
// bp_decode() produces better decodes, uses way less memory
uint8_t plain[ft8::N];
int n_errors = 0;
ft8::bp_decode(log174, kLDPC_iterations, plain, &n_errors);
//ldpc_decode(log174, kLDPC_iterations, plain, &n_errors);
//ft8::ldpc_decode(log174, kLDPC_iterations, plain, &n_errors);
if (n_errors > 0) {
LOG(LOG_DEBUG, "ldpc_decode() = %d (%.0f Hz)\n", n_errors, freq_hz);
continue;
}
int sum_plain = 0;
for (int i = 0; i < ft8::N; ++i) {
sum_plain += plain[i];
}
if (sum_plain == 0) {
// All zeroes message
continue;
}
// Extract payload + CRC (first ft8::K bits)
uint8_t a91[ft8::K_BYTES];
ft8::pack_bits(plain, ft8::K, a91);

117
ft8/decode.cpp
View file

@ -13,65 +13,87 @@ static void heapify_up(Candidate *heap, int heap_size);
static void decode_symbol(const uint8_t *power, const uint8_t *code_map, int bit_idx, float *log174);
static void decode_multi_symbols(const uint8_t *power, int num_bins, int n_syms, const uint8_t *code_map, int bit_idx, float *log174);
static int get_index(const MagArray *power, int block, int time_sub, int freq_sub, int bin) {
return ((((block * power->time_osr) + time_sub) * power->freq_osr + freq_sub) * power->num_bins) + bin;
}
// Localize top N candidates in frequency and time according to their sync strength (looking at Costas symbols)
// We treat and organize the candidate list as a min-heap (empty initially).
int find_sync(const uint8_t *power, int num_blocks, int num_bins, const uint8_t *sync_map, int num_candidates, Candidate *heap) {
int find_sync(const MagArray *power, const uint8_t *sync_map, int num_candidates, Candidate *heap) {
int heap_size = 0;
int num_alt = power->time_osr * power->freq_osr;
// 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
// sync symbols we included in the score, so the score is averaged.
for (int alt = 0; alt < 4; ++alt) {
for (int time_offset = -7; time_offset < num_blocks - ft8::NN + 7; ++time_offset) {
for (int freq_offset = 0; freq_offset < num_bins - 8; ++freq_offset) {
int score = 0;
for (int time_sub = 0; time_sub < power->time_osr; ++time_sub) {
for (int freq_sub = 0; freq_sub < power->freq_osr; ++freq_sub) {
for (int time_offset = -7; time_offset < power->num_blocks - ft8::NN + 7; ++time_offset) {
for (int freq_offset = 0; freq_offset < power->num_bins - 8; ++freq_offset) {
int score = 0;
// Compute average score over sync symbols (m+k = 0-7, 36-43, 72-79)
int num_symbols = 0;
for (int m = 0; m <= 72; m += 36) {
for (int k = 0; k < 7; ++k) {
// Check for time boundaries
if (time_offset + k + m < 0) continue;
if (time_offset + k + m >= num_blocks) break;
// Compute average score over sync symbols (m+k = 0-7, 36-43, 72-79)
int num_symbols = 0;
for (int m = 0; m <= 72; m += 36) {
for (int k = 0; k < 7; ++k) {
// Check for time boundaries
if (time_offset + k + m < 0) continue;
if (time_offset + k + m >= power->num_blocks) break;
int offset = ((time_offset + k + m) * 4 + alt) * num_bins + freq_offset;
const uint8_t *p8 = power + offset;
// int offset = ((time_offset + k + m) * num_alt + alt) * power->num_bins + freq_offset;
int offset = get_index(power, time_offset + k + m, time_sub, freq_sub, freq_offset);
const uint8_t *p8 = power->mag + offset;
score += 8 * p8[sync_map[k]] -
p8[0] - p8[1] - p8[2] - p8[3] -
p8[4] - p8[5] - p8[6] - p8[7];
// Weighted difference between the expected and all other symbols
// Does not work as well as the alternative score below
// score += 8 * p8[sync_map[k]] -
// p8[0] - p8[1] - p8[2] - p8[3] -
// p8[4] - p8[5] - p8[6] - p8[7];
// Check only the neighbors of the expected symbol frequency- and time-wise
int sm = sync_map[k]; // Index of the expected bin
if (sm > 0) {
// look at one frequency bin lower
score += p8[sm] - p8[sm - 1];
}
if (sm < 7) {
// look at one frequency bin higher
score += p8[sm] - p8[sm + 1];
}
if (k > 0) {
// look one symbol back in time
score += p8[sm] - p8[sm - num_alt * power->num_bins];
}
if (k < 6) {
// look one symbol forward in time
score += p8[sm] - p8[sm + num_alt * power->num_bins];
}
// int sm = sync_map[k];
// score += 4 * (int)p8[sm];
// if (sm > 0) score -= p8[sm - 1];
// if (sm < 7) score -= p8[sm + 1];
// if (k > 0) score -= p8[sm - 4 * num_bins];
// if (k < 6) score -= p8[sm + 4 * num_bins];
++num_symbols;
++num_symbols;
}
}
}
score /= num_symbols;
score /= num_symbols;
// If the heap is full AND the current candidate is better than
// the worst in the heap, we remove the worst and make space
if (heap_size == num_candidates && score > heap[0].score) {
heap[0] = heap[heap_size - 1];
--heap_size;
// If the heap is full AND the current candidate is better than
// the worst in the heap, we remove the worst and make space
if (heap_size == num_candidates && score > heap[0].score) {
heap[0] = heap[heap_size - 1];
--heap_size;
heapify_down(heap, heap_size);
}
heapify_down(heap, heap_size);
}
// If there's free space in the heap, we add the current candidate
if (heap_size < num_candidates) {
heap[heap_size].score = score;
heap[heap_size].time_offset = time_offset;
heap[heap_size].freq_offset = freq_offset;
heap[heap_size].time_sub = alt / 2;
heap[heap_size].freq_sub = alt % 2;
++heap_size;
// If there's free space in the heap, we add the current candidate
if (heap_size < num_candidates) {
heap[heap_size].score = score;
heap[heap_size].time_offset = time_offset;
heap[heap_size].freq_offset = freq_offset;
heap[heap_size].time_sub = time_sub;
heap[heap_size].freq_sub = freq_sub;
++heap_size;
heapify_up(heap, heap_size);
heapify_up(heap, heap_size);
}
}
}
}
@ -83,19 +105,22 @@ int find_sync(const uint8_t *power, int num_blocks, int num_bins, const uint8_t
// Compute log likelihood log(p(1) / p(0)) of 174 message bits
// for later use in soft-decision LDPC decoding
void extract_likelihood(const uint8_t *power, int num_bins, const Candidate & cand, const uint8_t *code_map, float *log174) {
int offset = (cand.time_offset * 4 + cand.time_sub * 2 + cand.freq_sub) * num_bins + cand.freq_offset;
void extract_likelihood(const MagArray *power, const Candidate & cand, const uint8_t *code_map, float *log174) {
int num_alt = power->time_osr * power->freq_osr;
// int offset = (cand.time_offset * num_alt + cand.time_sub * power->freq_osr + cand.freq_sub) * power->num_bins + cand.freq_offset;
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
const int n_syms = 1;
const int n_bits = 3 * n_syms;
const int n_tones = (1 << n_bits);
for (int k = 0; k < ft8::ND; k += n_syms) {
// Add either 7 or 14 extra symbols to account for sync
int sym_idx = (k < ft8::ND / 2) ? (k + 7) : (k + 14);
int bit_idx = 3 * k;
// Pointer to 8 bins of the current symbol
const uint8_t *ps = power + (offset + sym_idx * 4 * num_bins);
const uint8_t *ps = power->mag + (offset + sym_idx * num_alt * power->num_bins);
decode_symbol(ps, code_map, bit_idx, log174);
}

12
ft8/decode.h
View file

@ -4,6 +4,14 @@
namespace ft8 {
struct MagArray {
int num_blocks; // number of total blocks (symbols)
int num_bins; // number of FFT bins
int time_osr; // number of time subdivisions
int freq_osr; // number of frequency subdivisions
uint8_t * mag; // FFT magnitudes as [blocks][time_sub][freq_sub][num_bins]
};
struct Candidate {
int16_t score;
int16_t time_offset;
@ -15,11 +23,11 @@ struct Candidate {
// Localize top N candidates in frequency and time according to their sync strength (looking at Costas symbols)
// We treat and organize the candidate list as a min-heap (empty initially).
int find_sync(const uint8_t *power, int num_blocks, int num_bins, const uint8_t *sync_map, int num_candidates, Candidate *heap);
int find_sync(const MagArray * power, const uint8_t *sync_map, int num_candidates, Candidate *heap);
// Compute log likelihood log(p(1) / p(0)) of 174 message bits
// for later use in soft-decision LDPC decoding
void extract_likelihood(const uint8_t *power, int num_bins, const Candidate & cand, const uint8_t *code_map, float *log174);
void extract_likelihood(const MagArray *power, const Candidate & cand, const uint8_t *code_map, float *log174);
}

6
ft8/unpack.cpp
View file

@ -199,6 +199,8 @@ int unpack_text(const uint8_t *a71, char *text) {
carry = (a71[i] & 1) ? 0x80 : 0;
}
char c14[14];
c14[13] = 0;
for (int idx = 12; idx >= 0; --idx) {
// Divide the long integer in b71 by 42
uint16_t rem = 0;
@ -207,10 +209,10 @@ int unpack_text(const uint8_t *a71, char *text) {
b71[i] = rem / 42;
rem = rem % 42;
}
text[idx] = charn(rem, 0);
c14[idx] = charn(rem, 0);
}
text[13] = '\0';
strcpy(text, trim(c14));
return 0; // Success
}