VE7UWU
/
sdrangel
mirror of https://github.com/f4exb/sdrangel.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity
sdrangel/ft8/ft8.h

679 lines
21 KiB
C++
Raw Blame History

///////////////////////////////////////////////////////////////////////////////////
// Copyright (C) 2023 Edouard Griffiths, F4EXB. //
// //
// This is the code from ft8mon: https://github.com/rtmrtmrtmrtm/ft8mon //
// written by Robert Morris, AB1HL //
// reformatted and adapted to Qt and SDRangel context //
// //
// This program is free software; you can redistribute it and/or modify //
// it under the terms of the GNU General Public License as published by //
// the Free Software Foundation as version 3 of the License, or //
// (at your option) any later version. //
// //
// This program is distributed in the hope that it will be useful, //
// but WITHOUT ANY WARRANTY; without even the implied warranty of //
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the //
// GNU General Public License V3 for more details. //
// //
// You should have received a copy of the GNU General Public License //
// along with this program. If not, see <http://www.gnu.org/licenses/>. //
///////////////////////////////////////////////////////////////////////////////////
#ifndef ft8_h
#define ft8_h
#include <thread>
#include <mutex>
#include "fft.h"
#include "export.h"
namespace FT8 {
// Callback interface to get the results
class FT8_API CallbackInterface
{
public:
virtual int hcb(
int *a91,
float hz0,
float off,
const char *,
float snr,
int pass,
int correct_bits
) = 0; //!< virtual nathod called each time there is a result
};
//
// manage statistics for soft decoding, to help
// decide how likely each symbol is to be correct,
// to drive LDPC decoding.
//
// meaning of the how (problt_how) parameter:
// 0: gaussian
// 1: index into the actual distribution
// 2: do something complex for the tails.
// 3: index into the actual distribution plus gaussian for tails.
// 4: similar to 3.
// 5: laplace
//
class FT8_API Stats
{
public:
std::vector<float> a_;
float sum_;
bool finalized_;
float mean_; // cached
float stddev_; // cached
float b_; // cached
int how_;
public:
Stats(int how, float log_tail, float log_rate);
void add(float x);
void finalize();
float mean();
float stddev();
// fraction of distribution that's less than x.
// assumes normal distribution.
// this is PHI(x), or the CDF at x,
// or the integral from -infinity
// to x of the PDF.
float gaussian_problt(float x);
// https://en.wikipedia.org/wiki/Laplace_distribution
// m and b from page 116 of Mark Owen's Practical Signal Processing.
float laplace_problt(float x);
// look into the actual distribution.
float problt(float x);
private:
float log_tail_;
float log_rate_;
};
class FT8_API Strength
{
public:
float hz_;
int off_;
float strength_; // higher is better
};
// same as Python class CDECODE
//
struct FT8_API cdecode
{
float hz0;
float hz1;
float off;
int *bits; // 174
};
// 1920-point FFT at 12000 samples/second
// 6.25 Hz spacing, 0.16 seconds/symbol
// encode chain:
// 77 bits
// append 14 bits CRC (for 91 bits)
// LDPC(174,91) yields 174 bits
// that's 58 3-bit FSK-8 symbols
// gray code each 3 bits
// insert three 7-symbol Costas sync arrays
// at symbol #s 0, 36, 72 of final signal
// thus: 79 FSK-8 symbols
// total transmission time is 12.64 seconds
// tunable parameters
struct FT8_API FT8Params
{
int nthreads; // number of parallel threads, for multi-core
int npasses_one; // number of spectral subtraction passes
int npasses_two; // number of spectral subtraction passes
int ldpc_iters; // how hard LDPC decoding should work
int snr_win; // averaging window, in symbols, for SNR conversion
int snr_how; // technique to measure "N" for SNR. 0 means median of the 8 tones.
float shoulder200; // for 200 sps bandpass filter
float shoulder200_extra; // for bandpass filter
float second_hz_win; // +/- hz
int second_hz_n; // divide total window into this many pieces
float second_off_win; // +/- search window in symbol-times
int second_off_n;
int third_hz_n;
float third_hz_win;
int third_off_n;
float third_off_win;
float log_tail;
float log_rate;
int problt_how_noise;
int problt_how_sig;
int use_apriori;
int use_hints; // 1 means use all hints, 2 means just CQ hints
int win_type;
int osd_depth; // 6; // don't increase beyond 6, produces too much garbage
int osd_ldpc_thresh; // demand this many correct LDPC parity bits before OSD
int ncoarse; // number of offsets per hz produced by coarse()
int ncoarse_blocks;
float tminus; // start looking at 0.5 - tminus seconds
float tplus;
int coarse_off_n;
int coarse_hz_n;
float already_hz;
float overlap;
int overlap_edges;
float nyquist;
int oddrate;
float pass0_frac;
int reduce_how;
float go_extra;
int do_reduce;
int pass_threshold;
int strength_how;
int known_strength_how;
int coarse_strength_how;
float reduce_shoulder;
float reduce_factor;
float reduce_extra;
float coarse_all;
int second_count;
int soft_phase_win;
float subtract_ramp;
int soft_ones;
int soft_pairs;
int soft_triples;
int do_second;
int do_fine_hz;
int do_fine_off;
int do_third;
float fine_thresh;
int fine_max_off;
int fine_max_tone;
int known_sparse;
float c_soft_weight;
int c_soft_win;
int bayes_how;
FT8Params()
{
nthreads = 8; // number of parallel threads, for multi-core
npasses_one = 3; // number of spectral subtraction passes
npasses_two = 3; // number of spectral subtraction passes
ldpc_iters = 25; // how hard LDPC decoding should work
snr_win = 7; // averaging window, in symbols, for SNR conversion
snr_how = 3; // technique to measure "N" for SNR. 0 means median of the 8 tones.
shoulder200 = 10; // for 200 sps bandpass filter
shoulder200_extra = 0.0; // for bandpass filter
second_hz_win = 3.5; // +/- hz
second_hz_n = 8; // divide total window into this many pieces
second_off_win = 0.5; // +/- search window in symbol-times
second_off_n = 10;
third_hz_n = 3;
third_hz_win = 0.25;
third_off_n = 4;
third_off_win = 0.075;
log_tail = 0.1;
log_rate = 8.0;
problt_how_noise = 0;
problt_how_sig = 0;
use_apriori = 1;
use_hints = 2; // 1 means use all hints, 2 means just CQ hints
win_type = 1;
osd_depth = 0; // 6; // don't increase beyond 6, produces too much garbage
osd_ldpc_thresh = 70; // demand this many correct LDPC parity bits before OSD
ncoarse = 1; // number of offsets per hz produced by coarse()
ncoarse_blocks = 1;
tminus = 2.2; // start looking at 0.5 - tminus seconds
tplus = 2.4;
coarse_off_n = 4;
coarse_hz_n = 4;
already_hz = 27;
overlap = 20;
overlap_edges = 0;
nyquist = 0.925;
oddrate = 1;
pass0_frac = 1.0;
reduce_how = 2;
go_extra = 3.5;
do_reduce = 1;
pass_threshold = 1;
strength_how = 4;
known_strength_how = 7;
coarse_strength_how = 6;
reduce_shoulder = -1;
reduce_factor = 0.25;
reduce_extra = 0;
coarse_all = -1;
second_count = 3;
soft_phase_win = 2;
subtract_ramp = 0.11;
soft_ones = 2;
soft_pairs = 1;
soft_triples = 1;
do_second = 1;
do_fine_hz = 1;
do_fine_off = 1;
do_third = 2;
fine_thresh = 0.19;
fine_max_off = 2;
fine_max_tone = 4;
known_sparse = 1;
c_soft_weight = 7;
c_soft_win = 2;
bayes_how = 1;
}
}; // class FT8Params
// The FT8 worker
class FT8_API FT8
{
public:
std::thread *th_;
float min_hz_;
float max_hz_;
std::vector<float> samples_; // input to each pass
std::vector<float> nsamples_; // subtract from here
int start_; // sample number of 0.5 seconds into samples[]
int rate_; // samples/second
double deadline_; // start time + budget
double final_deadline_; // keep going this long if no decodes
std::vector<int> hints1_;
std::vector<int> hints2_;
int pass_;
float down_hz_;
std::mutex cb_mu_;
CallbackInterface *cb_; // call-back interface
std::mutex hack_mu_;
int hack_size_;
int hack_off_;
int hack_len_;
float hack_0_;
float hack_1_;
const float *hack_data_;
std::vector<std::complex<float>> hack_bins_;
std::vector<cdecode> prevdecs_;
FFTEngine::Plan *plan32_;
FT8(
const std::vector<float> &samples,
float min_hz,
float max_hz,
int start,
int rate,
int hints1[],
int hints2[],
double deadline,
double final_deadline,
CallbackInterface *cb,
std::vector<cdecode> prevdecs,
FFTEngine *fftEngine
);
~FT8();
// strength of costas block of signal with tone 0 at bi0,
// and symbol zero at si0.
float one_coarse_strength(const FFTEngine::ffts_t &bins, int bi0, int si0);
// return symbol length in samples at the given rate.
// insist on integer symbol lengths so that we can
// use whole FFT bins.
int blocksize(int rate);
//
// look for potential signals by searching FFT bins for Costas symbol
// blocks. returns a vector of candidate positions.
//
std::vector<Strength> coarse(const FFTEngine::ffts_t &bins, int si0, int si1);
//
// reduce the sample rate from arate to brate.
// center hz0..hz1 in the new nyquist range.
// but first filter to that range.
// sets delta_hz to hz moved down.
//
std::vector<float> reduce_rate(
const std::vector<float> &a,
float hz0,
float hz1,
int arate,
int brate,
float &delta_hz
);
// The actual main process
void go(int npasses);
//
// what's the strength of the Costas sync blocks of
// the signal starting at hz and off?
//
float one_strength(const std::vector<float> &samples200, float hz, int off);
//
// given a complete known signal's symbols in syms,
// how strong is it? used to look for the best
// offset and frequency at which to subtract a
// decoded signal.
//
float one_strength_known(
const std::vector<float> &samples,
int rate,
const std::vector<int> &syms,
float hz,
int off
);
int search_time_fine(
const std::vector<float> &samples200,
int off0,
int offN,
float hz,
int gran,
float &str
);
int search_time_fine_known(
const std::vector<std::complex<float>> &bins,
int rate,
const std::vector<int> &syms,
int off0,
int offN,
float hz,
int gran,
float &str
);
//
// search for costas blocks in an MxN time/frequency grid.
// hz0 +/- hz_win in hz_inc increments. hz0 should be near 25.
// off0 +/- off_win in off_inc incremenents.
//
std::vector<Strength> search_both(
const std::vector<float> &samples200,
float hz0,
int hz_n,
float hz_win,
int off0,
int off_n,
int off_win
);
void search_both_known(
const std::vector<float> &samples,
int rate,
const std::vector<int> &syms,
float hz0,
float off_secs0, // seconds
float &hz_out,
float &off_out
);
//
// shift frequency by shifting the bins of one giant FFT.
// so no problem with phase mismatch &c at block boundaries.
// surprisingly fast at 200 samples/second.
// shifts *down* by hz.
//
std::vector<float> fft_shift(
const std::vector<float> &samples,
int off,
int len,
int rate,
float hz
);
//
// shift down by hz.
//
std::vector<float> fft_shift_f(
const std::vector<std::complex<float>> &bins,
int rate,
float hz
);
// shift the frequency by a fraction of 6.25,
// to center hz on bin 4 (25 hz).
std::vector<float> shift200(
const std::vector<float> &samples200,
int off,
int len,
float hz
);
// returns a mini-FFT of 79 8-tone symbols.
FFTEngine::ffts_t extract(const std::vector<float> &samples200, float, int off);
//
// m79 is a 79x8 array of complex.
//
FFTEngine::ffts_t un_gray_code_c(const FFTEngine::ffts_t &m79);
//
// m79 is a 79x8 array of float.
//
std::vector<std::vector<float>> un_gray_code_r(const std::vector<std::vector<float>> &m79);
//
// normalize levels by windowed median.
// this helps, but why?
//
std::vector<std::vector<float>> convert_to_snr(const std::vector<std::vector<float>> &m79);
//
// normalize levels by windowed median.
// this helps, but why?
//
std::vector<std::vector<std::complex<float>>> c_convert_to_snr(
const std::vector<std::vector<std::complex<float>>> &m79
);
//
// statistics to decide soft probabilities,
// to drive LDPC decoder.
// distribution of strongest tones, and
// distribution of noise.
//
void make_stats(
const std::vector<std::vector<float>> &m79,
Stats &bests,
Stats &all
);
//
// convert 79x8 complex FFT bins to magnitudes.
//
// exploits local phase coherence by decreasing magnitudes of bins
// whose phase is far from the phases of nearby strongest tones.
//
// relies on each tone being reasonably well centered in its FFT bin
// (in time and frequency) so that each tone completes an integer
// number of cycles and thus preserves phase from one symbol to the
// next.
//
std::vector<std::vector<float>> soft_c2m(const FFTEngine::ffts_t &c79);
//
// guess the probability that a bit is zero vs one,
// based on strengths of strongest tones that would
// give it those values. for soft LDPC decoding.
//
// returns log-likelihood, zero is positive, one is negative.
//
float bayes(
float best_zero,
float best_one,
int lli,
Stats &bests,
Stats &all
);
//
// c79 is 79x8 complex tones, before un-gray-coding.
//
void soft_decode(const FFTEngine::ffts_t &c79, float ll174[]);
//
// c79 is 79x8 complex tones, before un-gray-coding.
//
void c_soft_decode(const FFTEngine::ffts_t &c79x, float ll174[]);
//
// turn 79 symbol numbers into 174 bits.
// strip out the three Costas sync blocks,
// leaving 58 symbol numbers.
// each represents three bits.
// (all post-un-gray-code).
// str is per-symbol strength; must be positive.
// each returned element is < 0 for 1, > 0 for zero,
// scaled by str.
//
std::vector<float> extract_bits(const std::vector<int> &syms, const std::vector<float> str);
// decode successive pairs of symbols. exploits the likelyhood
// that they have the same phase, by summing the complex
// correlations for each possible pair and using the max.
void soft_decode_pairs(
const FFTEngine::ffts_t &m79x,
float ll174[]
);
void soft_decode_triples(
const FFTEngine::ffts_t &m79x,
float ll174[]
);
//
// given log likelyhood for each bit, try LDPC and OSD decoders.
// on success, puts corrected 174 bits into a174[].
//
int decode(const float ll174[], int a174[], int use_osd, std::string &comment);
//
// bandpass filter some FFT bins.
// smooth transition from stop-band to pass-band,
// so that it's not a brick-wall filter, so that it
// doesn't ring.
//
std::vector<std::complex<float>> fbandpass(
const std::vector<std::complex<float>> &bins0,
float bin_hz,
float low_outer, // start of transition
float low_inner, // start of flat area
float high_inner, // end of flat area
float high_outer // end of transition
);
//
// move hz down to 25, filter+convert to 200 samples/second.
//
// like fft_shift(). one big FFT, move bins down and
// zero out those outside the band, then IFFT,
// then re-sample.
//
// XXX maybe merge w/ fft_shift() / shift200().
//
std::vector<float> down_v7(const std::vector<float> &samples, float hz);
std::vector<float> down_v7_f(const std::vector<std::complex<float>> &bins, int len, float hz);
//
// putative start of signal is at hz and symbol si0.
//
// return 2 if it decodes to a brand-new message.
// return 1 if it decodes but we've already seen it,
// perhaps in a different pass.
// return 0 if we could not decode.
//
// XXX merge with one_iter().
//
int one(const std::vector<std::complex<float>> &bins, int len, float hz, int off);
// return 2 if it decodes to a brand-new message.
// return 1 if it decodes but we've already seen it,
// perhaps in a different pass.
// return 0 if we could not decode.
int one_iter(const std::vector<float> &samples200, int best_off, float hz_for_cb);
//
// estimate SNR, yielding numbers vaguely similar to WSJT-X.
// m79 is a 79x8 complex FFT output.
//
float guess_snr(const FFTEngine::ffts_t &m79);
//
// compare phases of successive symbols to guess whether
// the starting offset is a little too high or low.
// we expect each symbol to have the same phase.
// an error in causes the phase to advance at a steady rate.
// so if hz is wrong, we expect the phase to advance
// or retard at a steady pace.
// an error in offset causes each symbol to start at
// a phase that depends on the symbol's frequency;
// a particular offset error causes a phase error
// that depends on frequency.
// hz0 is actual FFT bin number of m79[...][0] (always 4).
//
// the output adj_hz is relative to the FFT bin center;
// a positive number means the real signal seems to be
// a bit higher in frequency that the bin center.
//
// adj_off is the amount to change the offset, in samples.
// should be subtracted from offset.
//
void fine(const FFTEngine::ffts_t &m79, int, float &adj_hz, float &adj_off);
//
// subtract a corrected decoded signal from nsamples_,
// perhaps revealing a weaker signal underneath,
// to be decoded in a subsequent pass.
//
// re79[] holds the error-corrected symbol numbers.
//
void subtract(
const std::vector<int> re79,
float hz0,
float hz1,
float off_sec
);
//
// decode, give to callback, and subtract.
//
// return 2 if it decodes to a brand-new message.
// return 1 if it decodes but we've already seen it,
// perhaps in a different pass.
// return 0 if we could not decode.
//
int try_decode(
const std::vector<float> &samples200,
float ll174[174],
float best_hz,
int best_off_samples,
float hz0_for_cb,
float,
int use_osd,
const char *comment1,
const FFTEngine::ffts_t &m79
);
//
// given 174 bits corrected by LDPC, work
// backwards to the symbols that must have
// been sent.
// used to help ensure that subtraction subtracts
// at the right place.
//
std::vector<int> recode(int a174[]);
//
// the signal is at roughly 25 hz in samples200.
//
// return 2 if it decodes to a brand-new message.
// return 1 if it decodes but we've already seen it,
// perhaps in a different pass.
// return 0 if we could not decode.
//
int one_iter1(
const std::vector<float> &samples200x,
int best_off,
float best_hz,
float hz0_for_cb,
float hz1_for_cb
);
FT8Params& getParams() { return params; }
private:
FT8Params params;
FFTEngine *fftEngine_;
static const double apriori174[];
}; // class FT8
class FT8_API FT8Decoder {
public:
void entry(
float xsamples[],
int nsamples,
int start,
int rate,
float min_hz,
float max_hz,
int hints1[],
int hints2[],
double time_left,
double total_time_left,
CallbackInterface *cb,
int,
struct cdecode *
);
FT8Params& getParams() { return params; }
private:
FT8Params params;
}; // FT8Decoder
} // namespace FT8
#endif
Reference in New Issue View Git Blame Copy Permalink