diff --git a/common/monitor.c b/common/monitor.c index b2c596c..0904c2e 100644 --- a/common/monitor.c +++ b/common/monitor.c @@ -42,7 +42,7 @@ static void waterfall_init(waterfall_t* me, int max_blocks, int num_bins, int ti me->time_osr = time_osr; me->freq_osr = freq_osr; me->block_stride = (time_osr * freq_osr * num_bins); - me->mag = (uint8_t*)malloc(mag_size); + me->mag = (WF_ELEM_T*)malloc(mag_size); LOG(LOG_DEBUG, "Waterfall size = %zu\n", mag_size); } @@ -66,23 +66,35 @@ void monitor_init(monitor_t* me, const monitor_config_t* cfg) for (int i = 0; i < me->nfft; ++i) { // window[i] = 1; - me->window[i] = hann_i(i, me->nfft); + me->window[i] = me->fft_norm * hann_i(i, me->nfft); // me->window[i] = blackman_i(i, me->nfft); // me->window[i] = hamming_i(i, me->nfft); // me->window[i] = (i < len_window) ? hann_i(i, len_window) : 0; } me->last_frame = (float*)calloc(me->nfft, sizeof(me->last_frame[0])); - size_t fft_work_size; - kiss_fftr_alloc(me->nfft, 0, 0, &fft_work_size); - LOG(LOG_INFO, "Block size = %d\n", me->block_size); LOG(LOG_INFO, "Subblock size = %d\n", me->subblock_size); + + size_t fft_work_size = 0; + kiss_fftr_alloc(me->nfft, 0, 0, &fft_work_size); + me->fft_work = malloc(fft_work_size); + me->fft_cfg = kiss_fftr_alloc(me->nfft, 0, me->fft_work, &fft_work_size); + LOG(LOG_INFO, "N_FFT = %d\n", me->nfft); LOG(LOG_DEBUG, "FFT work area = %zu\n", fft_work_size); - me->fft_work = malloc(fft_work_size); - me->fft_cfg = kiss_fftr_alloc(me->nfft, 0, me->fft_work, &fft_work_size); +#ifdef WATERFALL_USE_PHASE + me->nifft = 64; // Gives 200 Hz sample rate for FT8 (160ms symbol period) + + size_t ifft_work_size = 0; + kiss_fft_alloc(me->nifft, 1, 0, &ifft_work_size); + me->ifft_work = malloc(ifft_work_size); + me->ifft_cfg = kiss_fft_alloc(me->nifft, 1, me->ifft_work, &ifft_work_size); + + LOG(LOG_INFO, "N_iFFT = %d\n", me->nifft); + LOG(LOG_DEBUG, "iFFT work area = %zu\n", ifft_work_size); +#endif // Allocate enough blocks to fit the entire FT8/FT4 slot in memory const int max_blocks = (int)(slot_time / symbol_period); @@ -140,15 +152,14 @@ void monitor_process(monitor_t* me, const float* frame) ++frame_pos; } - // Compute windowed analysis frame + // Do DFT of windowed analysis frame for (int pos = 0; pos < me->nfft; ++pos) { - timedata[pos] = me->fft_norm * me->window[pos] * me->last_frame[pos]; + timedata[pos] = me->window[pos] * me->last_frame[pos]; } - kiss_fftr(me->fft_cfg, timedata, freqdata); - // Loop over two possible frequency bin offsets (for averaging) + // Loop over possible frequency OSR offsets for (int freq_sub = 0; freq_sub < me->wf.freq_osr; ++freq_sub) { for (int bin = me->min_bin; bin < me->max_bin; ++bin) @@ -156,11 +167,18 @@ void monitor_process(monitor_t* me, const float* frame) int src_bin = (bin * me->wf.freq_osr) + freq_sub; float mag2 = (freqdata[src_bin].i * freqdata[src_bin].i) + (freqdata[src_bin].r * freqdata[src_bin].r); float db = 10.0f * log10f(1E-12f + mag2); + +#ifdef WATERFALL_USE_PHASE + // Save the magnitude in dB and phase in radians + float phase = atan2f(freqdata[src_bin].i, freqdata[src_bin].r); + me->wf.mag[offset].mag = db; + me->wf.mag[offset].phase = phase; +#else // Scale decibels to unsigned 8-bit range and clamp the value // Range 0-240 covers -120..0 dB in 0.5 dB steps int scaled = (int)(2 * db + 240); - me->wf.mag[offset] = (scaled < 0) ? 0 : ((scaled > 255) ? 255 : scaled); +#endif ++offset; if (db > me->max_mag) @@ -171,3 +189,74 @@ void monitor_process(monitor_t* me, const float* frame) ++me->wf.num_blocks; } + +#ifdef WATERFALL_USE_PHASE +void monitor_resynth(const monitor_t* me, const candidate_t* candidate, 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; + + // 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; + + // DFT frequency data - initialize to zero + kiss_fft_cpx freqdata[num_ifft]; + for (int i = 0; i < num_ifft; ++i) + { + freqdata[i].r = 0; + freqdata[i].i = 0; + } + + 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) + { + if ((i >= 0) && (i < me->wf.num_bins)) + { + int tgt_bin = (me->wf.freq_osr * (i - candidate->freq_offset) + num_ifft) % num_ifft; + float weight = 1.0f; + if (i < candidate->freq_offset) + { + weight = ((i - candidate->freq_offset) + taper_width) / (float)taper_width; + } + else if (i > candidate->freq_offset + 7) + { + weight = ((candidate->freq_offset + 7 - i) + taper_width) / (float)taper_width; + } + + // Convert (dB magnitude, phase) to (real, imaginary) + float mag = powf(10.0f, el[i].mag / 20) / 2 * weight; + freqdata[tgt_bin].r = mag * cosf(el[i].phase); + freqdata[tgt_bin].i = mag * sinf(el[i].phase); + + int i2 = i + me->wf.num_bins; + tgt_bin = (tgt_bin + 1) % num_ifft; + float mag2 = powf(10.0f, el[i2].mag / 20) / 2 * weight; + freqdata[tgt_bin].r = mag2 * cosf(el[i2].phase); + freqdata[tgt_bin].i = mag2 * sinf(el[i2].phase); + } + } + + // Compute inverse DFT and overlap-add the waveform + kiss_fft_cpx timedata[num_ifft]; + kiss_fft(me->ifft_cfg, freqdata, timedata); + for (int i = 0; i < num_ifft; ++i) + { + signal[pos + i] += timedata[i].i; + } + + // Move to the next symbol + el += me->wf.block_stride; + pos += num_shift; + } +} +#endif diff --git a/common/monitor.h b/common/monitor.h index c094285..7440362 100644 --- a/common/monitor.h +++ b/common/monitor.h @@ -39,6 +39,11 @@ typedef struct // KISS FFT housekeeping variables void* fft_work; ///< Work area required by Kiss FFT kiss_fftr_cfg fft_cfg; ///< Kiss FFT housekeeping object +#ifdef WATERFALL_USE_PHASE + int nifft; ///< iFFT size + void* ifft_work; ///< Work area required by inverse Kiss FFT + kiss_fft_cfg ifft_cfg; ///< Inverse Kiss FFT housekeeping object +#endif } monitor_t; void monitor_init(monitor_t* me, const monitor_config_t* cfg); @@ -46,6 +51,10 @@ void monitor_reset(monitor_t* me); 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); +#endif + #ifdef __cplusplus } #endif diff --git a/demo/decode_ft8.c b/demo/decode_ft8.c index 57fce1c..85c3364 100644 --- a/demo/decode_ft8.c +++ b/demo/decode_ft8.c @@ -132,7 +132,7 @@ void decode(const monitor_t* mon) const waterfall_t* wf = &mon->wf; // Find top candidates by Costas sync score and localize them in time and frequency candidate_t candidate_list[kMax_candidates]; - int num_candidates = ft8_find_sync(wf, kMax_candidates, candidate_list, kMin_score); + int num_candidates = ftx_find_sync(wf, kMax_candidates, candidate_list, kMin_score); // Hash table for decoded messages (to check for duplicates) int num_decoded = 0; @@ -153,9 +153,22 @@ void decode(const monitor_t* mon) 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; +#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); +#endif + ftx_message_t message; decode_status_t status; - if (!ft8_decode(wf, cand, kLDPC_iterations, &message, &status)) + if (!ftx_decode(wf, cand, kLDPC_iterations, &message, &status)) { // float snr = cand->score * 0.5f; // TODO: compute better approximation of SNR // printf("000000 %2.1f %+4.2f %4.0f ~ %s\n", snr, time_sec, freq_hz, "---"); @@ -297,7 +310,7 @@ int main(int argc, char** argv) int sample_rate = 12000; int num_samples = slot_time * sample_rate; float signal[num_samples]; - bool isContinuous = false; + bool is_live = false; if (wav_path != NULL) { @@ -314,7 +327,7 @@ int main(int argc, char** argv) audio_init(); audio_open(dev_name); num_samples = (slot_time - 0.4f) * sample_rate; - isContinuous = true; + is_live = true; } // Compute FFT over the whole signal and store it @@ -335,7 +348,7 @@ int main(int argc, char** argv) do { - if (dev_name != NULL) + if (is_live) { // Wait for the start of time slot while (true) @@ -374,7 +387,7 @@ int main(int argc, char** argv) // Reset internal variables for the next time slot monitor_reset(&mon); - } while (isContinuous); + } while (is_live); monitor_free(&mon); diff --git a/ft8/decode.c b/ft8/decode.c index 6125234..48d466e 100644 --- a/ft8/decode.c +++ b/ft8/decode.c @@ -31,11 +31,11 @@ static void heapify_down(candidate_t heap[], int heap_size); static void heapify_up(candidate_t heap[], int heap_size); static void ftx_normalize_logl(float* log174); -static void ft4_extract_symbol(const uint8_t* wf, float* logl); -static void ft8_extract_symbol(const uint8_t* wf, float* logl); -static void ft8_decode_multi_symbols(const uint8_t* wf, int num_bins, int n_syms, int bit_idx, float* log174); +static void ft4_extract_symbol(const WF_ELEM_T* wf, float* logl); +static void ft8_extract_symbol(const WF_ELEM_T* wf, float* logl); +static void ft8_decode_multi_symbols(const WF_ELEM_T* wf, int num_bins, int n_syms, int bit_idx, float* log174); -static const uint8_t* get_cand_mag(const waterfall_t* wf, const candidate_t* candidate) +static const WF_ELEM_T* get_cand_mag(const waterfall_t* wf, const candidate_t* candidate) { int offset = candidate->time_offset; offset = (offset * wf->time_osr) + candidate->time_sub; @@ -44,70 +44,13 @@ static const uint8_t* get_cand_mag(const waterfall_t* wf, const candidate_t* can return wf->mag + offset; } -int ft8_snr(const waterfall_t* wf, const candidate_t* candidate) -{ - int sum_signal = 0; - int sum_noise = 0; - int num_average = 0; - - // Get the pointer to symbol 0 of the candidate - const uint8_t* mag_cand = get_cand_mag(wf, candidate); - - if (wf->protocol == FTX_PROTOCOL_FT4) - { - } - - // Compute average score over sync symbols (m+k = 0-7, 36-43, 72-79) - for (int block = 0; block < FT8_NN; ++block) - { - int block_abs = candidate->time_offset + block; // relative to the captured signal - // Check for time boundaries - if (block_abs < 0) - continue; - if (block_abs >= wf->num_blocks) - break; - - // Get the pointer to symbol 'block' of the candidate - const uint8_t* p8 = mag_cand + (block * wf->block_stride); - - int k = block % FT8_SYNC_OFFSET; - int sm = -1; - if (k < 7) - { - // Check only the neighbors of the expected symbol frequency- and time-wise - sm = kFT8_Costas_pattern[k]; // Index of the expected bin - } - else - { - int max; - for (int m = 0; m < 8; ++m) - { - if ((sm == -1) || (p8[m] > max)) - { - sm = m; - max = p8[m]; - } - } - } - - if (sm != -1) - { - sum_signal += p8[sm]; - sum_noise += (3 + (int)p8[0] + (int)p8[1] + (int)p8[2] + (int)p8[3] + (int)p8[4] + (int)p8[5] + (int)p8[6] + (int)p8[7] - (int)p8[sm]) / 7; - ++num_average; - } - } - // return num_average; - return (sum_signal - sum_noise) / num_average; -} - static int ft8_sync_score(const waterfall_t* wf, const candidate_t* candidate) { int score = 0; int num_average = 0; // Get the pointer to symbol 0 of the candidate - const uint8_t* mag_cand = get_cand_mag(wf, candidate); + const WF_ELEM_T* mag_cand = get_cand_mag(wf, candidate); // Compute average score over sync symbols (m+k = 0-7, 36-43, 72-79) for (int m = 0; m < FT8_NUM_SYNC; ++m) @@ -123,7 +66,7 @@ static int ft8_sync_score(const waterfall_t* wf, const candidate_t* candidate) break; // Get the pointer to symbol 'block' of the candidate - const uint8_t* p8 = mag_cand + (block * wf->block_stride); + const WF_ELEM_T* p8 = mag_cand + (block * wf->block_stride); // Weighted difference between the expected and all other symbols // Does not work as well as the alternative score below @@ -137,25 +80,25 @@ static int ft8_sync_score(const waterfall_t* wf, const candidate_t* candidate) if (sm > 0) { // look at one frequency bin lower - score += p8[sm] - p8[sm - 1]; + score += WF_ELEM_MAG_INT(p8[sm]) - WF_ELEM_MAG_INT(p8[sm - 1]); ++num_average; } if (sm < 7) { // look at one frequency bin higher - score += p8[sm] - p8[sm + 1]; + score += WF_ELEM_MAG_INT(p8[sm]) - WF_ELEM_MAG_INT(p8[sm + 1]); ++num_average; } if ((k > 0) && (block_abs > 0)) { // look one symbol back in time - score += p8[sm] - p8[sm - wf->block_stride]; + score += WF_ELEM_MAG_INT(p8[sm]) - WF_ELEM_MAG_INT(p8[sm - wf->block_stride]); ++num_average; } if (((k + 1) < FT8_LENGTH_SYNC) && ((block_abs + 1) < wf->num_blocks)) { // look one symbol forward in time - score += p8[sm] - p8[sm + wf->block_stride]; + score += WF_ELEM_MAG_INT(p8[sm]) - WF_ELEM_MAG_INT(p8[sm + wf->block_stride]); ++num_average; } } @@ -173,7 +116,7 @@ static int ft4_sync_score(const waterfall_t* wf, const candidate_t* candidate) int num_average = 0; // Get the pointer to symbol 0 of the candidate - const uint8_t* mag_cand = get_cand_mag(wf, candidate); + const WF_ELEM_T* mag_cand = get_cand_mag(wf, candidate); // Compute average score over sync symbols (block = 1-4, 34-37, 67-70, 100-103) for (int m = 0; m < FT4_NUM_SYNC; ++m) @@ -189,7 +132,7 @@ static int ft4_sync_score(const waterfall_t* wf, const candidate_t* candidate) break; // Get the pointer to symbol 'block' of the candidate - const uint8_t* p4 = mag_cand + (block * wf->block_stride); + const WF_ELEM_T* p4 = mag_cand + (block * wf->block_stride); int sm = kFT4_Costas_pattern[m][k]; // Index of the expected bin @@ -200,25 +143,25 @@ static int ft4_sync_score(const waterfall_t* wf, const candidate_t* candidate) if (sm > 0) { // look at one frequency bin lower - score += p4[sm] - p4[sm - 1]; + score += WF_ELEM_MAG_INT(p4[sm]) - WF_ELEM_MAG_INT(p4[sm - 1]); ++num_average; } if (sm < 3) { // look at one frequency bin higher - score += p4[sm] - p4[sm + 1]; + score += WF_ELEM_MAG_INT(p4[sm]) - WF_ELEM_MAG_INT(p4[sm + 1]); ++num_average; } if ((k > 0) && (block_abs > 0)) { // look one symbol back in time - score += p4[sm] - p4[sm - wf->block_stride]; + score += WF_ELEM_MAG_INT(p4[sm]) - WF_ELEM_MAG_INT(p4[sm - wf->block_stride]); ++num_average; } if (((k + 1) < FT4_LENGTH_SYNC) && ((block_abs + 1) < wf->num_blocks)) { // look one symbol forward in time - score += p4[sm] - p4[sm + wf->block_stride]; + score += WF_ELEM_MAG_INT(p4[sm]) - WF_ELEM_MAG_INT(p4[sm + wf->block_stride]); ++num_average; } } @@ -230,8 +173,11 @@ static int ft4_sync_score(const waterfall_t* wf, const candidate_t* candidate) return score; } -int ft8_find_sync(const waterfall_t* wf, int num_candidates, candidate_t heap[], int min_score) +int ftx_find_sync(const waterfall_t* wf, int num_candidates, candidate_t heap[], int min_score) { + int (*sync_fun)(const waterfall_t*, const candidate_t*) = (wf->protocol == FTX_PROTOCOL_FT4) ? ft4_sync_score : ft8_sync_score; + int num_tones = (wf->protocol == FTX_PROTOCOL_FT4) ? 4 : 8; + int heap_size = 0; candidate_t candidate; @@ -242,19 +188,11 @@ int ft8_find_sync(const waterfall_t* wf, int num_candidates, candidate_t heap[], { for (candidate.freq_sub = 0; candidate.freq_sub < wf->freq_osr; ++candidate.freq_sub) { - for (candidate.time_offset = -12; candidate.time_offset < 24; ++candidate.time_offset) + for (candidate.time_offset = -10; candidate.time_offset < 20; ++candidate.time_offset) { - for (candidate.freq_offset = 0; (candidate.freq_offset + 7) < wf->num_bins; ++candidate.freq_offset) + for (candidate.freq_offset = 0; (candidate.freq_offset + num_tones - 1) < wf->num_bins; ++candidate.freq_offset) { - if (wf->protocol == FTX_PROTOCOL_FT4) - { - candidate.score = ft4_sync_score(wf, &candidate); - } - else - { - candidate.score = ft8_sync_score(wf, &candidate); - // candidate.score = ft8_snr(wf, &candidate); - } + candidate.score = sync_fun(wf, &candidate); if (candidate.score < min_score) continue; @@ -300,7 +238,7 @@ int ft8_find_sync(const waterfall_t* wf, int num_candidates, candidate_t heap[], static void ft4_extract_likelihood(const waterfall_t* wf, const candidate_t* cand, float* log174) { - const uint8_t* mag_cand = get_cand_mag(wf, cand); + const WF_ELEM_T* mag = get_cand_mag(wf, cand); // Pointer to 4 magnitude bins of the first symbol // Go over FSK tones and skip Costas sync symbols for (int k = 0; k < FT4_ND; ++k) @@ -319,17 +257,14 @@ static void ft4_extract_likelihood(const waterfall_t* wf, const candidate_t* can } else { - // Pointer to 4 bins of the current symbol - const uint8_t* ps = mag_cand + (sym_idx * wf->block_stride); - - ft4_extract_symbol(ps, log174 + bit_idx); + ft4_extract_symbol(mag + (sym_idx * wf->block_stride), log174 + bit_idx); } } } static void ft8_extract_likelihood(const waterfall_t* wf, const candidate_t* cand, float* log174) { - const uint8_t* mag_cand = get_cand_mag(wf, cand); + const WF_ELEM_T* mag = get_cand_mag(wf, cand); // Pointer to 8 magnitude bins of the first symbol // Go over FSK tones and skip Costas sync symbols for (int k = 0; k < FT8_ND; ++k) @@ -349,10 +284,7 @@ static void ft8_extract_likelihood(const waterfall_t* wf, const candidate_t* can } else { - // Pointer to 8 bins of the current symbol - const uint8_t* ps = mag_cand + (sym_idx * wf->block_stride); - - ft8_extract_symbol(ps, log174 + bit_idx); + ft8_extract_symbol(mag + (sym_idx * wf->block_stride), log174 + bit_idx); } } } @@ -378,7 +310,7 @@ static void ftx_normalize_logl(float* log174) } } -bool ft8_decode(const waterfall_t* wf, const candidate_t* cand, int max_iterations, ftx_message_t* message, decode_status_t* status) +bool ftx_decode(const waterfall_t* wf, const candidate_t* cand, int max_iterations, ftx_message_t* message, decode_status_t* status) { float log174[FTX_LDPC_N]; // message bits encoded as likelihood if (wf->protocol == FTX_PROTOCOL_FT4) @@ -505,14 +437,14 @@ static void heapify_up(candidate_t heap[], int heap_size) } // Compute unnormalized log likelihood log(p(1) / p(0)) of 2 message bits (1 FSK symbol) -static void ft4_extract_symbol(const uint8_t* wf, float* logl) +static void ft4_extract_symbol(const WF_ELEM_T* wf, float* logl) { // Cleaned up code for the simple case of n_syms==1 float s2[4]; for (int j = 0; j < 4; ++j) { - s2[j] = (float)wf[kFT4_Gray_map[j]]; + s2[j] = WF_ELEM_MAG(wf[kFT4_Gray_map[j]]); } logl[0] = max2(s2[2], s2[3]) - max2(s2[0], s2[1]); @@ -520,23 +452,49 @@ static void ft4_extract_symbol(const uint8_t* wf, float* logl) } // Compute unnormalized log likelihood log(p(1) / p(0)) of 3 message bits (1 FSK symbol) -static void ft8_extract_symbol(const uint8_t* wf, float* logl) +static void ft8_extract_symbol(const WF_ELEM_T* wf, float* logl) { // Cleaned up code for the simple case of n_syms==1 +#if 1 float s2[8]; for (int j = 0; j < 8; ++j) { - s2[j] = (float)wf[kFT8_Gray_map[j]]; + s2[j] = WF_ELEM_MAG(wf[kFT8_Gray_map[j]]); } logl[0] = max4(s2[4], s2[5], s2[6], s2[7]) - max4(s2[0], s2[1], s2[2], s2[3]); logl[1] = max4(s2[2], s2[3], s2[6], s2[7]) - max4(s2[0], s2[1], s2[4], s2[5]); logl[2] = max4(s2[1], s2[3], s2[5], s2[7]) - max4(s2[0], s2[2], s2[4], s2[6]); +#else + float a[7] = { + // (float)wf[7] - (float)wf[0], // 0: p(111) / p(000) + (float)wf[5] - (float)wf[2], // 0: p(100) / p(011) + (float)wf[3] - (float)wf[0], // 1: p(010) / p(000) + (float)wf[6] - (float)wf[3], // 2: p(101) / p(010) + (float)wf[6] - (float)wf[2], // 3: p(101) / p(011) + (float)wf[7] - (float)wf[4], // 4: p(111) / p(110) + (float)wf[4] - (float)wf[1], // 5: p(110) / p(001) + (float)wf[5] - (float)wf[1] // 6: p(100) / p(001) + }; + float k = 1.0f; + + // logl[0] = k * (a[0] + a[2] + a[3] + a[5] + a[6]) / 5; + // logl[1] = k * (a[0] / 4 + (a[1] - a[3]) * 5 / 24 + (a[5] - a[2]) / 6 + (a[4] - a[6]) / 24); + // logl[2] = k * (a[0] / 4 + (a[1] - a[3]) / 24 + (a[2] - a[5]) / 6 + (a[4] - a[6]) * 5 / 24); + logl[0] = k * (a[1] / 6 + a[2] / 3 + a[3] / 6 + a[4] / 6 + a[5] / 3 + a[6] / 6); + logl[1] = k * (-a[0] / 4 + a[1] * 7 / 24 + (a[4] - a[3]) / 8 + a[5] / 3 + a[6] / 24); + logl[2] = k * (-a[0] / 4 + (a[1] - a[6]) / 8 + a[2] / 3 + a[3] / 24 + a[4] * 7 / 24 - a[5] * 5 / 18); +#endif + // for (int i = 0; i < 8; ++i) + // printf("%d ", WF_ELEM_MAG_INT(wf[i])); + // for (int i = 0; i < 3; ++i) + // printf("%.1f ", logl[i]); + // printf("\n"); } // Compute unnormalized log likelihood log(p(1) / p(0)) of bits corresponding to several FSK symbols at once -static void ft8_decode_multi_symbols(const uint8_t* wf, int num_bins, int n_syms, int bit_idx, float* log174) +static void ft8_decode_multi_symbols(const WF_ELEM_T* wf, int num_bins, int n_syms, int bit_idx, float* log174) { const int n_bits = 3 * n_syms; const int n_tones = (1 << n_bits); @@ -548,20 +506,20 @@ static void ft8_decode_multi_symbols(const uint8_t* wf, int num_bins, int n_syms int j1 = j & 0x07; if (n_syms == 1) { - s2[j] = (float)wf[kFT8_Gray_map[j1]]; + s2[j] = WF_ELEM_MAG(wf[kFT8_Gray_map[j1]]); continue; } int j2 = (j >> 3) & 0x07; if (n_syms == 2) { - s2[j] = (float)wf[kFT8_Gray_map[j2]]; - s2[j] += (float)wf[kFT8_Gray_map[j1] + 4 * num_bins]; + s2[j] = WF_ELEM_MAG(wf[kFT8_Gray_map[j2]]); + s2[j] += WF_ELEM_MAG(wf[kFT8_Gray_map[j1] + 4 * num_bins]); continue; } int j3 = (j >> 6) & 0x07; - s2[j] = (float)wf[kFT8_Gray_map[j3]]; - s2[j] += (float)wf[kFT8_Gray_map[j2] + 4 * num_bins]; - s2[j] += (float)wf[kFT8_Gray_map[j1] + 8 * num_bins]; + s2[j] = WF_ELEM_MAG(wf[kFT8_Gray_map[j3]]); + s2[j] += WF_ELEM_MAG(wf[kFT8_Gray_map[j2] + 4 * num_bins]); + s2[j] += WF_ELEM_MAG(wf[kFT8_Gray_map[j1] + 8 * num_bins]); } // Extract bit significance (and convert them to float) diff --git a/ft8/decode.h b/ft8/decode.h index 81237c0..064b588 100644 --- a/ft8/decode.h +++ b/ft8/decode.h @@ -12,7 +12,25 @@ extern "C" { #endif -/// Input structure to ft8_find_sync() function. This structure describes stored waterfall data over the whole message slot. +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. @@ -25,12 +43,12 @@ typedef struct 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 - uint8_t* mag; ///< FFT magnitudes stored as uint8_t[blocks][time_osr][freq_osr][num_bins] + 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 } waterfall_t; -/// Output structure of ft8_find_sync() and input structure of ft8_decode(). +/// 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 { @@ -39,12 +57,13 @@ typedef struct 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 - int16_t snr; } candidate_t; /// Structure that contains the status of various steps during decoding of a message typedef struct { + float freq; + float time; int ldpc_errors; ///< Number of LDPC errors during decoding uint16_t crc_extracted; ///< CRC value recovered from the message uint16_t crc_calculated; ///< CRC value calculated over the payload @@ -59,7 +78,7 @@ typedef struct /// @param[in,out] heap Array of 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 ft8_find_sync(const waterfall_t* power, int num_candidates, candidate_t heap[], int min_score); +int ftx_find_sync(const waterfall_t* power, int num_candidates, candidate_t heap[], int min_score); /// 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 @@ -68,7 +87,7 @@ int ft8_find_sync(const waterfall_t* power, int num_candidates, candidate_t heap /// @param[out] message ftx_message_t structure that will receive the decoded message /// @param[out] status 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 ft8_decode(const waterfall_t* power, const candidate_t* cand, int max_iterations, ftx_message_t* message, decode_status_t* status); +bool ftx_decode(const waterfall_t* power, const candidate_t* cand, int max_iterations, ftx_message_t* message, decode_status_t* status); #ifdef __cplusplus }