From 0c4efc2d8792693c7f9b1ac3873c36e6e117f8d7 Mon Sep 17 00:00:00 2001 From: Karlis Goba Date: Thu, 27 Dec 2018 22:33:07 +0200 Subject: [PATCH] Moved decoding DSP routines to separate file --- Makefile | 2 +- decode_ft8.cpp | 266 +++++----------------------------------------- ft8/decode.cpp | 243 ++++++++++++++++++++++++++++++++++++++++++ ft8/decode.h | 21 ++++ ft8/ldpc.cpp | 28 ++++- ft8/ldpc.h | 4 + ft8/pack_v2.cpp | 26 ++--- ft8/text.cpp | 10 ++ ft8/text.h | 2 + ft8/unpack_v2.cpp | 81 +++++++------- test.cpp | 14 +++ 11 files changed, 399 insertions(+), 298 deletions(-) create mode 100644 ft8/decode.cpp create mode 100644 ft8/decode.h diff --git a/Makefile b/Makefile index bfa5a1b..1046037 100644 --- a/Makefile +++ b/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: diff --git a/decode_ft8.cpp b/decode_ft8.cpp index 9c13390..d43c340 100644 --- a/decode_ft8.cpp +++ b/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); diff --git a/ft8/decode.cpp b/ft8/decode.cpp new file mode 100644 index 0000000..13f0b2b --- /dev/null +++ b/ft8/decode.cpp @@ -0,0 +1,243 @@ +#include "decode.h" + +#include + +#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; + } +} \ No newline at end of file diff --git a/ft8/decode.h b/ft8/decode.h new file mode 100644 index 0000000..03ee73c --- /dev/null +++ b/ft8/decode.h @@ -0,0 +1,21 @@ +#pragma once + +#include + +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); diff --git a/ft8/ldpc.cpp b/ft8/ldpc.cpp index 8441e22..205f508 100644 --- a/ft8/ldpc.cpp +++ b/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; diff --git a/ft8/ldpc.h b/ft8/ldpc.h index 370e0f7..56866da 100644 --- a/ft8/ldpc.h +++ b/ft8/ldpc.h @@ -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[]); diff --git a/ft8/pack_v2.cpp b/ft8/pack_v2.cpp index 7681670..41e9dad 100644 --- a/ft8/pack_v2.cpp +++ b/ft8/pack_v2.cpp @@ -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 { diff --git a/ft8/text.cpp b/ft8/text.cpp index ec0fe60..7f52184 100644 --- a/ft8/text.cpp +++ b/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 diff --git a/ft8/text.h b/ft8/text.h index 2415eb6..7218050 100644 --- a/ft8/text.h +++ b/ft8/text.h @@ -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 diff --git a/ft8/unpack_v2.cpp b/ft8/unpack_v2.cpp index a8a96b7..653043f 100644 --- a/ft8/unpack_v2.cpp +++ b/ft8/unpack_v2.cpp @@ -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 -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 -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. 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. 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; diff --git a/test.cpp b/test.cpp index 85a5866..07b8afb 100644 --- a/test.cpp +++ b/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();