diff --git a/.gitignore b/.gitignore index e80131a..1600647 100644 --- a/.gitignore +++ b/.gitignore @@ -2,8 +2,9 @@ gen_ft8 decode_ft8 test_ft8 libft8.a -wsjtx2/ +wsjtx/ .build/ +.env/ .DS_Store .vscode/ __pycache__/ \ No newline at end of file diff --git a/common/common.h b/common/common.h index fc2a2f9..241115d 100644 --- a/common/common.h +++ b/common/common.h @@ -1,3 +1,5 @@ #ifndef M_PI #define M_PI 3.14159265358979323846 #endif + +#define TWO_PI 6.28318530717958647692 \ No newline at end of file diff --git a/common/monitor.c b/common/monitor.c index 6bd228a..734936f 100644 --- a/common/monitor.c +++ b/common/monitor.c @@ -12,27 +12,27 @@ static float hann_i(int i, int N) return x * x; } -// static float hamming_i(int i, int N) -// { -// const float a0 = (float)25 / 46; -// const float a1 = 1 - a0; +static float hamming_i(int i, int N) +{ + const float a0 = (float)25 / 46; + const float a1 = 1 - a0; -// float x1 = cosf(2 * (float)M_PI * i / N); -// return a0 - a1 * x1; -// } + float x1 = cosf((float)TWO_PI * i / N); + return a0 - a1 * x1; +} -// static float blackman_i(int i, int N) -// { -// const float alpha = 0.16f; // or 2860/18608 -// const float a0 = (1 - alpha) / 2; -// const float a1 = 1.0f / 2; -// const float a2 = alpha / 2; +static float blackman_i(int i, int N) +{ + const float alpha = 0.16f; // or 2860/18608 + const float a0 = (1 - alpha) / 2; + const float a1 = 1.0f / 2; + const float a2 = alpha / 2; -// float x1 = cosf(2 * (float)M_PI * i / N); -// float x2 = 2 * x1 * x1 - 1; // Use double angle formula + float x1 = cosf((float)TWO_PI * i / N); + float x2 = 2 * x1 * x1 - 1; // Use double angle formula -// return a0 - a1 * x1 + a2 * x2; -// } + return a0 - a1 * x1 + a2 * x2; +} static void waterfall_init(ftx_waterfall_t* me, int max_blocks, int num_bins, int time_osr, int freq_osr) { @@ -86,7 +86,7 @@ void monitor_init(monitor_t* me, const monitor_config_t* cfg) LOG(LOG_DEBUG, "FFT work area = %zu\n", fft_work_size); #ifdef WATERFALL_USE_PHASE - me->nifft = 64; // Gives 200 Hz sample rate for FT8 (160ms symbol period) + me->nifft = 64; // Gives 100 Hz sample rate for FT8 (160ms symbol period) size_t ifft_work_size = 0; kiss_fft_alloc(me->nifft, 1, 0, &ifft_work_size); @@ -192,20 +192,14 @@ void monitor_process(monitor_t* me, const float* frame) } #ifdef WATERFALL_USE_PHASE -void monitor_resynth(const monitor_t* me, const candidate_t* candidate, float* signal) +void monitor_resynth(const monitor_t* me, const ftx_candidate_t* cand, float* signal) { const int num_ifft = me->nifft; const int num_shift = num_ifft / 2; - const int taper_width = 4; - const int num_tones = 8; + const int taper_width = 1; - // Starting offset is 3 subblocks due to analysis buffer loading - int offset = 1; // candidate->time_offset; - offset = (offset * me->wf.time_osr) + 1; // + candidate->time_sub; - offset = (offset * me->wf.freq_osr); // + candidate->freq_sub; - offset = (offset * me->wf.num_bins); // + candidate->freq_offset; - - WF_ELEM_T* el = me->wf.mag + offset; + // Starting offset is 1 block and 1 subblock due to how data is shifted to analysis buffer + WF_ELEM_T* el = wfall_get_element(&me->wf, 0, 0, 0, cand->freq_sub); // DFT frequency data - initialize to zero kiss_fft_cpx freqdata[num_ifft]; @@ -218,20 +212,20 @@ void monitor_resynth(const monitor_t* me, const candidate_t* candidate, float* s int pos = 0; for (int num_block = 1; num_block < me->wf.num_blocks; ++num_block) { - // Extract frequency data around the selected candidate only - for (int i = candidate->freq_offset - taper_width - 1; i < candidate->freq_offset + 8 + taper_width - 1; ++i) + // Extract frequency data around the selected cand only + for (int i = cand->freq_offset - taper_width - 1; i < cand->freq_offset + 8 + taper_width - 1; ++i) { if ((i >= 0) && (i < me->wf.num_bins)) { - int tgt_bin = (me->wf.freq_osr * (i - candidate->freq_offset) + num_ifft) % num_ifft; + int tgt_bin = (me->wf.freq_osr * (i - cand->freq_offset) + num_ifft) % num_ifft; float weight = 1.0f; - if (i < candidate->freq_offset) + if (i < cand->freq_offset) { - weight = ((i - candidate->freq_offset) + taper_width) / (float)taper_width; + weight = ((i - cand->freq_offset) + taper_width) / (float)taper_width; } - else if (i > candidate->freq_offset + 7) + else if (i > cand->freq_offset + 7) { - weight = ((candidate->freq_offset + 7 - i) + taper_width) / (float)taper_width; + weight = ((cand->freq_offset + 7 - i) + taper_width) / (float)taper_width; } // Convert (dB magnitude, phase) to (real, imaginary) diff --git a/common/monitor.h b/common/monitor.h index dd70d8d..32d7234 100644 --- a/common/monitor.h +++ b/common/monitor.h @@ -52,7 +52,7 @@ void monitor_process(monitor_t* me, const float* frame); void monitor_free(monitor_t* me); #ifdef WATERFALL_USE_PHASE -void monitor_resynth(const monitor_t* me, const candidate_t* candidate, float* signal); +void monitor_resynth(const monitor_t* me, const ftx_candidate_t* candidate, float* signal); #endif #ifdef __cplusplus diff --git a/demo/decode_ft8.c b/demo/decode_ft8.c index 1f45967..2ea504e 100644 --- a/demo/decode_ft8.c +++ b/demo/decode_ft8.c @@ -22,6 +22,7 @@ const int kMax_candidates = 140; const int kLDPC_iterations = 25; const int kMax_decoded_messages = 50; +const int kMax_callsign_age = 10; const int kFreq_osr = 2; // Frequency oversampling rate (bin subdivision) const int kTime_osr = 2; // Time oversampling rate (symbol subdivision) @@ -152,22 +153,33 @@ void decode(const monitor_t* mon, struct tm* tm_slot_start) float freq_hz = (mon->min_bin + cand->freq_offset + (float)cand->freq_sub / wf->freq_osr) / mon->symbol_period; float time_sec = (cand->time_offset + (float)cand->time_sub / wf->time_osr) * mon->symbol_period; + float log174[FTX_LDPC_N]; + #ifdef WATERFALL_USE_PHASE - // int resynth_len = 12000 * 16; - // float resynth_signal[resynth_len]; - // for (int pos = 0; pos < resynth_len; ++pos) - // { - // resynth_signal[pos] = 0; - // } - // monitor_resynth(mon, cand, resynth_signal); - // char resynth_path[80]; - // sprintf(resynth_path, "resynth_%04f_%02.1f.wav", freq_hz, time_sec); - // save_wav(resynth_signal, resynth_len, 12000, resynth_path); + float slot_period = ((wf->protocol == FTX_PROTOCOL_FT8) ? FT8_SLOT_TIME : FT4_SLOT_TIME); + int resamp_len = wf->num_blocks * FT8_DS_SYM_LEN; + float resamp_wave[resamp_len]; + for (int pos = 0; pos < resamp_len; ++pos) + { + resamp_wave[pos] = 0; + } + monitor_resynth(mon, cand, resamp_wave); + + int offset_crude = cand->time_offset * FT8_DS_SYM_LEN + cand->time_sub * FT8_DS_SYM_LEN / 2; + int offset_fine; + float df_fine; + find_sync_fine(resamp_wave, resamp_len, FT8_DS_RATE, offset_crude, &offset_fine, &df_fine); + + float dfphi = (float)TWO_PI * df_fine / FT8_DS_RATE; + extract_likelihood_fine(resamp_wave, resamp_len, offset_crude + offset_fine, dfphi, log174); +#else + ftx_extract_likelihood(wf, cand, log174); #endif ftx_message_t message; ftx_decode_status_t status; - if (!ftx_decode_candidate(wf, cand, kLDPC_iterations, &message, &status)) + // if (!ftx_decode_candidate(wf, cand, kLDPC_iterations, &message, &status)) + if (!ftx_decode_candidate(log174, wf->protocol, kLDPC_iterations, &message, &status)) { if (status.ldpc_errors > 0) { @@ -223,10 +235,17 @@ void decode(const monitor_t* mon, struct tm* tm_slot_start) printf("%02d%02d%02d %+05.1f %+4.2f %4.0f ~ %s\n", tm_slot_start->tm_hour, tm_slot_start->tm_min, tm_slot_start->tm_sec, snr, time_sec, freq_hz, text); + +#ifdef WATERFALL_USE_PHASE + LOG(LOG_INFO, "offset_crude = %d, offset_fine = %d, df_fine = %f\n", offset_crude, offset_fine, df_fine); + char resamp_path[80]; + snprintf(resamp_path, 80, "test/resamp/resamp_%04.0f_%04.0f.wav", freq_hz, 1000 * time_sec); + save_wav(resamp_wave, resamp_len, FT8_DS_RATE, resamp_path); +#endif } } LOG(LOG_INFO, "Decoded %d messages, callsign hashtable size %d\n", num_decoded, callsign_hashtable_size); - hashtable_cleanup(10); + hashtable_cleanup(kMax_callsign_age); } int main(int argc, char** argv) diff --git a/demo/gen_ft8.c b/demo/gen_ft8.c index 0f56068..4976868 100644 --- a/demo/gen_ft8.c +++ b/demo/gen_ft8.c @@ -4,8 +4,8 @@ #include #include -#include "common/common.h" -#include "common/wave.h" +#include +#include #include "ft8/message.h" #include "ft8/encode.h" #include "ft8/constants.h" @@ -55,13 +55,13 @@ void synth_gfsk(const uint8_t* symbols, int n_sym, float f0, float symbol_bt, fl LOG(LOG_DEBUG, "n_spsym = %d\n", n_spsym); // Compute the smoothed frequency waveform. // Length = (nsym+2)*n_spsym samples, first and last symbols extended - float dphi_peak = 2 * M_PI * hmod / n_spsym; + float dphi_peak = (float)TWO_PI * hmod / n_spsym; float dphi[n_wave + 2 * n_spsym]; // Shift frequency up by f0 for (int i = 0; i < n_wave + 2 * n_spsym; ++i) { - dphi[i] = 2 * M_PI * f0 / signal_rate; + dphi[i] = (float)TWO_PI * f0 / signal_rate; } float pulse[3 * n_spsym]; @@ -88,14 +88,14 @@ void synth_gfsk(const uint8_t* symbols, int n_sym, float f0, float symbol_bt, fl for (int k = 0; k < n_wave; ++k) { // Don't include dummy symbols signal[k] = sinf(phi); - phi = fmodf(phi + dphi[k + n_spsym], 2 * M_PI); + phi = fmodf(phi + dphi[k + n_spsym], (float)TWO_PI); } // Apply envelope shaping to the first and last symbols int n_ramp = n_spsym / 8; for (int i = 0; i < n_ramp; ++i) { - float env = (1 - cosf(2 * M_PI * i / (2 * n_ramp))) / 2; + float env = (1 - cosf((float)TWO_PI * i / (2 * n_ramp))) / 2; signal[i] *= env; signal[n_wave - 1 - i] *= env; } diff --git a/ft8/constants.h b/ft8/constants.h index eb50fca..6e82e31 100644 --- a/ft8/constants.h +++ b/ft8/constants.h @@ -10,9 +10,11 @@ extern "C" #define FT8_SYMBOL_PERIOD (0.160f) ///< FT8 symbol duration, defines tone deviation in Hz and symbol rate #define FT8_SLOT_TIME (15.0f) ///< FT8 slot period +#define FT8_NUM_TONES (8) #define FT4_SYMBOL_PERIOD (0.048f) ///< FT4 symbol duration, defines tone deviation in Hz and symbol rate #define FT4_SLOT_TIME (7.5f) ///< FT4 slot period +#define FT4_NUM_TONES (4) // Define FT8 symbol counts // FT8 message structure: diff --git a/ft8/decode.c b/ft8/decode.c index efb1718..92ab1f2 100644 --- a/ft8/decode.c +++ b/ft8/decode.c @@ -5,13 +5,14 @@ #include #include +#include // #define LOG_LEVEL LOG_DEBUG // #include "debug.h" // Lookup table for y = 10*log10(1 + 10^(x/10)), where -// y - increase in signal level dB when adding a weaker independent signal -// x - specific relative strength of the weaker signal in dB +// y - increase in signal power dB when adding a weaker independent signal +// x - relative power of the weaker signal in dB // Table index corresponds to x in dB (index 0: 0 dB, index 1: -1 dB etc) static const float db_power_sum[40] = { 3.01029995663981f, 2.53901891043867f, 2.1244260279434f, 1.76434862436485f, 1.45540463109294f, @@ -51,11 +52,7 @@ static void ft8_decode_multi_symbols(const WF_ELEM_T* wf, int num_bins, int n_sy static const WF_ELEM_T* get_cand_mag(const ftx_waterfall_t* wf, const ftx_candidate_t* candidate) { - int offset = candidate->time_offset; - offset = (offset * wf->time_osr) + candidate->time_sub; - offset = (offset * wf->freq_osr) + candidate->freq_sub; - offset = (offset * wf->num_bins) + candidate->freq_offset; - return wf->mag + offset; + return wfall_get_element(wf, candidate->time_offset, candidate->time_sub, candidate->freq_offset, candidate->freq_sub); } static int ft8_sync_score(const ftx_waterfall_t* wf, const ftx_candidate_t* candidate) @@ -84,10 +81,6 @@ static int ft8_sync_score(const ftx_waterfall_t* wf, const ftx_candidate_t* cand // Weighted difference between the expected and all other symbols // Does not work as well as the alternative score below - // score += 8 * p8[kFT8_Costas_pattern[k]] - - // p8[0] - p8[1] - p8[2] - p8[3] - - // p8[4] - p8[5] - p8[6] - p8[7]; - // ++num_average; // Check only the neighbors of the expected symbol frequency- and time-wise int sm = kFT8_Costas_pattern[k]; // Index of the expected bin @@ -150,9 +143,6 @@ static int ft4_sync_score(const ftx_waterfall_t* wf, const ftx_candidate_t* cand int sm = kFT4_Costas_pattern[m][k]; // Index of the expected bin - // score += (4 * p4[sm]) - p4[0] - p4[1] - p4[2] - p4[3]; - // num_average += 4; - // Check only the neighbors of the expected symbol frequency- and time-wise if (sm > 0) { @@ -250,6 +240,20 @@ int ftx_find_candidates(const ftx_waterfall_t* wf, int num_candidates, ftx_candi return heap_size; } +void ftx_extract_likelihood(const ftx_waterfall_t* wf, const ftx_candidate_t* cand, float* log174) +{ + if (wf->protocol == FTX_PROTOCOL_FT4) + { + ft4_extract_likelihood(wf, cand, log174); + } + else + { + ft8_extract_likelihood(wf, cand, log174); + } + + ftx_normalize_logl(log174); +} + static void ft4_extract_likelihood(const ftx_waterfall_t* wf, const ftx_candidate_t* cand, float* log174) { const WF_ELEM_T* mag = get_cand_mag(wf, cand); // Pointer to 4 magnitude bins of the first symbol @@ -324,20 +328,8 @@ static void ftx_normalize_logl(float* log174) } } -bool ftx_decode_candidate(const ftx_waterfall_t* wf, const ftx_candidate_t* cand, int max_iterations, ftx_message_t* message, ftx_decode_status_t* status) +bool ftx_decode_candidate(const float* log174, ftx_protocol_t protocol, int max_iterations, ftx_message_t* message, ftx_decode_status_t* status) { - float log174[FTX_LDPC_N]; // message bits encoded as likelihood - if (wf->protocol == FTX_PROTOCOL_FT4) - { - ft4_extract_likelihood(wf, cand, log174); - } - else - { - ft8_extract_likelihood(wf, cand, log174); - } - - ftx_normalize_logl(log174); - uint8_t plain174[FTX_LDPC_N]; // message bits (0/1) bp_decode(log174, max_iterations, plain174, &status->ldpc_errors); // ldpc_decode(log174, max_iterations, plain174, &status->ldpc_errors); @@ -366,7 +358,7 @@ bool ftx_decode_candidate(const ftx_waterfall_t* wf, const ftx_candidate_t* cand // Reuse CRC value as a hash for the message (TODO: 14 bits only, should perhaps use full 16 or 32 bits?) message->hash = status->crc_calculated; - if (wf->protocol == FTX_PROTOCOL_FT4) + if (protocol == FTX_PROTOCOL_FT4) { // '[..] for FT4 only, in order to avoid transmitting a long string of zeros when sending CQ messages, // the assembled 77-bit message is bitwise exclusive-OR’ed with [a] pseudorandom sequence before computing the CRC and FEC parity bits' diff --git a/ft8/decode.h b/ft8/decode.h index 43cd376..11bcf79 100644 --- a/ft8/decode.h +++ b/ft8/decode.h @@ -5,6 +5,7 @@ #include #include "constants.h" +#include "waterfall.h" #include "message.h" #ifdef __cplusplus @@ -12,53 +13,6 @@ extern "C" { #endif -typedef struct -{ - float mag; - float phase; -} waterfall_cpx_t; - -// #define WATERFALL_USE_PHASE - -#ifdef WATERFALL_USE_PHASE -#define WF_ELEM_T waterfall_cpx_t -#define WF_ELEM_MAG(x) ((x).mag) -#define WF_ELEM_MAG_INT(x) (int)(2 * ((x).mag + 120.0f)) -#else -#define WF_ELEM_T uint8_t -#define WF_ELEM_MAG(x) ((float)(x)*0.5f - 120.0f) -#define WF_ELEM_MAG_INT(x) (int)(x) -#endif - -/// Input structure to ftx_find_sync() function. This structure describes stored waterfall data over the whole message slot. -/// Fields time_osr and freq_osr specify additional oversampling rate for time and frequency resolution. -/// If time_osr=1, FFT magnitude data is collected once for every symbol transmitted, i.e. every 1/6.25 = 0.16 seconds. -/// Values time_osr > 1 mean each symbol is further subdivided in time. -/// If freq_osr=1, each bin in the FFT magnitude data corresponds to 6.25 Hz, which is the tone spacing. -/// Values freq_osr > 1 mean the tone spacing is further subdivided by FFT analysis. -typedef struct -{ - int max_blocks; ///< number of blocks (symbols) allocated in the mag array - int num_blocks; ///< number of blocks (symbols) stored in the mag array - int num_bins; ///< number of FFT bins in terms of 6.25 Hz - int time_osr; ///< number of time subdivisions - int freq_osr; ///< number of frequency subdivisions - WF_ELEM_T* mag; ///< FFT magnitudes stored as uint8_t[blocks][time_osr][freq_osr][num_bins] - int block_stride; ///< Helper value = time_osr * freq_osr * num_bins - ftx_protocol_t protocol; ///< Indicate if using FT4 or FT8 -} ftx_waterfall_t; - -/// Output structure of ftx_find_sync() and input structure of ftx_decode(). -/// Holds the position of potential start of a message in time and frequency. -typedef struct -{ - int16_t score; ///< Candidate score (non-negative number; higher score means higher likelihood) - int16_t time_offset; ///< Index of the time block - int16_t freq_offset; ///< Index of the frequency bin - uint8_t time_sub; ///< Index of the time subdivision used - uint8_t freq_sub; ///< Index of the frequency subdivision used -} ftx_candidate_t; - /// Structure that contains the status of various steps during decoding of a message typedef struct { @@ -72,22 +26,25 @@ typedef struct /// 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). -/// @param[in] power Waterfall data collected during message slot +/// @param[in] wf Waterfall data collected during message slot /// @param[in] sync_pattern Synchronization pattern /// @param[in] num_candidates Number of maximum candidates (size of heap array) /// @param[in,out] heap Array of ftx_candidate_t type entries (with num_candidates allocated entries) /// @param[in] min_score Minimal score allowed for pruning unlikely candidates (can be zero for no effect) /// @return Number of candidates filled in the heap -int ftx_find_candidates(const ftx_waterfall_t* power, int num_candidates, ftx_candidate_t heap[], int min_score); +int ftx_find_candidates(const ftx_waterfall_t* wf, int num_candidates, ftx_candidate_t heap[], int min_score); + +void ftx_extract_likelihood(const ftx_waterfall_t* wf, const ftx_candidate_t* cand, float* log174); /// Attempt to decode a message candidate. Extracts the bit probabilities, runs LDPC decoder, checks CRC and unpacks the message in plain text. -/// @param[in] power Waterfall data collected during message slot +/// @param[in] wf Waterfall data collected during message slot /// @param[in] cand Candidate to decode /// @param[in] max_iterations Maximum allowed LDPC iterations (lower number means faster decode, but less precise) /// @param[out] message ftx_message_t structure that will receive the decoded message /// @param[out] status ftx_decode_status_t structure that will be filled with the status of various decoding steps /// @return True if the decoding was successful, false otherwise (check status for details) -bool ftx_decode_candidate(const ftx_waterfall_t* power, const ftx_candidate_t* cand, int max_iterations, ftx_message_t* message, ftx_decode_status_t* status); +// bool ftx_decode_candidate(const ftx_waterfall_t* wf, const ftx_candidate_t* cand, int max_iterations, ftx_message_t* message, ftx_decode_status_t* status); +bool ftx_decode_candidate(const float* log174, ftx_protocol_t protocol, int max_iterations, ftx_message_t* message, ftx_decode_status_t* status); #ifdef __cplusplus } diff --git a/ft8/ldpc.c b/ft8/ldpc.c index 539a1df..e8c0d48 100644 --- a/ft8/ldpc.c +++ b/ft8/ldpc.c @@ -25,7 +25,7 @@ static float fast_atanh(float x); // 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) +void ldpc_decode(const float codeword[], int max_iters, uint8_t plain[], int* ok) { float m[FTX_LDPC_M][FTX_LDPC_N]; // ~60 kB float e[FTX_LDPC_M][FTX_LDPC_N]; // ~60 kB @@ -127,7 +127,7 @@ static int ldpc_check(uint8_t codeword[]) return errors; } -void bp_decode(float codeword[], int max_iters, uint8_t plain[], int* ok) +void bp_decode(const float codeword[], int max_iters, uint8_t plain[], int* ok) { float tov[FTX_LDPC_N][3]; float toc[FTX_LDPC_M][7]; diff --git a/ft8/ldpc.h b/ft8/ldpc.h index 02f4a99..6f25823 100644 --- a/ft8/ldpc.h +++ b/ft8/ldpc.h @@ -12,9 +12,9 @@ extern "C" // 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, uint8_t plain[], int* ok); +void ldpc_decode(const float codeword[], int max_iters, uint8_t plain[], int* ok); -void bp_decode(float codeword[], int max_iters, uint8_t plain[], int* ok); +void bp_decode(const float codeword[], int max_iters, uint8_t plain[], int* ok); #ifdef __cplusplus } diff --git a/ft8/waterfall.c b/ft8/waterfall.c new file mode 100644 index 0000000..4aa4f71 --- /dev/null +++ b/ft8/waterfall.c @@ -0,0 +1,180 @@ +#include "waterfall.h" + +#include +#include +#include + +#define LOG_LEVEL LOG_INFO +#include "debug.h" + +WF_ELEM_T* wfall_get_element(const ftx_waterfall_t* wf, + int16_t time_offset, uint8_t time_sub, int16_t freq_offset, uint8_t freq_sub) +{ + int offset = time_offset; + offset = (offset * wf->time_osr) + time_sub; + offset = (offset * wf->freq_osr) + freq_sub; + offset = (offset * wf->num_bins) + freq_offset; + return wf->mag + offset; +} + +static float cabs2f(complex float z) +{ + return crealf(z) * crealf(z) + cimagf(z) * cimagf(z); +} + +static float score_sync_fine(const float* wave, int num_samples, int i0, float dfphi) +{ + float sync = 0; + int num_sync = 0; + + // Sum over 7 Costas frequencies ... + for (int i = 0; i < 7; ++i) + { + complex float csync2[FT8_DS_SYM_LEN]; + float dphi = dfphi + (float)TWO_PI * kFT8_Costas_pattern[i] / FT8_DS_SYM_LEN; + float phi = 0; + for (int k = 0; k < FT8_DS_SYM_LEN; ++k) + { + csync2[k] = cosf(phi) - I * sinf(phi); + phi = fmodf(phi + dphi, (float)TWO_PI); + } + + // ... and three Costas arrays + for (int j = 0; j <= 72; j += 36) + { + int i1 = i0 + (i + j) * FT8_DS_SYM_LEN; + + complex float z = 0; + for (int k = 0; k < FT8_DS_SYM_LEN; ++k) + { + if (i1 + k >= 0 && i1 + k < num_samples) + z += wave[i1 + k] * csync2[k]; + } + sync += cabs2f(z); + ++num_sync; + } + } + return sync; +} + +void find_sync_fine(const float* wave, int num_samples, float sample_rate, int offset_crude, int *offset_fine, float *df_fine) +{ + float score_max = 0; + int i0_max = -1; + for (int i0 = 0; i0 < FT8_DS_SYM_LEN; ++i0) + { + float score = score_sync_fine(wave, num_samples, offset_crude + i0, 0.0f); + if ((i0_max < 0) || (score > score_max)) + { + i0_max = i0; + score_max = score; + } + } + LOG(LOG_DEBUG, "i0_max = %d, score_max = %f\n", i0_max, score_max); + + float df_max = 0; + for (float df = -3.2f; df <= 3.2f; df += 0.1f) + { + float dfphi = (float)TWO_PI * df / sample_rate; + float score = score_sync_fine(wave, num_samples, offset_crude + i0_max, dfphi); + if (score > score_max) + { + df_max = df; + score_max = score; + } + } + LOG(LOG_DEBUG, "df_max = %+.1f, score_max = %f\n", df_max, score_max); + + float dfphi = (float)TWO_PI * df_max / sample_rate; + for (int i0 = 0; i0 < FT8_DS_SYM_LEN; ++i0) + { + float score = score_sync_fine(wave, num_samples, offset_crude + i0, dfphi); + if (score > score_max) + { + i0_max = i0; + score_max = score; + } + } + LOG(LOG_DEBUG, "i0_max = %d, score_max = %f\n", i0_max, score_max); + + *offset_fine = i0_max; + *df_fine = df_max; +} + +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 ftx_normalize_logl(float* log174) +{ + // Compute the variance of log174 + float sum = 0; + float sum2 = 0; + for (int i = 0; i < FTX_LDPC_N; ++i) + { + sum += log174[i]; + sum2 += log174[i] * log174[i]; + } + float inv_n = 1.0f / FTX_LDPC_N; + float variance = (sum2 - (sum * sum * inv_n)) * inv_n; + + // Normalize log174 distribution and scale it with experimentally found coefficient + float norm_factor = sqrtf(16.0f / variance); + for (int i = 0; i < FTX_LDPC_N; ++i) + { + log174[i] *= norm_factor; + } +} + +void extract_likelihood_fine(const float* wave, int num_samples, int start_pos, float dfphi, float* log174) +{ + complex float dft[FT8_NUM_TONES][FT8_DS_SYM_LEN]; + + for (int i = 0; i < FT8_NUM_TONES; ++i) + { + float dphi = dfphi + (float)TWO_PI * kFT8_Gray_map[i] / FT8_DS_SYM_LEN; + float phi = 0; + for (int k = 0; k < FT8_DS_SYM_LEN; ++k) + { + dft[i][k] = cosf(phi) - I * sinf(phi); + phi = fmodf(phi + dphi, (float)TWO_PI); + } + } + + for (int k = 0; k < FT8_ND; ++k) + { + // Skip either 7 or 14 sync symbols + // TODO: replace magic numbers with constants + int sym_idx = k + ((k < 29) ? 7 : 14); + int bit_idx = 3 * k; + + float s2[FT8_NUM_TONES]; + int wave_pos = start_pos + sym_idx * FT8_DS_SYM_LEN; + for (int j = 0; j < FT8_NUM_TONES; ++j) + { + complex float sum = 0; + for (int k = 0; k < FT8_DS_SYM_LEN; ++k) + { + if ((wave_pos + k >= 0) && (wave_pos + k < num_samples)) + { + sum += wave[wave_pos + k] * dft[j][k]; + } + } + float mag2 = cabs2f(sum); + s2[j] = 10.0f * log10f(1E-12f + mag2); + } + + 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]); + // LOG(LOG_INFO, "log174[%d] @%d = %f, %f, %f\n", bit_idx, sym_idx, log174[bit_idx + 0], log174[bit_idx + 1], log174[bit_idx + 2]); + } + + ftx_normalize_logl(log174); +} diff --git a/ft8/waterfall.h b/ft8/waterfall.h new file mode 100644 index 0000000..a89064c --- /dev/null +++ b/ft8/waterfall.h @@ -0,0 +1,72 @@ +#ifndef _INCLUDE_WATERFALL_H_ +#define _INCLUDE_WATERFALL_H_ + +#include "constants.h" + +#ifdef __cplusplus +extern "C" +{ +#endif + +typedef struct +{ + float mag; + float phase; +} waterfall_cpx_t; + +#define WATERFALL_USE_PHASE + +#ifdef WATERFALL_USE_PHASE +#define WF_ELEM_T waterfall_cpx_t +#define WF_ELEM_MAG(x) ((x).mag) +#define WF_ELEM_MAG_INT(x) (int)(2 * ((x).mag + 120.0f)) +#else +#define WF_ELEM_T uint8_t +#define WF_ELEM_MAG(x) ((float)(x) * 0.5f - 120.0f) +#define WF_ELEM_MAG_INT(x) (int)(x) +#endif + +/// Input structure to ftx_find_sync() function. This structure describes stored waterfall data over the whole message slot. +/// Fields time_osr and freq_osr specify additional oversampling rate for time and frequency resolution. +/// If time_osr=1, FFT magnitude data is collected once for every symbol transmitted, i.e. every 1/6.25 = 0.16 seconds. +/// Values time_osr > 1 mean each symbol is further subdivided in time. +/// If freq_osr=1, each bin in the FFT magnitude data corresponds to 6.25 Hz, which is the tone spacing. +/// Values freq_osr > 1 mean the tone spacing is further subdivided by FFT analysis. +typedef struct +{ + int max_blocks; ///< number of blocks (symbols) allocated in the mag array + int num_blocks; ///< number of blocks (symbols) stored in the mag array + int num_bins; ///< number of FFT bins in terms of 6.25 Hz + int time_osr; ///< number of time subdivisions + int freq_osr; ///< number of frequency subdivisions + WF_ELEM_T* mag; ///< FFT magnitudes stored as uint8_t[blocks][time_osr][freq_osr][num_bins] + int block_stride; ///< Helper value = time_osr * freq_osr * num_bins + ftx_protocol_t protocol; ///< Indicate if using FT4 or FT8 +} ftx_waterfall_t; + +/// Output structure of ftx_find_sync() and input structure of ftx_decode(). +/// Holds the position of potential start of a message in time and frequency. +typedef struct +{ + int16_t score; ///< Candidate score (non-negative number; higher score means higher likelihood) + int16_t time_offset; ///< Index of the time block + int16_t freq_offset; ///< Index of the frequency bin + uint8_t time_sub; ///< Index of the time subdivision used + uint8_t freq_sub; ///< Index of the frequency subdivision used +} ftx_candidate_t; + + +#define FT8_DS_SYM_LEN (32) +#define FT8_DS_RATE (200) + +WF_ELEM_T* wfall_get_element(const ftx_waterfall_t* wf, + int16_t time_offset, uint8_t time_sub, int16_t freq_offset, uint8_t freq_sub); + +void find_sync_fine(const float* wave, int num_samples, float sample_rate, int offset_crude, int *offset_fine, float *df_fine); +void extract_likelihood_fine(const float* wave, int num_samples, int start_pos, float dfphi, float* log174); + +#ifdef __cplusplus +} +#endif + +#endif // _INCLUDE_WATERFALL_H_ \ No newline at end of file