VE7UWU/sdrangel
1
0
Fork 0
mirror of https://github.com/f4exb/sdrangel.git synced 2024-11-21 23:55:13 -05:00
Code Issues Projects Releases Wiki Activity

FT8 demod: reorganized code (1)

Browse Source
This commit is contained in:
f4exb 2023-01-10 00:15:46 +01:00
parent 65e9d1a167
commit 7ec86b09de
5 changed files with 3606 additions and 3262 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

36
ft8/fft.cpp
View File

@ -35,39 +35,6 @@ namespace FT8 {
// MEASURE=0, ESTIMATE=64, PATIENT=32
int fftw_type = FFTW_ESTIMATE;
// a cached fftw plan, for both of:
// fftwf_plan_dft_r2c_1d(n, m_in, m_out, FFTW_ESTIMATE);
// fftwf_plan_dft_c2r_1d(n, m_in, m_out, FFTW_ESTIMATE);
class Plan
{
public:
int n_;
int type_;
//
// real -> complex
//
fftwf_complex *c_; // (n_ / 2) + 1 of these
float *r_; // n_ of these
fftwf_plan fwd_; // forward plan
fftwf_plan rev_; // reverse plan
//
// complex -> complex
//
fftwf_complex *cc1_; // n
fftwf_complex *cc2_; // n
fftwf_plan cfwd_; // forward plan
fftwf_plan crev_; // reverse plan
// how much CPU time spent in FFTs that use this plan.
#if TIMING
double time_;
#endif
const char *why_;
int uses_;
};
static std::mutex plansmu;
static Plan *plans[1000];
static int nplans;
@ -622,8 +589,7 @@ std::vector<float> hilbert_shift(const std::vector<float> &x, float hz0, float h
return ret;
}
void
fft_stats()
void fft_stats()
{
for (int i = 0; i < nplans; i++)
{

54
ft8/fft.h
View File

@ -28,18 +28,50 @@
namespace FT8
{
class Plan;
Plan *get_plan(int n, const char *why);
// a cached fftw plan, for both of:
// fftwf_plan_dft_r2c_1d(n, m_in, m_out, FFTW_ESTIMATE);
// fftwf_plan_dft_c2r_1d(n, m_in, m_out, FFTW_ESTIMATE);
class Plan
{
public:
int n_;
int type_;
std::vector<std::complex<float>> one_fft(const std::vector<float> &samples, int i0, int block, const char *why, Plan *p);
std::vector<float> one_ifft(const std::vector<std::complex<float>> &bins, const char *why);
typedef std::vector<std::vector<std::complex<float>>> ffts_t;
ffts_t ffts(const std::vector<float> &samples, int i0, int block, const char *why);
std::vector<std::complex<float>> one_fft_c(const std::vector<float> &samples, int i0, int block, const char *why);
std::vector<std::complex<float>> one_fft_cc(const std::vector<std::complex<float>> &samples, int i0, int block, const char *why);
std::vector<std::complex<float>> one_ifft_cc(const std::vector<std::complex<float>> &bins, const char *why);
std::vector<std::complex<float>> analytic(const std::vector<float> &x, const char *why);
std::vector<float> hilbert_shift(const std::vector<float> &x, float hz0, float hz1, int rate);
//
// real -> complex
//
fftwf_complex *c_; // (n_ / 2) + 1 of these
float *r_; // n_ of these
fftwf_plan fwd_; // forward plan
fftwf_plan rev_; // reverse plan
//
// complex -> complex
//
fftwf_complex *cc1_; // n
fftwf_complex *cc2_; // n
fftwf_plan cfwd_; // forward plan
fftwf_plan crev_; // reverse plan
// how much CPU time spent in FFTs that use this plan.
#if TIMING
double time_;
#endif
const char *why_;
int uses_;
};
Plan *get_plan(int n, const char *why);
std::vector<std::complex<float>> one_fft(const std::vector<float> &samples, int i0, int block, const char *why, Plan *p);
std::vector<float> one_ifft(const std::vector<std::complex<float>> &bins, const char *why);
typedef std::vector<std::vector<std::complex<float>>> ffts_t;
ffts_t ffts(const std::vector<float> &samples, int i0, int block, const char *why);
std::vector<std::complex<float>> one_fft_c(const std::vector<float> &samples, int i0, int block, const char *why);
std::vector<std::complex<float>> one_fft_cc(const std::vector<std::complex<float>> &samples, int i0, int block, const char *why);
std::vector<std::complex<float>> one_ifft_cc(const std::vector<std::complex<float>> &bins, const char *why);
std::vector<std::complex<float>> analytic(const std::vector<float> &x, const char *why);
std::vector<float> hilbert_shift(const std::vector<float> &x, float hz0, float hz1, int rate);
} // namespace FT8

6135
ft8/ft8.cpp
View File

File diff suppressed because it is too large Load Diff

637
ft8/ft8.h
View File

@ -21,6 +21,10 @@
#ifndef ft8_h
#define ft8_h
#include <thread>
#include <mutex>
#include "fft.h"
namespace FT8 {
// Callback interface to get the results
class CallbackInterface
@ -37,6 +41,62 @@ public:
) = 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 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 Strength
{
public:
float hz_;
int off_;
float strength_; // higher is better
};
// same as Python class CDECODE
//
struct cdecode
@ -47,23 +107,568 @@ struct cdecode
int *bits; // 174
};
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 *
);
// 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 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
{
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_;
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
);
~FT8();
// strength of costas block of signal with tone 0 at bi0,
// and symbol zero at si0.
float one_coarse_strength(const 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 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.
ffts_t extract(const std::vector<float> &samples200, float, int off);
//
// m79 is a 79x8 array of complex.
//
ffts_t un_gray_code_c(const 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 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 ffts_t &c79, float ll174[]);
//
// c79 is 79x8 complex tones, before un-gray-coding.
//
void c_soft_decode(const 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 ffts_t &m79x,
float ll174[]
);
void soft_decode_triples(
const 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 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 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 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;
static const double apriori174[];
}; // class FT8
class 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;
};
float set(char *param, char *val);
} // namespace FT8
#endif

6
sdrbench/test_ft8.cpp
View File

@ -167,7 +167,9 @@ void MainBench::testFT8(const QString& wavFile)
wfile.close();
FT8::entry(
FT8::FT8Decoder decoder;
decoder.entry(
samples.data(),
samples.size(),
0.5 * header.m_sampleRate,
@ -180,7 +182,7 @@ void MainBench::testFT8(const QString& wavFile)
budget,
&testft8Callback,
0,
(struct FT8::cdecode *) 0
(struct FT8::cdecode *) nullptr
);
qDebug("MainBench::testFT8: done");
const std::map<std::string, bool>& msgMap = testft8Callback.getMsgMap();