kopia lustrzana https://github.com/F5OEO/ft8_lib
Moved decoding DSP routines to separate file
Karlis Goba
rodzic
0c2e35b998
commit
0c4efc2d87
2
Makefile
2
Makefile
|
|
@ -16,7 +16,7 @@ gen_ft8: gen_ft8.o ft8/constants.o ft8/text.o ft8/pack_v2.o ft8/encode_v2.o comm
|
|||
test: test.o ft8/v1/encode.o ft8/v1/pack.o ft8/v1/unpack.o ft8/pack_v2.o ft8/encode_v2.o ft8/text.o ft8/constants.o
|
||||
$(CXX) $(LDFLAGS) -o $@ $^
|
||||
|
||||
decode_ft8: decode_ft8.o fft/kiss_fftr.o fft/kiss_fft.o ft8/ldpc.o ft8/unpack_v2.o ft8/text.o ft8/constants.o common/wave.o
|
||||
decode_ft8: decode_ft8.o fft/kiss_fftr.o fft/kiss_fft.o ft8/decode.o ft8/ldpc.o ft8/unpack_v2.o ft8/text.o ft8/constants.o common/wave.o
|
||||
$(CXX) $(LDFLAGS) -o $@ $^
|
||||
|
||||
clean:
|
||||
|
|
|
|||
266
decode_ft8.cpp
266
decode_ft8.cpp
|
|
@ -5,6 +5,7 @@
|
|||
|
||||
#include "ft8/unpack_v2.h"
|
||||
#include "ft8/ldpc.h"
|
||||
#include "ft8/decode.h"
|
||||
#include "ft8/constants.h"
|
||||
|
||||
#include "common/wave.h"
|
||||
|
|
@ -43,121 +44,15 @@ float blackman_i(int i, int N) {
|
|||
const float a2 = alpha / 2;
|
||||
|
||||
float x1 = cosf(2 * (float)M_PI * i / (N - 1));
|
||||
float x2 = cosf(4 * (float)M_PI * i / (N - 1));
|
||||
//float x2 = cosf(4 * (float)M_PI * i / (N - 1));
|
||||
float x2 = 2*x1*x1 - 1; // Use double angle formula
|
||||
|
||||
return a0 - a1*x1 + a2*x2;
|
||||
}
|
||||
|
||||
|
||||
struct Candidate {
|
||||
int16_t score;
|
||||
int16_t time_offset;
|
||||
int16_t freq_offset;
|
||||
uint8_t time_sub;
|
||||
uint8_t freq_sub;
|
||||
};
|
||||
|
||||
|
||||
void heapify_down(Candidate *heap, int heap_size) {
|
||||
// heapify from the root down
|
||||
int current = 0;
|
||||
while (true) {
|
||||
int largest = current;
|
||||
int left = 2 * current + 1;
|
||||
int right = left + 1;
|
||||
|
||||
if (left < heap_size && heap[left].score < heap[largest].score) {
|
||||
largest = left;
|
||||
}
|
||||
if (right < heap_size && heap[right].score < heap[largest].score) {
|
||||
largest = right;
|
||||
}
|
||||
if (largest == current) {
|
||||
break;
|
||||
}
|
||||
|
||||
Candidate tmp = heap[largest];
|
||||
heap[largest] = heap[current];
|
||||
heap[current] = tmp;
|
||||
current = largest;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
void heapify_up(Candidate *heap, int heap_size) {
|
||||
// heapify from the last node up
|
||||
int current = heap_size - 1;
|
||||
while (current > 0) {
|
||||
int parent = (current - 1) / 2;
|
||||
if (heap[current].score >= heap[parent].score) {
|
||||
break;
|
||||
}
|
||||
|
||||
Candidate tmp = heap[parent];
|
||||
heap[parent] = heap[current];
|
||||
heap[current] = tmp;
|
||||
current = parent;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// Find 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 heap_size = 0;
|
||||
|
||||
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;
|
||||
|
||||
// Compute score over Costas symbols (0-7, 36-43, 72-79)
|
||||
int num_scores = 0;
|
||||
for (int m = 0; m <= 72; m += 36) {
|
||||
for (int k = 0; k < 7; ++k) {
|
||||
if (time_offset + k + m < 0) continue;
|
||||
if (time_offset + k + m >= num_blocks) break;
|
||||
int offset = ((time_offset + k + m) * 4 + alt) * num_bins + freq_offset;
|
||||
score += 8 * (int)power[offset + sync_map[k]] -
|
||||
power[offset + 0] - power[offset + 1] -
|
||||
power[offset + 2] - power[offset + 3] -
|
||||
power[offset + 4] - power[offset + 5] -
|
||||
power[offset + 6] - power[offset + 7];
|
||||
++num_scores;
|
||||
}
|
||||
}
|
||||
score /= num_scores;
|
||||
|
||||
// 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);
|
||||
}
|
||||
|
||||
// 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;
|
||||
|
||||
heapify_up(heap, heap_size);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
return heap_size;
|
||||
}
|
||||
|
||||
|
||||
// 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) {
|
||||
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
|
||||
|
||||
|
|
@ -172,11 +67,12 @@ void extract_power(const float *signal, int num_blocks, int num_bins, uint8_t *p
|
|||
printf("N_FFT = %d\n", nfft);
|
||||
printf("FFT work area = %lu\n", fft_work_size);
|
||||
|
||||
void * fft_work = malloc(fft_work_size);
|
||||
void *fft_work = malloc(fft_work_size);
|
||||
kiss_fftr_cfg fft_cfg = kiss_fftr_alloc(nfft, 0, fft_work, &fft_work_size);
|
||||
|
||||
// Currently bit unsure about the scaling factor of kiss FFT
|
||||
int offset = 0;
|
||||
float fft_norm = 1.0f / nfft;
|
||||
float fft_norm = 1.0f / nfft / sqrtf(nfft);
|
||||
float max_mag = -100.0f;
|
||||
for (int i = 0; i < num_blocks; ++i) {
|
||||
// Loop over two possible time offsets (0 and block_size/2)
|
||||
|
|
@ -194,8 +90,8 @@ void extract_power(const float *signal, int num_blocks, int num_bins, uint8_t *p
|
|||
|
||||
// Compute log magnitude in decibels
|
||||
for (int j = 0; j < nfft/2 + 1; ++j) {
|
||||
float mag2 = (freqdata[j].i * freqdata[j].i + freqdata[j].r * freqdata[j].r);
|
||||
mag_db[j] = 10.0f * log10f(1.0E-10f + mag2 * fft_norm);
|
||||
float mag2 = fft_norm * (freqdata[j].i * freqdata[j].i + freqdata[j].r * freqdata[j].r);
|
||||
mag_db[j] = 10.0f * log10f(1E-10f + mag2);
|
||||
}
|
||||
|
||||
// Loop over two possible frequency bin offsets (for averaging)
|
||||
|
|
@ -221,100 +117,6 @@ void extract_power(const float *signal, int num_blocks, int num_bins, uint8_t *p
|
|||
}
|
||||
|
||||
|
||||
float max2(float a, float b) {
|
||||
return (a >= b) ? a : b;
|
||||
}
|
||||
|
||||
|
||||
float max4(float a, float b, float c, float d) {
|
||||
return max2(max2(a, b), max2(c, d));
|
||||
}
|
||||
|
||||
|
||||
// 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;
|
||||
|
||||
// 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) {
|
||||
int sym_idx = (k < FT8_ND / 2) ? (k + 7) : (k + 14);
|
||||
|
||||
// Pointer to 8 bins of the current symbol
|
||||
const uint8_t *ps = power + (offset + sym_idx * 4 * num_bins);
|
||||
float s2[n_tones];
|
||||
|
||||
for (int j = 0; j < n_tones; ++j) {
|
||||
int j1 = j & 0x07;
|
||||
s2[j] = (float)ps[code_map[j1]];
|
||||
//int j2 = (j >> 3) & 0x07;
|
||||
//s2[j] = (float)ps[code_map[j2]];
|
||||
//s2[j] += (float)ps[code_map[j1] + 4 * num_bins];
|
||||
}
|
||||
|
||||
// Extract bit significance (and convert them to float)
|
||||
// 8 FSK tones = 3 bits
|
||||
int bit_idx = 3 * k;
|
||||
for (int i = 0; i < n_bits; ++i) {
|
||||
uint16_t mask = (n_tones >> (i + 1));
|
||||
|
||||
float max_zero = -1000, max_one = -1000;
|
||||
for (int n = 0; n < n_tones; ++n) {
|
||||
if (n & mask) {
|
||||
max_one = max2(max_one, s2[n]);
|
||||
}
|
||||
else {
|
||||
max_zero = max2(max_zero, s2[n]);
|
||||
}
|
||||
}
|
||||
if (bit_idx + i >= 174) break;
|
||||
log174[bit_idx + i] = max_one - max_zero;
|
||||
}
|
||||
// log174[bit_idx + 0] = max4(s2[4], s2[5], s2[6], s2[7]) - max4(s2[0], s2[1], s2[2], s2[3]);
|
||||
// log174[bit_idx + 1] = max4(s2[2], s2[3], s2[6], s2[7]) - max4(s2[0], s2[1], s2[4], s2[5]);
|
||||
// log174[bit_idx + 2] = max4(s2[1], s2[3], s2[5], s2[7]) - max4(s2[0], s2[2], s2[4], s2[6]);
|
||||
}
|
||||
|
||||
// Compute the variance of log174
|
||||
float sum = 0;
|
||||
float sum2 = 0;
|
||||
float inv_n = 1.0f / FT8_N;
|
||||
for (int i = 0; i < FT8_N; ++i) {
|
||||
sum += log174[i];
|
||||
sum2 += log174[i] * log174[i];
|
||||
}
|
||||
float variance = (sum2 - sum * sum * inv_n) * inv_n;
|
||||
|
||||
// Normalize log174 such that sigma = 2.83 (Why? It's in WSJT-X)
|
||||
float norm_factor = 3.83f / sqrtf(variance);
|
||||
for (int i = 0; i < FT8_N; ++i) {
|
||||
log174[i] *= norm_factor;
|
||||
//printf("%.1f ", log174[i]);
|
||||
}
|
||||
//printf("\n");
|
||||
}
|
||||
|
||||
|
||||
void test_tones(float *log174) {
|
||||
for (int i = 0; i < FT8_ND; ++i) {
|
||||
const uint8_t inv_map[8] = {0, 1, 3, 2, 6, 4, 5, 7};
|
||||
uint8_t tone = ("0000000011721762454112705354533170166234757420515470163426"[i]) - '0';
|
||||
uint8_t b3 = inv_map[tone];
|
||||
log174[3 * i] = (b3 & 4) ? +1.0 : -1.0;
|
||||
log174[3 * i + 1] = (b3 & 2) ? +1.0 : -1.0;
|
||||
log174[3 * i + 2] = (b3 & 1) ? +1.0 : -1.0;
|
||||
}
|
||||
// 3140652 00000000117217624541127053545 3140652 33170166234757420515470163426 3140652
|
||||
// 0000000011721762454112705354533170166234757420515470163426
|
||||
// 0000000011721762454112705454544170166344757430515470073537
|
||||
// 0000000011711761444111704343433170166233747320414370072427
|
||||
// 0000000011711761454111705353533170166233757320515370072527
|
||||
}
|
||||
|
||||
|
||||
void print_tones(const uint8_t *code_map, const float *log174) {
|
||||
for (int k = 0; k < 3 * FT8_ND; k += 3) {
|
||||
uint8_t max = 0;
|
||||
|
|
@ -358,14 +160,19 @@ int main(int argc, char **argv) {
|
|||
uint8_t power[num_blocks * 4 * num_bins];
|
||||
extract_power(signal, num_blocks, num_bins, power);
|
||||
|
||||
Candidate heap[kMax_candidates];
|
||||
// Find top candidates by Costas sync score and localize them in time and frequency
|
||||
Candidate candidate_list[kMax_candidates];
|
||||
int num_candidates = find_sync(power, num_blocks, num_bins, kCostas_map, kMax_candidates, candidate_list);
|
||||
|
||||
// TODO: sort the candidates by strongest sync first?
|
||||
|
||||
// Go over candidates and attempt to decode messages
|
||||
char decoded[kMax_decoded_messages][kMax_message_length];
|
||||
int num_decoded = 0;
|
||||
|
||||
int num_candidates = find_sync(power, num_blocks, num_bins, kCostas_map, kMax_candidates, heap);
|
||||
|
||||
for (int idx = 0; idx < num_candidates; ++idx) {
|
||||
Candidate &cand = heap[idx];
|
||||
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 log174[FT8_N];
|
||||
extract_likelihood(power, num_bins, cand, kGray_map, log174);
|
||||
|
|
@ -374,41 +181,22 @@ int main(int argc, char **argv) {
|
|||
uint8_t plain[FT8_N];
|
||||
int n_errors = 0;
|
||||
bp_decode(log174, kLDPC_iterations, plain, &n_errors);
|
||||
//ldpc_decode(log174, num_iters, plain, &n_errors);
|
||||
//ldpc_decode(log174, kLDPC_iterations, plain, &n_errors);
|
||||
|
||||
if (n_errors > 0) {
|
||||
//printf("ldpc_decode() = %d\n", n_errors);
|
||||
continue;
|
||||
}
|
||||
|
||||
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;
|
||||
|
||||
// printf("%03d: score = %d freq = %.1f time = %.2f\n", idx,
|
||||
// cand.score, freq_hz, time_sec);
|
||||
|
||||
//print_tones(kGray_map, log174);
|
||||
|
||||
// Extract payload + CRC
|
||||
// Extract payload + CRC (first FT8_K bits)
|
||||
uint8_t a91[12];
|
||||
uint8_t mask = 0x80;
|
||||
int byte_idx = 0;
|
||||
for (int i = 0; i < 12; ++i) {
|
||||
a91[i] = 0;
|
||||
}
|
||||
for (int i = 0; i < FT8_K; ++i) {
|
||||
if (plain[i]) {
|
||||
a91[byte_idx] |= mask;
|
||||
}
|
||||
mask >>= 1;
|
||||
if (!mask) {
|
||||
mask = 0x80;
|
||||
++byte_idx;
|
||||
}
|
||||
}
|
||||
pack_bits(plain, FT8_K, a91);
|
||||
|
||||
// TODO: check CRC
|
||||
|
||||
// printf("%03d: score = %d freq = %.1f time = %.2f\n", idx,
|
||||
// cand.score, freq_hz, time_sec);
|
||||
// print_tones(kGray_map, log174);
|
||||
// for (int i = 0; i < 12; ++i) {
|
||||
// printf("%02x ", a91[i]);
|
||||
// }
|
||||
|
|
@ -417,7 +205,7 @@ int main(int argc, char **argv) {
|
|||
char message[kMax_message_length];
|
||||
unpack77(a91, message);
|
||||
|
||||
// Check for duplicate messages
|
||||
// Check for duplicate messages (TODO: use hashing)
|
||||
bool found = false;
|
||||
for (int i = 0; i < num_decoded; ++i) {
|
||||
if (0 == strcmp(decoded[i], message)) {
|
||||
|
|
@ -432,7 +220,7 @@ int main(int argc, char **argv) {
|
|||
|
||||
// Fake WSJT-X-like output for now
|
||||
int snr = 0; // TODO: compute SNR
|
||||
printf("000000 %3d %4.1f %4d ~ %s\n", snr, time_sec, (int)(freq_hz + 0.5f), message);
|
||||
printf("000000 %3d %4.1f %4d ~ %s\n", idx, time_sec, (int)(freq_hz + 0.5f), message);
|
||||
}
|
||||
}
|
||||
printf("Decoded %d messages\n", num_decoded);
|
||||
|
|
|
|||
|
|
@ -0,0 +1,243 @@
|
|||
#include "decode.h"
|
||||
|
||||
#include <math.h>
|
||||
|
||||
#include "constants.h"
|
||||
|
||||
static float max2(float a, float b);
|
||||
static float max4(float a, float b, float c, float d);
|
||||
static void heapify_down(Candidate *heap, int heap_size);
|
||||
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);
|
||||
|
||||
|
||||
// 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 heap_size = 0;
|
||||
|
||||
// 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;
|
||||
|
||||
// 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;
|
||||
|
||||
int offset = ((time_offset + k + m) * 4 + alt) * num_bins + freq_offset;
|
||||
const uint8_t *p8 = power + offset;
|
||||
|
||||
score += 8 * p8[sync_map[k]] -
|
||||
p8[0] - p8[1] - p8[2] - p8[3] -
|
||||
p8[4] - p8[5] - p8[6] - p8[7];
|
||||
|
||||
// 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;
|
||||
}
|
||||
}
|
||||
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;
|
||||
|
||||
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;
|
||||
|
||||
heapify_up(heap, heap_size);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
return heap_size;
|
||||
}
|
||||
|
||||
|
||||
// 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;
|
||||
|
||||
// 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) {
|
||||
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);
|
||||
|
||||
decode_symbol(ps, code_map, bit_idx, log174);
|
||||
}
|
||||
|
||||
// Compute the variance of log174
|
||||
float sum = 0;
|
||||
float sum2 = 0;
|
||||
float inv_n = 1.0f / FT8_N;
|
||||
for (int i = 0; i < FT8_N; ++i) {
|
||||
sum += log174[i];
|
||||
sum2 += log174[i] * log174[i];
|
||||
}
|
||||
float variance = (sum2 - sum * sum * inv_n) * inv_n;
|
||||
|
||||
// Normalize log174 such that sigma = 2.83 (Why? It's in WSJT-X, ft8b.f90)
|
||||
// Seems to be 2.83 = sqrt(8). Experimentally sqrt(16) works better.
|
||||
float norm_factor = sqrtf(16.0f / variance);
|
||||
for (int i = 0; i < FT8_N; ++i) {
|
||||
log174[i] *= norm_factor;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
static float max2(float a, float b) {
|
||||
return (a >= b) ? a : b;
|
||||
}
|
||||
|
||||
|
||||
static float max4(float a, float b, float c, float d) {
|
||||
return max2(max2(a, b), max2(c, d));
|
||||
}
|
||||
|
||||
|
||||
static void heapify_down(Candidate *heap, int heap_size) {
|
||||
// heapify from the root down
|
||||
int current = 0;
|
||||
while (true) {
|
||||
int largest = current;
|
||||
int left = 2 * current + 1;
|
||||
int right = left + 1;
|
||||
|
||||
if (left < heap_size && heap[left].score < heap[largest].score) {
|
||||
largest = left;
|
||||
}
|
||||
if (right < heap_size && heap[right].score < heap[largest].score) {
|
||||
largest = right;
|
||||
}
|
||||
if (largest == current) {
|
||||
break;
|
||||
}
|
||||
|
||||
Candidate tmp = heap[largest];
|
||||
heap[largest] = heap[current];
|
||||
heap[current] = tmp;
|
||||
current = largest;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
static void heapify_up(Candidate *heap, int heap_size) {
|
||||
// heapify from the last node up
|
||||
int current = heap_size - 1;
|
||||
while (current > 0) {
|
||||
int parent = (current - 1) / 2;
|
||||
if (heap[current].score >= heap[parent].score) {
|
||||
break;
|
||||
}
|
||||
|
||||
Candidate tmp = heap[parent];
|
||||
heap[parent] = heap[current];
|
||||
heap[current] = tmp;
|
||||
current = parent;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// Compute unnormalized log likelihood log(p(1) / p(0)) of 3 message bits (1 FSK symbol)
|
||||
static void decode_symbol(const uint8_t *power, const uint8_t *code_map, int bit_idx, float *log174) {
|
||||
// Cleaned up code for the simple case of n_syms==1
|
||||
float s2[8];
|
||||
|
||||
for (int j = 0; j < 8; ++j) {
|
||||
s2[j] = (float)power[code_map[j]];
|
||||
}
|
||||
|
||||
log174[bit_idx + 0] = max4(s2[4], s2[5], s2[6], s2[7]) - max4(s2[0], s2[1], s2[2], s2[3]);
|
||||
log174[bit_idx + 1] = max4(s2[2], s2[3], s2[6], s2[7]) - max4(s2[0], s2[1], s2[4], s2[5]);
|
||||
log174[bit_idx + 2] = max4(s2[1], s2[3], s2[5], s2[7]) - max4(s2[0], s2[2], s2[4], s2[6]);
|
||||
}
|
||||
|
||||
|
||||
// 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, const uint8_t *code_map, 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);
|
||||
|
||||
float s2[n_tones];
|
||||
|
||||
for (int j = 0; j < n_tones; ++j) {
|
||||
int j1 = j & 0x07;
|
||||
if (n_syms == 1) {
|
||||
s2[j] = (float)power[code_map[j1]];
|
||||
continue;
|
||||
}
|
||||
int j2 = (j >> 3) & 0x07;
|
||||
if (n_syms == 2) {
|
||||
s2[j] = (float)power[code_map[j2]];
|
||||
s2[j] += (float)power[code_map[j1] + 4 * num_bins];
|
||||
continue;
|
||||
}
|
||||
int j3 = (j >> 6) & 0x07;
|
||||
s2[j] = (float)power[code_map[j3]];
|
||||
s2[j] += (float)power[code_map[j2] + 4 * num_bins];
|
||||
s2[j] += (float)power[code_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
|
||||
for (int i = 0; i < n_bits; ++i) {
|
||||
if (bit_idx + i >= FT8_N) {
|
||||
// Respect array size
|
||||
break;
|
||||
}
|
||||
|
||||
uint16_t mask = (n_tones >> (i + 1));
|
||||
float max_zero = -1000, max_one = -1000;
|
||||
for (int n = 0; n < n_tones; ++n) {
|
||||
if (n & mask) {
|
||||
max_one = max2(max_one, s2[n]);
|
||||
}
|
||||
else {
|
||||
max_zero = max2(max_zero, s2[n]);
|
||||
}
|
||||
}
|
||||
|
||||
log174[bit_idx + i] = max_one - max_zero;
|
||||
}
|
||||
}
|
||||
|
|
@ -0,0 +1,21 @@
|
|||
#pragma once
|
||||
|
||||
#include <stdint.h>
|
||||
|
||||
struct Candidate {
|
||||
int16_t score;
|
||||
int16_t time_offset;
|
||||
int16_t freq_offset;
|
||||
uint8_t time_sub;
|
||||
uint8_t freq_sub;
|
||||
};
|
||||
|
||||
|
||||
// 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);
|
||||
|
||||
|
||||
// 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);
|
||||
28
ft8/ldpc.cpp
28
ft8/ldpc.cpp
|
|
@ -19,6 +19,29 @@ float fast_tanh(float x);
|
|||
float fast_atanh(float x);
|
||||
|
||||
|
||||
// Packs a string of bits each represented as a zero/non-zero byte in plain[],
|
||||
// as a string of packed bits starting from the MSB of the first byte of packed[]
|
||||
void pack_bits(const uint8_t plain[], int num_bits, uint8_t packed[]) {
|
||||
int num_bytes = (num_bits + 7) / 8;
|
||||
for (int i = 0; i < num_bytes; ++i) {
|
||||
packed[i] = 0;
|
||||
}
|
||||
|
||||
uint8_t mask = 0x80;
|
||||
int byte_idx = 0;
|
||||
for (int i = 0; i < num_bits; ++i) {
|
||||
if (plain[i]) {
|
||||
packed[byte_idx] |= mask;
|
||||
}
|
||||
mask >>= 1;
|
||||
if (!mask) {
|
||||
mask = 0x80;
|
||||
++byte_idx;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// codeword is 174 log-likelihoods.
|
||||
// plain is a return value, 174 ints, to be 0 or 1.
|
||||
// max_iters is how hard to try.
|
||||
|
|
@ -95,11 +118,10 @@ void ldpc_decode(float codeword[], int max_iters, uint8_t plain[], int *ok) {
|
|||
int ldpc_check(uint8_t codeword[]) {
|
||||
int errors = 0;
|
||||
|
||||
// kNm[87][7]
|
||||
for (int j = 0; j < FT8_M; ++j) {
|
||||
uint8_t x = 0;
|
||||
for (int ii1 = 0; ii1 < kNrw[j]; ++ii1) {
|
||||
x ^= codeword[kNm[j][ii1] - 1];
|
||||
for (int i = 0; i < kNrw[j]; ++i) {
|
||||
x ^= codeword[kNm[j][i] - 1];
|
||||
}
|
||||
if (x != 0) {
|
||||
++errors;
|
||||
|
|
|
|||
|
|
@ -7,3 +7,7 @@
|
|||
void ldpc_decode(float codeword[], int max_iters, uint8_t plain[], int *ok);
|
||||
|
||||
void bp_decode(float codeword[], int max_iters, uint8_t plain[], int *ok);
|
||||
|
||||
// Packs a string of bits each represented as a zero/non-zero byte in plain[],
|
||||
// as a string of packed bits starting from the MSB of the first byte of packed[]
|
||||
void pack_bits(const uint8_t plain[], int num_bits, uint8_t packed[]);
|
||||
|
|
|
|||
|
|
@ -16,16 +16,6 @@ const char A3[] = "0123456789";
|
|||
const char A4[] = " ABCDEFGHIJKLMNOPQRSTUVWXYZ";
|
||||
|
||||
|
||||
int index(const char *string, char c) {
|
||||
for (int i = 0; *string; ++i, ++string) {
|
||||
if (c == *string) {
|
||||
return i;
|
||||
}
|
||||
}
|
||||
return -1; // Not found
|
||||
}
|
||||
|
||||
|
||||
// Pack a special token, a 22-bit hash code, or a valid base call
|
||||
// into a 28-bit integer.
|
||||
int32_t pack28(const char *callsign) {
|
||||
|
|
@ -81,9 +71,9 @@ int32_t pack28(const char *callsign) {
|
|||
|
||||
// Check for standard callsign
|
||||
int i0, i1, i2, i3, i4, i5;
|
||||
if ((i0 = index(A1, c6[0])) >= 0 && (i1 = index(A2, c6[1])) >= 0 &&
|
||||
(i2 = index(A3, c6[2])) >= 0 && (i3 = index(A4, c6[3])) >= 0 &&
|
||||
(i4 = index(A4, c6[4])) >= 0 && (i5 = index(A4, c6[5])) >= 0)
|
||||
if ((i0 = char_index(A1, c6[0])) >= 0 && (i1 = char_index(A2, c6[1])) >= 0 &&
|
||||
(i2 = char_index(A3, c6[2])) >= 0 && (i3 = char_index(A4, c6[3])) >= 0 &&
|
||||
(i4 = char_index(A4, c6[4])) >= 0 && (i5 = char_index(A4, c6[5])) >= 0)
|
||||
{
|
||||
//printf("Pack28: idx=[%d, %d, %d, %d, %d, %d]\n", i0, i1, i2, i3, i4, i5);
|
||||
// This is a standard callsign
|
||||
|
|
@ -125,7 +115,7 @@ bool chkcall(const char *call, char *bc) {
|
|||
|
||||
// TODO: implement suffix parsing (or rework?)
|
||||
//bc=w(1:6)
|
||||
//i0=index(w,'/')
|
||||
//i0=char_index(w,'/')
|
||||
//if(max(i0-1,n1-i0).gt.6) go to 100 !Base call must be < 7 characters
|
||||
//if(i0.ge.2 .and. i0.le.n1-1) then !Extract base call from compound call
|
||||
// if(i0-1.le.n1-i0) bc=w(i0+1:n1)//' '
|
||||
|
|
@ -206,9 +196,9 @@ int pack77_1(const char *msg, uint8_t *b77) {
|
|||
uint8_t i3 = 1; // No suffix or /R
|
||||
|
||||
// TODO: check for suffixes
|
||||
// if(index(w(1),'/P').ge.4 .or. index(w(2),'/P').ge.4) i3=2 !Type 2, with "/P"
|
||||
// if(index(w(1),'/P').ge.4 .or. index(w(1),'/R').ge.4) ipa=1
|
||||
// if(index(w(2),'/P').ge.4 .or. index(w(2),'/R').ge.4) ipb=1
|
||||
// if(char_index(w(1),'/P').ge.4 .or. char_index(w(2),'/P').ge.4) i3=2 !Type 2, with "/P"
|
||||
// if(char_index(w(1),'/P').ge.4 .or. char_index(w(1),'/R').ge.4) ipa=1
|
||||
// if(char_index(w(2),'/P').ge.4 .or. char_index(w(2),'/R').ge.4) ipb=1
|
||||
|
||||
// Shift in ipa and ipb bits into n28a and n28b
|
||||
n28a <<= 1; // ipa = 0
|
||||
|
|
@ -260,7 +250,7 @@ void packtext77(const char *text, uint8_t *b71) {
|
|||
|
||||
// Add index of the current char
|
||||
if (j < length) {
|
||||
int q = index(A0, text[j]);
|
||||
int q = char_index(A0, text[j]);
|
||||
x = (q > 0) ? q : 0;
|
||||
}
|
||||
else {
|
||||
|
|
|
|||
10
ft8/text.cpp
10
ft8/text.cpp
|
|
@ -33,6 +33,16 @@ bool equals(const char *string1, const char *string2) {
|
|||
}
|
||||
|
||||
|
||||
int char_index(const char *string, char c) {
|
||||
for (int i = 0; *string; ++i, ++string) {
|
||||
if (c == *string) {
|
||||
return i;
|
||||
}
|
||||
}
|
||||
return -1; // Not found
|
||||
}
|
||||
|
||||
|
||||
// Text message formatting:
|
||||
// - replaces lowercase letters with uppercase
|
||||
// - merges consecutive spaces into single space
|
||||
|
|
|
|||
|
|
@ -8,6 +8,8 @@ bool in_range(char c, char min, char max);
|
|||
bool starts_with(const char *string, const char *prefix);
|
||||
bool equals(const char *string1, const char *string2);
|
||||
|
||||
int char_index(const char *string, char c);
|
||||
|
||||
// Text message formatting:
|
||||
// - replaces lowercase letters with uppercase
|
||||
// - merges consecutive spaces into single space
|
||||
|
|
|
|||
|
|
@ -38,7 +38,7 @@ char charn(int c, int table_idx) {
|
|||
|
||||
// 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, char *result) {
|
||||
int unpack28(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) {
|
||||
if (n28 <= 2) {
|
||||
|
|
@ -117,8 +117,17 @@ int unpack28(uint32_t n28, char *result) {
|
|||
while (callsign[ws_len] == ' ') {
|
||||
ws_len++;
|
||||
}
|
||||
|
||||
strcpy(result, callsign + ws_len);
|
||||
|
||||
// Check if we should append /R or /P suffix
|
||||
if (ip) {
|
||||
if (i3 == 1) {
|
||||
strcat(result, "/R");
|
||||
}
|
||||
else if (i3 == 2) {
|
||||
strcat(result, "/P");
|
||||
}
|
||||
}
|
||||
return 0; // Success
|
||||
}
|
||||
|
||||
|
|
@ -148,10 +157,10 @@ int unpack_type1(const uint8_t *a77, uint8_t i3, char *message) {
|
|||
// Unpack both callsigns
|
||||
char field_1[14];
|
||||
char field_2[14];
|
||||
if (unpack28(n28a >> 1, field_1) < 0) {
|
||||
if (unpack28(n28a >> 1, n28a & 0x01, i3, field_1) < 0) {
|
||||
return -1;
|
||||
}
|
||||
if (unpack28(n28b >> 1, field_2) < 0) {
|
||||
if (unpack28(n28b >> 1, n28b & 0x01, i3, field_2) < 0) {
|
||||
return -2;
|
||||
}
|
||||
// Fix "CQ_" to "CQ " -> already done in unpack28()
|
||||
|
|
@ -164,19 +173,13 @@ int unpack_type1(const uint8_t *a77, uint8_t i3, char *message) {
|
|||
strcat(message, " ");
|
||||
strcat(message, field_2);
|
||||
|
||||
// TODO: add /P and /R suffixes
|
||||
// if(index(field_1,'<').le.0) then
|
||||
// ws_len=index(field_1,' ')
|
||||
// if(ws_len.ge.4 .and. ipa.eq.1 .and. i3.eq.1) field_1(ws_len:ws_len+1)='/R'
|
||||
// if(ws_len.ge.4 .and. ipa.eq.1 .and. i3.eq.2) field_1(ws_len:ws_len+1)='/P'
|
||||
// if(ws_len.ge.4) call add_call_to_recent_calls(field_1)
|
||||
// endif
|
||||
// if(index(field_2,'<').le.0) then
|
||||
// ws_len=index(field_2,' ')
|
||||
// if(ws_len.ge.4 .and. ipb.eq.1 .and. i3.eq.1) field_2(ws_len:ws_len+1)='/R'
|
||||
// if(ws_len.ge.4 .and. ipb.eq.1 .and. i3.eq.2) field_2(ws_len:ws_len+1)='/P'
|
||||
// if(ws_len.ge.4) call add_call_to_recent_calls(field_2)
|
||||
// endif
|
||||
// TODO: add to recent calls
|
||||
if (field_1[0] != '<' && strlen(field_1) >= 4) {
|
||||
// add_call_to_recent_calls(field_1)
|
||||
}
|
||||
if (field_2[0] != '<' && strlen(field_2) >= 4) {
|
||||
// add_call_to_recent_calls(field_2)
|
||||
}
|
||||
|
||||
char field_3[5];
|
||||
if (igrid4 <= MAXGRID4) {
|
||||
|
|
@ -287,16 +290,16 @@ int unpack77(const uint8_t *a77, char *message) {
|
|||
// 0.0 Free text
|
||||
return unpack_text(a77, message);
|
||||
}
|
||||
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 == 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, message);
|
||||
|
|
@ -305,17 +308,21 @@ int unpack77(const uint8_t *a77, char *message) {
|
|||
// Type 1 (standard message) or Type 2 ("/P" form for EU VHF contest)
|
||||
return unpack_type1(a77, i3, message);
|
||||
}
|
||||
else if (i3 == 3) {
|
||||
// Type 3: ARRL RTTY Contest
|
||||
}
|
||||
else if (i3 == 4) {
|
||||
// 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.
|
||||
// else if (i3 == 3) {
|
||||
// // Type 3: ARRL RTTY Contest
|
||||
// }
|
||||
// else if (i3 == 4) {
|
||||
// // 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.
|
||||
|
||||
// TODO: implement
|
||||
// read(c77,1050) n12,n58,iflip,nrpt,icq
|
||||
// 1050 format(b12,b58,b1,b2,b1)
|
||||
// // TODO: implement
|
||||
// // read(c77,1050) n12,n58,iflip,nrpt,icq
|
||||
// // 1050 format(b12,b58,b1,b2,b1)
|
||||
// }
|
||||
else {
|
||||
// unknown type
|
||||
message[0] = '\0';
|
||||
}
|
||||
|
||||
return 0;
|
||||
|
|
|
|||
14
test.cpp
14
test.cpp
|
|
@ -9,6 +9,7 @@
|
|||
#include "ft8/v1/encode.h"
|
||||
#include "ft8/pack_v2.h"
|
||||
#include "ft8/encode_v2.h"
|
||||
#include "ft8/constants.h"
|
||||
|
||||
#include "common/debug.h"
|
||||
|
||||
|
|
@ -95,6 +96,19 @@ void test3() {
|
|||
}
|
||||
|
||||
|
||||
void test_tones(float *log174) {
|
||||
// Just a test case
|
||||
for (int i = 0; i < FT8_ND; ++i) {
|
||||
const uint8_t inv_map[8] = {0, 1, 3, 2, 6, 4, 5, 7};
|
||||
uint8_t tone = ("0000000011721762454112705354533170166234757420515470163426"[i]) - '0';
|
||||
uint8_t b3 = inv_map[tone];
|
||||
log174[3 * i] = (b3 & 4) ? +1.0 : -1.0;
|
||||
log174[3 * i + 1] = (b3 & 2) ? +1.0 : -1.0;
|
||||
log174[3 * i + 2] = (b3 & 1) ? +1.0 : -1.0;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
int main() {
|
||||
test1();
|
||||
|
||||
|
|
|
|||
Ładowanie…
Reference in New Issue