diff --git a/decode_ft8.cpp b/decode_ft8.cpp index 724805f..abd1b1b 100644 --- a/decode_ft8.cpp +++ b/decode_ft8.cpp @@ -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); diff --git a/ft8/decode.cpp b/ft8/decode.cpp index 02ae2c5..42ef2a4 100644 --- a/ft8/decode.cpp +++ b/ft8/decode.cpp @@ -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); } diff --git a/ft8/decode.h b/ft8/decode.h index 52211d9..05f05e6 100644 --- a/ft8/decode.h +++ b/ft8/decode.h @@ -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); } \ No newline at end of file diff --git a/ft8/unpack.cpp b/ft8/unpack.cpp index 8ab061a..35eb783 100644 --- a/ft8/unpack.cpp +++ b/ft8/unpack.cpp @@ -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 }