From e6cbd139ff0510c6004b2657b6e1579c06175993 Mon Sep 17 00:00:00 2001 From: Karlis Goba Date: Mon, 24 Dec 2018 14:22:26 +0200 Subject: [PATCH] Updated decoding to v2 --- Makefile | 11 +- decode_ft8.cpp | 136 ++++++++++------- ft8/ldpc.cpp | 373 +++++++++++++++++++++++----------------------- ft8/ldpc.h | 4 +- ft8/pack_v2.cpp | 64 ++++---- ft8/text.cpp | 4 +- ft8/text.h | 2 +- ft8/unpack_v2.cpp | 322 +++++++++++++++++++++++++++++++++++++++ ft8/unpack_v2.h | 7 + ft8/v1/pack.cpp | 2 +- ft8/v1/unpack.cpp | 2 +- gen_ft8.cpp | 10 +- test.cpp | 6 +- 13 files changed, 658 insertions(+), 285 deletions(-) create mode 100644 ft8/unpack_v2.cpp create mode 100644 ft8/unpack_v2.h diff --git a/Makefile b/Makefile index 14700de..bfa5a1b 100644 --- a/Makefile +++ b/Makefile @@ -1,8 +1,7 @@ CXXFLAGS = -std=c++14 -I. LDFLAGS = -lm -TARGETS = gen_ft8 -#TARGETS = gen_ft8 decode_ft8 test +TARGETS = gen_ft8 decode_ft8 test .PHONY: run_tests all clean @@ -14,11 +13,11 @@ run_tests: test gen_ft8: gen_ft8.o ft8/constants.o ft8/text.o ft8/pack_v2.o ft8/encode_v2.o common/wave.o $(CXX) $(LDFLAGS) -o $@ $^ -#test: test.o ft8/encode.o ft8/pack.o ft8/text.o ft8/pack_v2.o ft8/encode_v2.o ft8/unpack.o -# $(CXX) $(LDFLAGS) -o $@ $^ +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.o ft8/text.o common/wave.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 + $(CXX) $(LDFLAGS) -o $@ $^ clean: rm -f *.o ft8/*.o common/*.o fft/*.o $(TARGETS) diff --git a/decode_ft8.cpp b/decode_ft8.cpp index 84b8bc7..5c850f3 100644 --- a/decode_ft8.cpp +++ b/decode_ft8.cpp @@ -3,16 +3,12 @@ #include #include -#include "common/wave.h" -#include "ft8/pack.h" -#include "ft8/encode.h" -#include "ft8/pack_v2.h" -#include "ft8/encode_v2.h" -#include "ft8/unpack.h" - +#include "ft8/unpack_v2.h" #include "ft8/ldpc.h" -#include "fft/kiss_fftr.h" +#include "ft8/constants.h" +#include "common/wave.h" +#include "fft/kiss_fftr.h" void usage() { printf("Decode a 15-second WAV file.\n"); @@ -79,11 +75,11 @@ void heapify_up(Candidate * heap, int heap_size) { // 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). -void 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 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 = 0; time_offset < num_blocks - NN; ++time_offset) { + for (int time_offset = 0; time_offset < num_blocks - FT8_NN; ++time_offset) { for (int freq_offset = 0; freq_offset < num_bins - 8; ++freq_offset) { int score = 0; @@ -122,6 +118,8 @@ void find_sync(const uint8_t *power, int num_blocks, int num_bins, const uint8_t } } } + + return heap_size; } @@ -134,6 +132,9 @@ void extract_power(const float * signal, int num_blocks, int num_bins, uint8_t * for (int i = 0; i < nfft; ++i) { window[i] = hann_i(i, nfft); } + // for (int i = 0; i < nfft; ++i) { + // window[i] = (i < block_size) ? 2 * hann_i(i, block_size) : 0.0f; + // } size_t fft_work_size; kiss_fftr_alloc(nfft, 0, 0, &fft_work_size); @@ -141,11 +142,12 @@ void extract_power(const float * signal, int num_blocks, int num_bins, uint8_t * 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); int offset = 0; float fft_norm = 1.0f / 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) for (int time_sub = 0; time_sub <= block_size/2; time_sub += block_size/2) { @@ -162,8 +164,8 @@ void extract_power(const float * signal, int num_blocks, int num_bins, uint8_t * // Compute log magnitude in decibels for (int j = 0; j < nfft/2 + 1; ++j) { - float mag2 = fft_norm * (freqdata[j].i * freqdata[j].i + freqdata[j].r * freqdata[j].r); - mag_db[j] = 10.0f * log10f(1.0E-10f + mag2); + 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); } // Loop over two possible frequency bin offsets (for averaging) @@ -173,15 +175,19 @@ void extract_power(const float * signal, int num_blocks, int num_bins, uint8_t * float db2 = mag_db[j * 2 + freq_sub + 1]; float db = (db1 + db2) / 2; - // Scale decibels to unsigned 8-bit range + // Scale decibels to unsigned 8-bit range and clamp the value int scaled = (int)(2 * (db + 100)); power[offset] = (scaled < 0) ? 0 : ((scaled > 255) ? 255 : scaled); ++offset; + + if (db > max_mag) max_mag = db; } } } } + printf("Max magnitude: %.1f dB\n", max_mag); + free(fft_work); } @@ -203,7 +209,7 @@ void extract_likelihood(const uint8_t *power, int num_bins, const Candidate & ca int k = 0; // Go over FSK tones and skip Costas sync symbols - for (int i = 7; i < NN - 7; ++i) { + for (int i = 7; i < FT8_NN - 7; ++i) { if (i == 36) i += 7; // Pointer to 8 bins of the current symbol @@ -222,14 +228,15 @@ void extract_likelihood(const uint8_t *power, int num_bins, const Candidate & ca // printf("%d %d %d %d %d %d %d %d : %.0f %.0f %.0f\n", // ps[0], ps[1], ps[2], ps[3], ps[4], ps[5], ps[6], ps[7], // log174[k + 0], log174[k + 1], log174[k + 2]); + k += 3; } // Compute the variance of log174 float sum = 0; float sum2 = 0; - float inv_n = 1.0f / (3 * ND); - for (int i = 0; i < 3 * ND; ++i) { + float inv_n = 1.0f / (3 * FT8_ND); + for (int i = 0; i < 3 * FT8_ND; ++i) { sum += log174[i]; sum2 += log174[i] * log174[i]; } @@ -238,7 +245,7 @@ void extract_likelihood(const uint8_t *power, int num_bins, const Candidate & ca // Normalize log174 such that sigma = 2.83 (Why? It's in WSJT-X) float norm_factor = 2.83f / sqrtf(variance); - for (int i = 0; i < 3 * ND; ++i) { + for (int i = 0; i < 3 * FT8_ND; ++i) { log174[i] *= norm_factor; //printf("%.1f ", log174[i]); } @@ -246,6 +253,35 @@ void extract_likelihood(const uint8_t *power, int num_bins, const Candidate & ca } +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; + if (log174[k + 0] > 0) max |= 4; + if (log174[k + 1] > 0) max |= 2; + if (log174[k + 2] > 0) max |= 1; + printf("%d", code_map[max]); + } + printf("\n"); +} + + int main(int argc, char ** argv) { // Expect one command-line argument if (argc < 2) { @@ -264,13 +300,6 @@ int main(int argc, char ** argv) { return -1; } - // Costas 7x7 tone pattern - const uint8_t kCostas_map_v1[] = { 2,5,6,0,4,1,3 }; - const uint8_t kCostas_map_v2[] = { 3,1,4,0,6,5,2 }; - // Gray maps (used only in v2) - const uint8_t kGray_map_v1[8] = { 0,1,2,3,4,5,6,7 }; // identity map - const uint8_t kGray_map_v2[8] = { 0,1,3,2,5,6,4,7 }; - const float fsk_dev = 6.25f; const int num_bins = (int)(sample_rate / (2 * fsk_dev)); @@ -282,41 +311,45 @@ int main(int argc, char ** argv) { extract_power(signal, num_blocks, num_bins, power); - int num_candidates = 250; + int num_candidates = 100; Candidate heap[num_candidates]; - find_sync(power, num_blocks, num_bins, kCostas_map_v1, num_candidates, heap); + find_sync(power, num_blocks, num_bins, kCostas_map, num_candidates, heap); for (int idx = 0; idx < num_candidates; ++idx) { Candidate &cand = heap[idx]; - float log174[3 * ND]; - extract_likelihood(power, num_bins, cand, kGray_map_v1, log174); + float log174[3 * FT8_ND]; + extract_likelihood(power, num_bins, cand, kGray_map, log174); - const int num_iters = 20; - int plain[3 * ND]; - int ok = 0; + const int num_iters = 25; + uint8_t plain[3 * FT8_ND]; + int n_errors = 0; - bp_decode(log174, num_iters, plain, &ok); - //ldpc_decode(log174, num_iters, plain, &ok); - //printf("ldpc_decode() = %d\n", ok); - - if (ok == 87) { 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); - uint8_t a87[11]; + printf("%03d: score = %d freq = %.1f time = %.2f\n", idx, + cand.score, freq_hz, time_sec); + + bp_decode(log174, num_iters, plain, &n_errors); + //ldpc_decode(log174, num_iters, plain, &n_errors); + printf("ldpc_decode() = %d\n", n_errors); + + if (n_errors == 0) { + + //print_tones(kGray_map, log174); + + // Extract payload + CRC + uint8_t a91[12]; uint8_t mask = 0x80; uint8_t position = 0; - for (int i = 0; i < 11; ++i) { - a87[i] = 0; + for (int i = 0; i < 12; ++i) { + a91[i] = 0; } - // Extract payload + CRC (last 87 bits) - for (int i = 174 - 87; i < 174; ++i) { + for (int i = 0; i < FT8_K; ++i) { if (plain[i]) { - a87[position] |= mask; + a91[position] |= mask; } mask >>= 1; if (!mask) { @@ -325,14 +358,17 @@ int main(int argc, char ** argv) { } } - for (int i = 0; i < 11; ++i) { - //printf("%02x ", a87[i]); - } - //printf("\n"); + // TODO: check CRC + + // for (int i = 0; i < 12; ++i) { + // printf("%02x ", a91[i]); + // } + // printf("\n"); char message[20]; - unpack(a87, message); + unpack77(a91, message); + // Fake WSJT-X-like output for now printf("000000 0 %4.1f %4d ~ %s\n", time_sec, (int)(freq_hz + 0.5f), message); } } diff --git a/ft8/ldpc.cpp b/ft8/ldpc.cpp index 6f89897..f4bc2e4 100644 --- a/ft8/ldpc.cpp +++ b/ft8/ldpc.cpp @@ -14,9 +14,190 @@ #include #include "constants.h" -int ldpc_check(int codeword[]); +int ldpc_check(uint8_t codeword[]); +float fast_tanh(float x); +float fast_atanh(float x); +// codeword is 174 log-likelihoods. +// plain is a return value, 174 ints, to be 0 or 1. +// max_iters is how hard to try. +// ok == 87 means success. +void ldpc_decode(float codeword[], int max_iters, uint8_t plain[], int *ok) { + float m[FT8_M][FT8_N]; // ~60 kB + float e[FT8_M][FT8_N]; // ~60 kB + int min_errors = FT8_M; + + for (int j = 0; j < FT8_M; j++) { + for (int i = 0; i < FT8_N; i++) { + m[j][i] = codeword[i]; + e[j][i] = 0.0f; + } + } + + for (int iter = 0; iter < max_iters; iter++) { + for (int j = 0; j < FT8_M; j++) { + for (int ii1 = 0; ii1 < kNrw[j]; ii1++) { + int i1 = kNm[j][ii1] - 1; + float a = 1.0f; + for (int ii2 = 0; ii2 < kNrw[j]; ii2++) { + int i2 = kNm[j][ii2] - 1; + if (i2 != i1) { + a *= fast_tanh(-m[j][i2] / 2.0f); + } + } + e[j][i1] = logf((1 - a) / (1 + a)); + } + } + + uint8_t cw[FT8_N]; + for (int i = 0; i < FT8_N; i++) { + float l = codeword[i]; + for (int j = 0; j < 3; j++) + l += e[kMn[i][j] - 1][i]; + cw[i] = (l > 0) ? 1 : 0; + } + + int errors = ldpc_check(cw); + + if (errors < min_errors) { + // Update the current best result + for (int i = 0; i < FT8_N; i++) { + plain[i] = cw[i]; + } + min_errors = errors; + + if (errors == 0) { + break; // Found a perfect answer + } + } + + for (int i = 0; i < FT8_N; i++) { + for (int ji1 = 0; ji1 < 3; ji1++) { + int j1 = kMn[i][ji1] - 1; + float l = codeword[i]; + for (int ji2 = 0; ji2 < 3; ji2++) { + if (ji1 != ji2) { + int j2 = kMn[i][ji2] - 1; + l += e[j2][i]; + } + } + m[j1][i] = l; + } + } + } + + *ok = min_errors; +} + + +// +// does a 174-bit codeword pass the FT8's LDPC parity checks? +// returns the number of parity errors. +// 0 means total success. +// +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]; + } + if (x != 0) { + ++errors; + } + } + return errors; +} + + +void bp_decode(float codeword[], int max_iters, uint8_t plain[], int *ok) { + float tov[FT8_N][3]; + float toc[FT8_M][7]; + + int min_errors = FT8_M; + + int nclast = 0; + int ncnt = 0; + + // initialize messages to checks + for (int i = 0; i < FT8_M; ++i) { + for (int j = 0; j < kNrw[i]; ++j) { + toc[i][j] = codeword[kNm[i][j] - 1]; + } + } + + for (int i = 0; i < FT8_N; ++i) { + for (int j = 0; j < 3; ++j) { + tov[i][j] = 0; + } + } + + for (int iter = 0; iter < max_iters; ++iter) { + float zn[FT8_N]; + uint8_t cw[FT8_N]; + + // Update bit log likelihood ratios (tov=0 in iter 0) + for (int i = 0; i < FT8_N; ++i) { + zn[i] = codeword[i] + tov[i][0] + tov[i][1] + tov[i][2]; + cw[i] = (zn[i] > 0) ? 1 : 0; + } + + // Check to see if we have a codeword (check before we do any iter) + int errors = ldpc_check(cw); + + if (errors < min_errors) { + // we have a better guess - update the result + for (int i = 0; i < FT8_N; i++) { + plain[i] = cw[i]; + } + min_errors = errors; + + if (errors == FT8_M) { + break; // Found a perfect answer + } + } + + // Send messages from bits to check nodes + for (int i = 0; i < FT8_M; ++i) { + for (int j = 0; j < kNrw[i]; ++j) { + int ibj = kNm[i][j] - 1; + toc[i][j] = zn[ibj]; + for (int kk = 0; kk < 3; ++kk) { + // subtract off what the bit had received from the check + if (kMn[ibj][kk] - 1 == i) { + toc[i][j] -= tov[ibj][kk]; + } + } + } + } + + // send messages from check nodes to variable nodes + for (int i = 0; i < FT8_M; ++i) { + for (int j = 0; j < kNrw[i]; ++j) { + toc[i][j] = fast_tanh(-toc[i][j] / 2); + } + } + + for (int i = 0; i < FT8_N; ++i) { + for (int j = 0; j < 3; ++j) { + int ichk = kMn[i][j] - 1; // kMn(:,j) are the checks that include bit j + float Tmn = 1.0f; + for (int k = 0; k < kNrw[ichk]; ++k) { + if (kNm[ichk][k] - 1 != i) { + Tmn *= toc[ichk][k]; + } + } + tov[i][j] = 2 * fast_atanh(-Tmn); + } + } + } + + *ok = min_errors; +} + // https://varietyofsound.wordpress.com/2011/02/14/efficient-tanh-computation-using-lamberts-continued-fraction/ // http://functions.wolfram.com/ElementaryFunctions/ArcTanh/10/0001/ // https://mathr.co.uk/blog/2017-09-06_approximating_hyperbolic_tangent.html @@ -95,192 +276,4 @@ float platanh(float x) { return isign * (x - 0.9914f) / 0.0012f; } return isign * 7.0f; -} - - -// codeword is 174 log-likelihoods. -// plain is a return value, 174 ints, to be 0 or 1. -// max_iters is how hard to try. -// ok == 87 means success. -void ldpc_decode(float codeword[], int max_iters, int plain[], int *ok) { - float m[FT8_M][FT8_N]; // ~60 kB - float e[FT8_M][FT8_N]; // ~60 kB - int best_score = -1; - - for (int j = 0; j < FT8_M; j++) { - for (int i = 0; i < FT8_N; i++) { - m[j][i] = codeword[i]; - e[j][i] = 0.0f; - } - } - - for (int iter = 0; iter < max_iters; iter++) { - for (int j = 0; j < FT8_M; j++) { - for (int ii1 = 0; ii1 < kNrw[j]; ii1++) { - int i1 = kNm[j][ii1] - 1; - float a = 1.0f; - for (int ii2 = 0; ii2 < kNrw[j]; ii2++) { - int i2 = kNm[j][ii2] - 1; - if (i2 != i1) { - a *= fast_tanh(-m[j][i2] / 2.0f); - } - } - e[j][i1] = logf((1 - a) / (1 + a)); - } - } - - int cw[FT8_N]; - for (int i = 0; i < FT8_N; i++) { - float l = codeword[i]; - for (int j = 0; j < 3; j++) - l += e[kMn[i][j] - 1][i]; - cw[i] = (l > 0) ? 1 : 0; - } - - int score = ldpc_check(cw); - - if (score > best_score) { - // Update the current best result - for (int i = 0; i < FT8_N; i++) { - //plain[i] = cw[colorder[i]]; // Reverse the column permutation - plain[i] = cw[i]; - } - best_score = score; - - if (score == FT8_M) { - // Found a perfect answer - break; - } - } - - - for (int i = 0; i < FT8_N; i++) { - for (int ji1 = 0; ji1 < 3; ji1++) { - int j1 = kMn[i][ji1] - 1; - float l = codeword[i]; - for (int ji2 = 0; ji2 < 3; ji2++) { - if (ji1 != ji2) { - int j2 = kMn[i][ji2] - 1; - l += e[j2][i]; - } - } - m[j1][i] = l; - } - } - } - - *ok = best_score; -} - - -// -// does a 174-bit codeword pass the FT8's LDPC parity checks? -// returns the number of parity checks that passed. -// 87 means total success. -// -int ldpc_check(int codeword[]) { - int score = 0; - - // kNm[87][7] - for (int j = 0; j < FT8_M; ++j) { - int x = 0; - for (int ii1 = 0; ii1 < kNrw[j]; ++ii1) { - x ^= codeword[kNm[j][ii1] - 1]; - } - if (x == 0) { - ++score; - } - } - return score; -} - - -void bp_decode(float codeword[], int max_iters, int plain[], int *ok) { - float tov[FT8_N][3]; - float toc[FT8_M][7]; - - int best_score = -1; - - int nclast = 0; - int ncnt = 0; - - // initialize messages to checks - for (int i = 0; i < FT8_M; ++i) { - for (int j = 0; j < kNrw[i]; ++j) { - toc[i][j] = codeword[kNm[i][j] - 1]; - } - } - - for (int i = 0; i < FT8_N; ++i) { - for (int j = 0; j < 3; ++j) { - tov[i][j] = 0; - } - } - - for (int iter = 0; iter < max_iters; ++iter) { - float zn[FT8_N]; - int cw[FT8_N]; - - // Update bit log likelihood ratios (tov=0 in iter 0) - for (int i = 0; i < FT8_N; ++i) { - zn[i] = codeword[i] + tov[i][0] + tov[i][1] + tov[i][2]; - cw[i] = (zn[i] > 0) ? 1 : 0; - } - - // Check to see if we have a codeword (check before we do any iter) - int score = ldpc_check(cw); - - if (score > best_score) { - // we have a better guess - update the result - for (int i = 0; i < FT8_N; i++) { - //plain[i] = cw[colorder[i]]; // Reverse the column permutation - plain[i] = cw[i]; - } - best_score = score; - - if (score == FT8_M) { - break; - } - } - - // Send messages from bits to check nodes - for (int i = 0; i < FT8_M; ++i) { - for (int j = 0; j < kNrw[i]; ++j) { - int ibj = kNm[i][j] - 1; - toc[i][j] = zn[ibj]; - for (int kk = 0; kk < 3; ++kk) { - // subtract off what the bit had received from the check - if (kMn[ibj][kk] - 1 == i) { - toc[i][j] -= tov[ibj][kk]; - } - } - } - } - - // send messages from check nodes to variable nodes - for (int i = 0; i < FT8_M; ++i) { - for (int j = 0; j < kNrw[i]; ++j) { - //toc[i][j] = pltanh(-toc[i][j] / 2); - toc[i][j] = fast_tanh(-toc[i][j] / 2); - //toc[i][j] = tanhf(-toc[i][j] / 2); - } - } - - for (int i = 0; i < FT8_N; ++i) { - for (int j = 0; j < 3; ++j) { - int ichk = kMn[i][j] - 1; // kMn(:,j) are the checks that include bit j - float Tmn = 1.0f; - for (int k = 0; k < kNrw[ichk]; ++k) { - if (kNm[ichk][k] - 1 != i) { - Tmn *= toc[ichk][k]; - } - } - //tov[i][j] = 2 * platanh(-Tmn); - tov[i][j] = 2 * fast_atanh(-Tmn); - //tov[i][j] = 2 * atanhf(-Tmn); - } - } - } - - *ok = best_score; -} +} \ No newline at end of file diff --git a/ft8/ldpc.h b/ft8/ldpc.h index 5928b3e..370e0f7 100644 --- a/ft8/ldpc.h +++ b/ft8/ldpc.h @@ -4,6 +4,6 @@ // plain is a return value, 174 ints, to be 0 or 1. // iters is how hard to try. // ok == 87 means success. -void ldpc_decode(float codeword[], int max_iters, int plain[], int *ok); +void ldpc_decode(float codeword[], int max_iters, uint8_t plain[], int *ok); -void bp_decode(float codeword[], int max_iters, int plain[], int *ok); +void bp_decode(float codeword[], int max_iters, uint8_t plain[], int *ok); diff --git a/ft8/pack_v2.cpp b/ft8/pack_v2.cpp index 1db6d98..7681670 100644 --- a/ft8/pack_v2.cpp +++ b/ft8/pack_v2.cpp @@ -9,11 +9,12 @@ namespace ft8_v2 { // TODO: This is wasteful, should figure out something more elegant -const char *A0 = " 0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ+-./?"; -const char *A1 = " 0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ"; -const char *A2 = "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ"; -const char *A3 = "0123456789"; -const char *A4 = " ABCDEFGHIJKLMNOPQRSTUVWXYZ"; +const char A0[] = " 0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ+-./?"; +const char A1[] = " 0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ"; +const char A2[] = "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ"; +const char A3[] = "0123456789"; +const char A4[] = " ABCDEFGHIJKLMNOPQRSTUVWXYZ"; + int index(const char *string, char c) { for (int i = 0; *string; ++i, ++string) { @@ -24,6 +25,7 @@ int index(const char *string, char c) { 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) { @@ -44,8 +46,8 @@ int32_t pack28(const char *callsign) { } // TODO: Check for <...> callsign - // if(c13(1:1).eq.'<')then - // call save_hash_call(c13,n10,n12,n22) !Save callsign in hash table + // if(text(1:1).eq.'<')then + // call save_hash_call(text,n10,n12,n22) !Save callsign in hash table // n28=NTOKENS + n22 char c6[6] = {' ', ' ', ' ', ' ', ' ', ' '}; @@ -55,21 +57,24 @@ int32_t pack28(const char *callsign) { length++; } - // Work-around for Swaziland prefix: + // Copy callsign to 6 character buffer if (starts_with(callsign, "3DA0") && length <= 7) { + // Work-around for Swaziland prefix: 3DA0XYZ -> 3D0XYZ memcpy(c6, "3D0", 3); memcpy(c6 + 3, callsign + 4, length - 4); } - // Work-around for Guinea prefixes: else if (starts_with(callsign, "3X") && is_letter(callsign[2]) && length <= 7) { + // Work-around for Guinea prefixes: 3XA0XYZ -> QA0XYZ memcpy(c6, "Q", 1); memcpy(c6 + 1, callsign + 2, length - 2); } else { if (is_digit(callsign[2]) && length <= 6) { + // AB0XYZ memcpy(c6, callsign, length); } else if (is_digit(callsign[1]) && length <= 5) { + // A0XYZ -> " A0XYZ" memcpy(c6 + 1, callsign, length); } } @@ -80,7 +85,7 @@ int32_t pack28(const char *callsign) { (i2 = index(A3, c6[2])) >= 0 && (i3 = index(A4, c6[3])) >= 0 && (i4 = index(A4, c6[4])) >= 0 && (i5 = index(A4, c6[5])) >= 0) { - printf("Pack28: idx=[%d, %d, %d, %d, %d, %d]\n", i0, i1, i2, i3, i4, i5); + //printf("Pack28: idx=[%d, %d, %d, %d, %d, %d]\n", i0, i1, i2, i3, i4, i5); // This is a standard callsign int32_t n28 = i0; n28 = n28 * 36 + i1; @@ -88,17 +93,17 @@ int32_t pack28(const char *callsign) { n28 = n28 * 27 + i3; n28 = n28 * 27 + i4; n28 = n28 * 27 + i5; - printf("Pack28: n28=%d (%04xh)\n", n28, n28); + //printf("Pack28: n28=%d (%04xh)\n", n28, n28); return NTOKENS + MAX22 + n28; } - //char c13[13]; + //char text[13]; //if (length > 13) return -1; // TODO: // Treat this as a nonstandard callsign: compute its 22-bit hash - // call save_hash_call(c13,n10,n12,n22) !Save callsign in hash table + // call save_hash_call(text,n10,n12,n22) !Save callsign in hash table // n28=NTOKENS + n22 // n28=iand(n28,ishft(1,28)-1) @@ -213,10 +218,6 @@ int pack77_1(const char *msg, uint8_t *b77) { // write(c77,1000) n28a,ipa,n28b,ipb,ir,igrid4,i3 // 1000 format(2(b28.28,b1),b1,b15.15,b3.3) - // 00 00 00 27 b3 00 01 0b 27 8f b9 e0 - // (00 00 00 2) 0 (f0 85 ab e) (0) - // 7842D5F - b77[0] = (n28a >> 21); b77[1] = (n28a >> 13); b77[2] = (n28a >> 5); @@ -232,30 +233,43 @@ int pack77_1(const char *msg, uint8_t *b77) { } -void packtext77(const char *c13, uint8_t *b71) { - // TODO: w=adjustr(c13) +void packtext77(const char *text, uint8_t *b71) { + int length = strlen(text); + + // Skip leading and trailing spaces + while (*text == ' ' && *text != 0) { + ++text; + --length; + } + while (length > 0 && text[length - 1] == ' ') { + --length; + } for (int i = 0; i < 9; ++i) { b71[i] = 0; } for (int j = 0; j < 13; ++j) { - int q = index(A0, c13[j]); - - // Multiply b71 by 42 + // Multiply the long integer in b71 by 42 uint16_t x = 0; for (int i = 8; i >= 0; --i) { x += b71[i] * (uint16_t)42; - b71[i] = x; + b71[i] = (x & 0xFF); x >>= 8; } // Add index of the current char - x = (q > 0) ? q : 0; + if (j < length) { + int q = index(A0, text[j]); + x = (q > 0) ? q : 0; + } + else { + x = 0; + } for (int i = 8; i >= 0; --i) { if (x == 0) break; x += b71[i]; - b71[i] = x; + b71[i] = (x & 0xFF); x >>= 8; } } diff --git a/ft8/text.cpp b/ft8/text.cpp index 3e15ef9..ec0fe60 100644 --- a/ft8/text.cpp +++ b/ft8/text.cpp @@ -78,13 +78,13 @@ int dd_to_int(const char *str, int length) { // Convert a 2 digit integer to string -void int_to_dd(char *str, int value, int width) { +void int_to_dd(char *str, int value, int width, bool full_sign) { if (value < 0) { *str = '-'; ++str; value = -value; } - else { + else if (full_sign) { *str = '+'; ++str; } diff --git a/ft8/text.h b/ft8/text.h index 2891668..2415eb6 100644 --- a/ft8/text.h +++ b/ft8/text.h @@ -17,4 +17,4 @@ void fmtmsg(char *msg_out, const char *msg_in); int dd_to_int(const char *str, int length); // Convert a 2 digit integer to string -void int_to_dd(char *str, int value, int width); +void int_to_dd(char *str, int value, int width, bool full_sign = false); diff --git a/ft8/unpack_v2.cpp b/ft8/unpack_v2.cpp new file mode 100644 index 0000000..a8a96b7 --- /dev/null +++ b/ft8/unpack_v2.cpp @@ -0,0 +1,322 @@ +#include "unpack_v2.h" +#include "text.h" + +#include + +//const uint32_t NBASE = 37L*36L*10L*27L*27L*27L; +const uint32_t MAX22 = 4194304L; +const uint32_t NTOKENS = 2063592L; +const uint16_t MAXGRID4 = 32400L; + + +// convert integer index to ASCII character according to one of 5 tables: +// table 0: " 0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ+-./?" +// table 1: " 0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ" +// table 2: "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ" +// table 3: "0123456789" +// table 4: " ABCDEFGHIJKLMNOPQRSTUVWXYZ" +char charn(int c, int table_idx) { + if (table_idx == 0 || table_idx == 1 || table_idx == 4) { + if (c == 0) return ' '; + c -= 1; + } + if (table_idx != 4) { + if (c < 10) return '0' + c; + c -= 10; + } + if (table_idx != 3) { + if (c < 26) return 'A' + c; + c -= 26; + } + if (table_idx == 0) { + if (c < 5) return "+-./?" [c]; + } + + return '_'; // unknown character, should never get here +} + + +// 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) { + // Check for special tokens DE, QRZ, CQ, CQ_nnn, CQ_aaaa + if (n28 < NTOKENS) { + if (n28 <= 2) { + if (n28 == 0) strcpy(result, "DE"); + if (n28 == 1) strcpy(result, "QRZ"); + if (n28 == 2) strcpy(result, "CQ"); + return 0; // Success + } + if (n28 <= 1002) { + // CQ_nnn with 3 digits + strcpy(result, "CQ "); + int_to_dd(result + 3, n28 - 3, 3); + return 0; // Success + } + if (n28 <= 532443L) { + // CQ_aaaa with 4 alphanumeric symbols + uint32_t n = n28 - 1003; + char aaaa[5]; + + aaaa[4] = '\0'; + aaaa[3] = charn(n % 27, 4); + n /= 27; + aaaa[2] = charn(n % 27, 4); + n /= 27; + aaaa[1] = charn(n % 27, 4); + n /= 27; + aaaa[0] = charn(n % 27, 4); + + // Skip leading whitespace + int ws_len = 0; + while (aaaa[ws_len] == ' ') { + ws_len++; + } + + strcpy(result, "CQ "); + strcat(result, aaaa + ws_len); + return 0; // Success + } + // ? TODO: unspecified in the WSJT-X code + return -1; + } + + n28 = n28 - NTOKENS; + if (n28 < MAX22) { + // This is a 22-bit hash of a result + //n22=n28 + //call hash22(n22,c13) !Retrieve result from hash table + // TODO: implement + return -2; + } + + // Standard callsign + uint32_t n = n28 - MAX22; + + char callsign[7]; + callsign[6] = '\0'; + callsign[5] = charn(n % 27, 4); + n /= 27; + callsign[4] = charn(n % 27, 4); + n /= 27; + callsign[3] = charn(n % 27, 4); + n /= 27; + callsign[2] = charn(n % 10, 3); + n /= 10; + callsign[1] = charn(n % 36, 2); + n /= 36; + callsign[0] = charn(n % 37, 1); + + // Skip trailing and leading whitespace in case of a short callsign + int ws_len = 0; + while (ws_len <= 5 && callsign[5 - ws_len] == ' ') { + callsign[5 - ws_len] = '\0'; + ws_len++; + } + ws_len = 0; + while (callsign[ws_len] == ' ') { + ws_len++; + } + + strcpy(result, callsign + ws_len); + return 0; // Success +} + + +int unpack_type1(const uint8_t *a77, uint8_t i3, char *message) { + uint32_t n28a, n28b; + uint16_t igrid4; + uint8_t ir; + + // Extract packed fields + // read(c77,1000) n28a,ipa,n28b,ipb,ir,igrid4,i3 + // 1000 format(2(b28,b1),b1,b15,b3) + n28a = (a77[0] << 21); + n28a |= (a77[1] << 13); + n28a |= (a77[2] << 5); + n28a |= (a77[3] >> 3); + n28b = ((a77[3] & 0x07) << 26); + n28b |= (a77[4] << 18); + n28b |= (a77[5] << 10); + n28b |= (a77[6] << 2); + n28b |= (a77[7] >> 6); + ir = ((a77[7] & 0x20) >> 5); + igrid4 = ((a77[7] & 0x1F) << 10); + igrid4 |= (a77[8] << 2); + igrid4 |= (a77[9] >> 6); + + // Unpack both callsigns + char field_1[14]; + char field_2[14]; + if (unpack28(n28a >> 1, field_1) < 0) { + return -1; + } + if (unpack28(n28b >> 1, field_2) < 0) { + return -2; + } + // Fix "CQ_" to "CQ " -> already done in unpack28() + // if (starts_with(field_1, "CQ_")) { + // field_1[2] = ' '; + // } + + // Append first two fields to the result + strcpy(message, field_1); + 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 + + char field_3[5]; + if (igrid4 <= MAXGRID4) { + // Extract 4 symbol grid locator + field_3[4] = '\0'; + + uint16_t n = igrid4; + field_3[3] = '0' + (n % 10); + n /= 10; + field_3[2] = '0' + (n % 10); + n /= 10; + field_3[1] = 'A' + (n % 18); + n /= 18; + field_3[0] = 'A' + (n % 18); + + if (ir != 0) { + // In case of ir=1 add an " R " before grid + strcat(message, " R "); + } + } + else { + // Extract report + int irpt = igrid4 - MAXGRID4; + + // Check special cases first + if (irpt == 1) field_3[0] = '\0'; + else if (irpt == 2) strcpy(field_3, "RRR"); + else if (irpt == 3) strcpy(field_3, "RR73"); + else if (irpt == 4) strcpy(field_3, "73"); + else if (irpt >= 5) { + // Extract signal report as a two digit number with a + or - sign + if (ir == 0) { + int_to_dd(field_3, irpt - 35, 2, true); + } + else { + field_3[0] = 'R'; + int_to_dd(field_3 + 1, irpt - 35, 2, true); + } + } + } + + // Append the last field to the result + if (strlen(field_3) > 0) { + strcat(message, " "); + strcat(message, field_3); + } + + return 0; // Success +} + + +int unpack_text(const uint8_t *a71, char *text) { + // TODO: test + uint8_t b71[9]; + + for (int i = 0; i < 9; ++i) { + b71[i] = a71[i]; + } + + for (int idx = 0; idx < 13; ++idx) { + // Divide the long integer in b71 by 42 + uint16_t rem = 0; + for (int i = 8; i >= 0; --i) { + rem = (rem << 8) | b71[i]; + b71[i] = rem / 42; + rem = rem % 42; + } + text[idx] = charn(rem, 0); + } + + text[13] = '\0'; + return 0; // Success +} + + +int unpack_telemetry(const uint8_t *a71, char *telemetry) { + uint8_t b71[9]; + + // Shift bits in a71 right by 1 + uint8_t carry = 0; + for (int i = 0; i < 9; ++i) { + b71[i] = (carry << 7) | (a71[i] >> 1); + carry = (a71[i] & 0x01); + } + + // Convert b71 to hexadecimal string + for (int i = 0; i < 9; ++i) { + uint8_t nibble1 = (b71[i] >> 4); + uint8_t nibble2 = (b71[i] & 0x0F); + char c1 = (nibble1 > 9) ? (nibble1 - 10 + 'A') : nibble1 + '0'; + char c2 = (nibble2 > 9) ? (nibble2 - 10 + 'A') : nibble2 + '0'; + telemetry[i * 2] = c1; + telemetry[i * 2 + 1] = c2; + } + + telemetry[18] = '\0'; + return 0; +} + + +int unpack77(const uint8_t *a77, char *message) { + uint8_t n3, i3; + + n3 = (a77[9] >> 6) & 0x07; + i3 = (a77[9] >> 3) & 0x07; + + if (i3 == 0 && n3 == 0) { + // 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 == 5) { + // 0.5 0123456789abcdef01 71 71 Telemetry (18 hex) + return unpack_telemetry(a77, message); + } + else if (i3 == 1 || i3 == 2) { + // 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. + + // TODO: implement + // read(c77,1050) n12,n58,iflip,nrpt,icq + // 1050 format(b12,b58,b1,b2,b1) + } + + return 0; +} diff --git a/ft8/unpack_v2.h b/ft8/unpack_v2.h new file mode 100644 index 0000000..2771d4c --- /dev/null +++ b/ft8/unpack_v2.h @@ -0,0 +1,7 @@ +#pragma once + +#include + + +// message should have at least 19 bytes allocated (18 characters + zero terminator) +int unpack77(const uint8_t *a77, char *message); diff --git a/ft8/v1/pack.cpp b/ft8/v1/pack.cpp index 48d495b..f4a10f2 100644 --- a/ft8/v1/pack.cpp +++ b/ft8/v1/pack.cpp @@ -2,7 +2,7 @@ #include #include "pack.h" -#include "text.h" +#include "../text.h" constexpr int32_t NBASE = 37*36*10*27*27*27L; constexpr int32_t NGBASE = 180*180L; diff --git a/ft8/v1/unpack.cpp b/ft8/v1/unpack.cpp index 5529048..dc8f24e 100644 --- a/ft8/v1/unpack.cpp +++ b/ft8/v1/unpack.cpp @@ -1,5 +1,5 @@ #include "unpack.h" -#include "text.h" +#include "../text.h" #include diff --git a/gen_ft8.cpp b/gen_ft8.cpp index 17b073f..3be8621 100644 --- a/gen_ft8.cpp +++ b/gen_ft8.cpp @@ -84,15 +84,17 @@ int main(int argc, char **argv) { printf("\n"); // Third, convert the FSK tones into an audio signal - const int num_samples = (int)(0.5 + FT8_NN / 6.25 * 12000); - const int num_silence = (15 * 12000 - num_samples) / 2; + const int sample_rate = 12000; + const float symbol_rate = 6.25f; + const int num_samples = (int)(0.5f + FT8_NN / symbol_rate * sample_rate); + const int num_silence = (15 * sample_rate - num_samples) / 2; float signal[num_silence + num_samples + num_silence]; for (int i = 0; i < num_silence + num_samples + num_silence; i++) { signal[i] = 0; } - synth_fsk(tones, FT8_NN, 1000, 6.25, 6.25, 12000, signal + num_silence); - save_wav(signal, num_silence + num_samples + num_silence, 12000, wav_path); + synth_fsk(tones, FT8_NN, 1000, symbol_rate, symbol_rate, sample_rate, signal + num_silence); + save_wav(signal, num_silence + num_samples + num_silence, sample_rate, wav_path); return 0; } diff --git a/test.cpp b/test.cpp index 03e99aa..85a5866 100644 --- a/test.cpp +++ b/test.cpp @@ -4,11 +4,11 @@ #include #include "ft8/text.h" -#include "ft8/pack.h" +#include "ft8/v1/pack.h" +#include "ft8/v1/unpack.h" +#include "ft8/v1/encode.h" #include "ft8/pack_v2.h" -#include "ft8/encode.h" #include "ft8/encode_v2.h" -#include "ft8/unpack.h" #include "common/debug.h"