Reorganise and tune downsampling

yl3jg/downsampling-wip
Karlis Goba 2025-03-15 07:44:23 +02:00
rodzic 50ee0c0636
commit d46d46edaf
13 zmienionych plików z 357 dodań i 138 usunięć

3
.gitignore vendored
Wyświetl plik

@ -2,8 +2,9 @@ gen_ft8
decode_ft8
test_ft8
libft8.a
wsjtx2/
wsjtx/
.build/
.env/
.DS_Store
.vscode/
__pycache__/

Wyświetl plik

@ -1,3 +1,5 @@
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif
#define TWO_PI 6.28318530717958647692

Wyświetl plik

@ -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)

Wyświetl plik

@ -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

Wyświetl plik

@ -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)

Wyświetl plik

@ -4,8 +4,8 @@
#include <math.h>
#include <stdbool.h>
#include "common/common.h"
#include "common/wave.h"
#include <common/common.h>
#include <common/wave.h>
#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;
}

Wyświetl plik

@ -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:

Wyświetl plik

@ -5,13 +5,14 @@
#include <stdbool.h>
#include <math.h>
#include <complex.h>
// #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-ORed with [a] pseudorandom sequence before computing the CRC and FEC parity bits'

Wyświetl plik

@ -5,6 +5,7 @@
#include <stdbool.h>
#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
}

Wyświetl plik

@ -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];

Wyświetl plik

@ -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
}

180
ft8/waterfall.c 100644
Wyświetl plik

@ -0,0 +1,180 @@
#include "waterfall.h"
#include <math.h>
#include <complex.h>
#include <common/common.h>
#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);
}

72
ft8/waterfall.h 100644
Wyświetl plik

@ -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_