From fbfa04a06dec9c33665d996cd3a8bcdf0cbd52d7 Mon Sep 17 00:00:00 2001 From: f4exb Date: Sat, 7 Mar 2026 11:57:08 +0100 Subject: [PATCH] FT4 demod: copy original FT8 worker to FT4. Isolate AI generated classes --- ft8/ft4.cpp | 3720 ++++++++++++++++++- ft8/ft4.h | 632 +++- ft8/ft4worker.cpp | 14 +- ft8/ft4worker.h | 4 +- ft8/ft8.cpp | 72 - ft8/util.cpp | 72 + ft8/util.h | 5 + plugins/channelrx/demodft8/ft8demodworker.h | 2 +- sdrbench/test_ft4.cpp | 2 +- 9 files changed, 4433 insertions(+), 90 deletions(-) diff --git a/ft8/ft4.cpp b/ft8/ft4.cpp index 404acaafc..33ffc8a06 100644 --- a/ft8/ft4.cpp +++ b/ft8/ft4.cpp @@ -7,21 +7,3596 @@ #include #include #include +#include #include #include +#include "ft4.h" #include "ft4worker.h" #include "fft.h" +#include "osd.h" +#include "libldpc.h" +#include "util.h" +#include "arrays.h" namespace FT8 { -FT4Decoder::~FT4Decoder() +// a-priori probability of each of the 174 LDPC codeword +// bits being one. measured from reconstructed correct +// codewords, into ft8bits, then python bprob.py. +// from ft8-n4 +const double FT4::apriori174[] = { + 0.47, 0.32, 0.29, 0.37, 0.52, 0.36, 0.40, 0.42, 0.42, 0.53, 0.44, + 0.44, 0.39, 0.46, 0.39, 0.38, 0.42, 0.43, 0.45, 0.51, 0.42, 0.48, + 0.31, 0.45, 0.47, 0.53, 0.59, 0.41, 0.03, 0.50, 0.30, 0.26, 0.40, + 0.65, 0.34, 0.49, 0.46, 0.49, 0.69, 0.40, 0.45, 0.45, 0.60, 0.46, + 0.43, 0.49, 0.56, 0.45, 0.55, 0.51, 0.46, 0.37, 0.55, 0.52, 0.56, + 0.55, 0.50, 0.01, 0.19, 0.70, 0.88, 0.75, 0.75, 0.74, 0.73, 0.18, + 0.71, 0.35, 0.60, 0.58, 0.36, 0.60, 0.38, 0.50, 0.02, 0.01, 0.98, + 0.48, 0.49, 0.54, 0.50, 0.49, 0.53, 0.50, 0.49, 0.49, 0.51, 0.51, + 0.51, 0.47, 0.50, 0.53, 0.51, 0.46, 0.51, 0.51, 0.48, 0.51, 0.52, + 0.50, 0.52, 0.51, 0.50, 0.49, 0.53, 0.52, 0.50, 0.46, 0.47, 0.48, + 0.52, 0.50, 0.49, 0.51, 0.49, 0.49, 0.50, 0.50, 0.50, 0.50, 0.51, + 0.50, 0.49, 0.49, 0.55, 0.49, 0.51, 0.48, 0.55, 0.49, 0.48, 0.50, + 0.51, 0.50, 0.51, 0.50, 0.51, 0.53, 0.49, 0.54, 0.50, 0.48, 0.49, + 0.46, 0.51, 0.51, 0.52, 0.49, 0.51, 0.49, 0.51, 0.50, 0.49, 0.50, + 0.50, 0.47, 0.49, 0.52, 0.49, 0.51, 0.49, 0.48, 0.52, 0.48, 0.49, + 0.47, 0.50, 0.48, 0.50, 0.49, 0.51, 0.51, 0.51, 0.49, +}; + +FT4::FT4( + const std::vector &samples, + float min_hz, + float max_hz, + int start, + int rate, + int hints1[], + int hints2[], + double deadline, + double final_deadline, + CallbackInterface *cb, + std::vector prevdecs, + FFTEngine *fftEngine +) +{ + samples_ = samples; + min_hz_ = min_hz; + max_hz_ = max_hz; + prevdecs_ = prevdecs; + start_ = start; + rate_ = rate; + deadline_ = deadline; + final_deadline_ = final_deadline; + cb_ = cb; + down_hz_ = 0; + + for (int i = 0; hints1[i]; i++) { + hints1_.push_back(hints1[i]); + } + + for (int i = 0; hints2[i]; i++) { + hints2_.push_back(hints2[i]); + } + + hack_size_ = -1; + hack_data_ = nullptr; + hack_off_ = -1; + hack_len_ = -1; + fftEngine_ = fftEngine; + npasses_ = 1; +} + +FT4::~FT4() +{ +} + +void FT4::start_work() +{ + go(npasses_); + emit finished(); +} + +// strength of costas block of signal with tone 0 at bi0, +// and symbol zero at si0. +float FT4::one_coarse_strength(const FFTEngine::ffts_t &bins, int bi0, int si0) +{ + int costas[] = {3, 1, 4, 0, 6, 5, 2}; + + // assert(si0 >= 0 && si0 + 72 + 8 <= (int)bins.size()); + // assert(bi0 >= 0 && bi0 + 8 <= (int)bins[0].size()); + + float sig = 0.0; + float noise = 0.0; + + if (params.coarse_all >= 0) + { + for (int si = 0; si < 79; si++) + { + float mx; + int mxi = -1; + float sum = 0; + + for (int i = 0; i < 8; i++) + { + float x = std::abs(bins[si0 + si][bi0 + i]); + sum += x; + + if (mxi < 0 || x > mx) + { + mxi = i; + mx = x; + } + } + + if (si >= 0 && si < 7) + { + float x = std::abs(bins[si0 + si][bi0 + costas[si - 0]]); + sig += x; + noise += sum - x; + } + else if (si >= 36 && si < 36 + 7) + { + float x = std::abs(bins[si0 + si][bi0 + costas[si - 36]]); + sig += x; + noise += sum - x; + } + else if (si >= 72 && si < 72 + 7) + { + float x = std::abs(bins[si0 + si][bi0 + costas[si - 72]]); + sig += x; + noise += sum - x; + } + else + { + sig += params.coarse_all * mx; + noise += params.coarse_all * (sum - mx); + } + } + } + else + { + // coarse_all = -1 + // just costas symbols + for (int si = 0; si < 7; si++) + { + for (int bi = 0; bi < 8; bi++) + { + float x = 0; + x += std::abs(bins[si0 + si][bi0 + bi]); + x += std::abs(bins[si0 + 36 + si][bi0 + bi]); + x += std::abs(bins[si0 + 72 + si][bi0 + bi]); + + if (bi == costas[si]) { + sig += x; + } else { + noise += x; + } + } + } + } + + if (params.coarse_strength_how == 0) { + return sig - noise; + } else if (params.coarse_strength_how == 1) { + return sig - noise / 7; + } else if (params.coarse_strength_how == 2) { + return sig / (noise / 7); + } else if (params.coarse_strength_how == 3) { + return sig / (sig + (noise / 7)); + } else if (params.coarse_strength_how == 4) { + return sig; + } else if (params.coarse_strength_how == 5) { + return sig / (sig + noise); + } else if (params.coarse_strength_how == 6) { + // this is it. + return sig / noise; + } else { + return 0; + } +} + +// return symbol length in samples at the given rate. +// insist on integer symbol lengths so that we can +// use whole FFT bins. +int FT4::blocksize(int rate) +{ + // FT8 symbol length is 1920 at 12000 samples/second. + int xblock = (1920*rate) / 12000; + int block = xblock; + return block; +} + + + +// +// look for potential psignals by searching FFT bins for Costas symbol +// blocks. returns a vector of candidate positions. +// +std::vector FT4::coarse(const FFTEngine::ffts_t &bins, int si0, int si1) +{ + int block = blocksize(rate_); + int nbins = bins[0].size(); + float bin_hz = rate_ / (float)block; + int min_bin = min_hz_ / bin_hz; + int max_bin = max_hz_ / bin_hz; + std::vector strengths; + + for (int bi = min_bin; bi < max_bin && bi + 8 <= nbins; bi++) + { + std::vector sv; + + for (int si = si0; si < si1 && si + 79 < (int)bins.size(); si++) + { + float s = one_coarse_strength(bins, bi, si); + Strength st; + st.strength_ = s; + st.hz_ = bi * 6.25; + st.off_ = si * block; + sv.push_back(st); + } + + if (sv.size() < 1) { + break; + } + + // save best ncoarse offsets, but require that they be separated + // by at least one symbol time. + + std::sort( + sv.begin(), + sv.end(), + [](const Strength &a, const Strength &b) -> bool { + return a.strength_ > b.strength_; + } + ); + + strengths.push_back(sv[0]); + int nn = 1; + + for (int i = 1; nn < params.ncoarse && i < (int)sv.size(); i++) + { + if (std::abs(sv[i].off_ - sv[0].off_) > params.ncoarse_blocks * block) + { + strengths.push_back(sv[i]); + nn++; + } + } + } + + return strengths; +} + +// +// 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 FT4::reduce_rate( + const std::vector &a, + float hz0, + float hz1, + int arate, + int brate, + float &delta_hz +) +{ + // assert(brate < arate); + // assert(hz1 - hz0 <= brate / 2); + + // the pass band is hz0..hz1 + // stop bands are 0..hz00 and hz11..nyquist. + float hz00, hz11; + hz0 = std::max(0.0f, hz0 - params.reduce_extra); + hz1 = std::min(arate / 2.0f, hz1 + params.reduce_extra); + + if (params.reduce_shoulder > 0) + { + hz00 = hz0 - params.reduce_shoulder; + hz11 = hz1 + params.reduce_shoulder; + } + else + { + float mid = (hz0 + hz1) / 2; + hz00 = mid - (brate * params.reduce_factor); + hz00 = std::min(hz00, hz0); + hz11 = mid + (brate * params.reduce_factor); + hz11 = std::max(hz11, hz1); + } + + int alen = a.size(); + std::vector> bins1 = fftEngine_->one_fft(a, 0, alen); + int nbins1 = bins1.size(); + float bin_hz = arate / (float)alen; + + if (params.reduce_how == 2) + { + // band-pass filter the FFT output. + bins1 = fbandpass( + bins1, + bin_hz, + hz00, + hz0, + hz1, + hz11 + ); + } + + if (params.reduce_how == 3) + { + for (int i = 0; i < nbins1; i++) + { + if (i < (hz0 / bin_hz)) { + bins1[i] = 0; + } else if (i > (hz1 / bin_hz)) { + bins1[i] = 0; + } + } + } + + // shift down. + int omid = ((hz0 + hz1) / 2) / bin_hz; + int nmid = (brate / 4.0) / bin_hz; + int delta = omid - nmid; // amount to move down + // assert(delta < nbins1); + int blen = round(alen * (brate / (float)arate)); + std::vector> bbins(blen / 2 + 1); + + for (int i = 0; i < (int)bbins.size(); i++) + { + if (delta > 0) { + bbins[i] = bins1[i + delta]; + } else { + bbins[i] = bins1[i]; + } + } + + // use ifft to reduce the rate. + std::vector vvv = fftEngine_->one_ifft(bbins); + delta_hz = delta * bin_hz; + + return vvv; +} + +void FT4::go(int npasses) +{ + if (0) + { + fprintf(stderr, "go: %.0f .. %.0f, %.0f, rate=%d\n", + min_hz_, max_hz_, max_hz_ - min_hz_, rate_); + } + + // trim to make samples_ a good size for FFTW. + int nice_sizes[] = {18000, 18225, 36000, 36450, + 54000, 54675, 72000, 72900, + 144000, 145800, 216000, 218700, + 0}; + int nice = -1; + + for (int i = 0; nice_sizes[i]; i++) + { + int sz = nice_sizes[i]; + + if (fabs(samples_.size() - sz) < 0.05 * samples_.size()) + { + nice = sz; + break; + } + } + + if (nice != -1) { + samples_.resize(nice); + } + + // assert(min_hz_ >= 0 && max_hz_ + 50 <= rate_ / 2); + + // can we reduce the sample rate? + int nrate = -1; + for (int xrate = 100; xrate < rate_; xrate += 100) + { + if (xrate < rate_ && (params.oddrate || (rate_ % xrate) == 0)) + { + if (((max_hz_ - min_hz_) + 50 + 2 * params.go_extra) < params.nyquist * (xrate / 2)) + { + nrate = xrate; + break; + } + } + } + + if (params.do_reduce && nrate > 0 && nrate < rate_ * 0.75) + { + // filter and reduce the sample rate from rate_ to nrate. + + double t0 = now(); + int osize = samples_.size(); + + float delta_hz; // how much it moved down + samples_ = reduce_rate( + samples_, + min_hz_ - 3.1 - params.go_extra, + max_hz_ + 50 - 3.1 + params.go_extra, + rate_, + nrate, + delta_hz + ); + + double t1 = now(); + + if (t1 - t0 > 0.1) + { + fprintf(stderr, "reduce oops, size %d -> %d, rate %d -> %d, took %.2f\n", + osize, + (int)samples_.size(), + rate_, + nrate, + t1 - t0); + } + + if (0) + { + fprintf(stderr, "%.0f..%.0f, range %.0f, rate %d -> %d, delta hz %.0f, %.6f sec\n", + min_hz_, max_hz_, + max_hz_ - min_hz_, + rate_, nrate, delta_hz, t1 - t0); + } + + if (delta_hz > 0) + { + down_hz_ = delta_hz; // to adjust hz for Python. + min_hz_ -= down_hz_; + max_hz_ -= down_hz_; + + for (int i = 0; i < (int)prevdecs_.size(); i++) + { + prevdecs_[i].hz0 -= delta_hz; + prevdecs_[i].hz1 -= delta_hz; + } + } + + // assert(max_hz_ + 50 < nrate / 2); + // assert(min_hz_ >= 0); + + float ratio = nrate / (float)rate_; + rate_ = nrate; + start_ = round(start_ * ratio); + } + + int block = blocksize(rate_); + + // start_ is the sample number of 0.5 seconds, the nominal start time. + + // make sure there's at least tplus*rate_ samples after the end. + if (start_ + params.tplus * rate_ + 79 * block + block > samples_.size()) + { + int need = start_ + params.tplus * rate_ + 79 * block - samples_.size(); + + // round up to a whole second, to ease fft plan caching. + if ((need % rate_) != 0) { + need += rate_ - (need % rate_); + } + + std::default_random_engine generator; + std::uniform_int_distribution distribution(0, samples_.size() - 1); + auto rnd = std::bind(distribution, generator); + std::vector v(need); + + for (int i = 0; i < need; i++) + { + // v[i] = 0; + v[i] = samples_[rnd()]; + } + + samples_.insert(samples_.end(), v.begin(), v.end()); + } + + int si0 = (start_ - params.tminus * rate_) / block; + + if (si0 < 0) { + si0 = 0; + } + + int si1 = (start_ + params.tplus * rate_) / block; + + // a copy from which to subtract. + nsamples_ = samples_; + int any = 0; + + for (int i = 0; i < (int)prevdecs_.size(); i++) + { + auto d = prevdecs_[i]; + + if (d.hz0 >= min_hz_ && d.hz0 <= max_hz_) + { + // reconstruct correct 79 symbols from LDPC output. + std::vector re79 = recode(d.bits); + + // fine up hz/off again now that we have more samples + float best_hz = (d.hz0 + d.hz1) / 2.0; + float best_off = d.off; // seconds + search_both_known( + samples_, + rate_, + re79, + best_hz, + best_off, + best_hz, + best_off + ); + + // subtract from nsamples_. + subtract(re79, best_hz, best_hz, best_off); + any += 1; + } + } + + if (any) { + samples_ = nsamples_; + } + + for (pass_ = 0; pass_ < npasses; pass_++) + { + double total_remaining = deadline_ - now(); + double remaining = total_remaining / (npasses - pass_); + + if (pass_ == 0) { + remaining *= params.pass0_frac; + } + + double deadline = now() + remaining; + int new_decodes = 0; + samples_ = nsamples_; + + std::vector order; + + // + // search coarsely for Costas blocks. + // in fractions of bins in off and hz. + // + + // just do this once, reuse for every fractional fft_shift + // and down_v7_f() to 200 sps. + std::vector> bins = fftEngine_->one_fft( + samples_, 0, samples_.size()); + + for (int hz_frac_i = 0; hz_frac_i < params.coarse_hz_n; hz_frac_i++) + { + // shift down by hz_frac + float hz_frac = hz_frac_i * (6.25 / params.coarse_hz_n); + std::vector samples1; + + if (hz_frac_i == 0) { + samples1 = samples_; + } else { + samples1 = fft_shift_f(bins, rate_, hz_frac); + } + + for (int off_frac_i = 0; off_frac_i < params.coarse_off_n; off_frac_i++) + { + int off_frac = off_frac_i * (block / params.coarse_off_n); + FFTEngine::ffts_t bins = fftEngine_->ffts(samples1, off_frac, block); + std::vector oo = coarse(bins, si0, si1); + + for (int i = 0; i < (int)oo.size(); i++) + { + oo[i].hz_ += hz_frac; + oo[i].off_ += off_frac; + } + + order.insert(order.end(), oo.begin(), oo.end()); + } + } + + // + // sort strongest-first. + // + std::sort( + order.begin(), + order.end(), + [](const Strength &a, const Strength &b) -> bool { + return a.strength_ > b.strength_; + } + ); + + char already[2000]; // XXX + + for (int i = 0; i < (int)(sizeof(already) / sizeof(already[0])); i++) { + already[i] = 0; + } + + for (int ii = 0; ii < (int)order.size(); ii++) + { + double tt = now(); + + if (ii > 0 && + tt > deadline && + (tt > deadline_ || new_decodes >= params.pass_threshold) && + (pass_ < npasses - 1 || tt > final_deadline_) + ) { + break; + } + + float hz = order[ii].hz_; + + if (already[(int)round(hz / params.already_hz)]) { + continue; + } + + int off = order[ii].off_; + int ret = one_merge(bins, samples_.size(), hz, off); + + if (ret) + { + if (ret == 2) { + new_decodes++; + } + + already[(int)round(hz / params.already_hz)] = 1; + } + } + } // pass +} + +// +// what's the strength of the Costas sync blocks of +// the signal starting at hz and off? +// +float FT4::one_strength(const std::vector &samples200, float hz, int off) +{ + int bin0 = round(hz / 6.25); + int costas[] = {3, 1, 4, 0, 6, 5, 2}; + int starts[] = {0, 36, 72}; + float sig = 0; + float noise = 0; + + for (int which = 0; which < 3; which++) + { + int start = starts[which]; + + for (int si = 0; si < 7; si++) + { + auto fft = fftEngine_->one_fft(samples200, off + (si + start) * 32, 32); + for (int bi = 0; bi < 8; bi++) + { + float x = std::abs(fft[bin0 + bi]); + + if (bi == costas[si]) { + sig += x; + } else { + noise += x; + } + } + } + } + + if (params.strength_how == 0) { + return sig - noise; + } else if (params.strength_how == 1) { + return sig - noise / 7; + } else if (params.strength_how == 2) { + return sig / (noise / 7); + } else if (params.strength_how == 3) { + return sig / (sig + (noise / 7)); + } else if (params.strength_how == 4) { + return sig; + } else if (params.strength_how == 5) { + return sig / (sig + noise); + } else if (params.strength_how == 6) { + return sig / noise; + } else { + return 0; + } +} + +// +// 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 FT4::one_strength_known( + const std::vector &samples, + int rate, + const std::vector &syms, + float hz, + int off +) +{ + int block = blocksize(rate); + // assert(syms.size() == 79); + int bin0 = round(hz / 6.25); + float sig = 0; + float noise = 0; + float sum7 = 0; + std::complex prev = 0; + + for (int si = 0; si < 79; si += params.known_sparse) + { + auto fft = fftEngine_->one_fft(samples, off + si * block, block); + + if (params.known_strength_how == 7) + { + std::complex c = fft[bin0 + syms[si]]; + + if (si > 0) { + sum7 += std::abs(c - prev); + } + + prev = c; + } + else + { + for (int bi = 0; bi < 8; bi++) + { + float x = std::abs(fft[bin0 + bi]); + + if (bi == syms[si]) { + sig += x; + } else { + noise += x; + } + } + } + } + + if (params.known_strength_how == 0) { + return sig - noise; + } else if (params.known_strength_how == 1) { + return sig - noise / 7; + } else if (params.known_strength_how == 2) { + return sig / (noise / 7); + } else if (params.known_strength_how == 3) { + return sig / (sig + (noise / 7)); + } else if (params.known_strength_how == 4) { + return sig; + } else if (params.known_strength_how == 5) { + return sig / (sig + noise); + } else if (params.known_strength_how == 6) { + return sig / noise; + } else if (params.known_strength_how == 7) { + return -sum7; + } else { + return 0; + } +} + +int FT4::search_time_fine( + const std::vector &samples200, + int off0, + int offN, + float hz, + int gran, + float &str +) +{ + if (off0 < 0) { + off0 = 0; + } + + // + // shift in frequency to put hz at 25. + // only shift the samples we need, both for speed, + // and try to always shift down the same number of samples + // to make it easier to cache fftw plans. + // + int len = (offN - off0) + 79 * 32 + 32; + + if (off0 + len > (int)samples200.size()) + { + // len = samples200.size() - off0; + // don't provoke random-length FFTs. + return -1; + } + + std::vector downsamples200 = shift200(samples200, off0, len, hz); + int best_off = -1; + float best_sum = 0.0; + + for (int g = 0; g <= (offN - off0) && g + 79 * 32 <= len; g += gran) + { + float sum = one_strength(downsamples200, 25, g); + + if (sum > best_sum || best_off == -1) + { + best_off = g; + best_sum = sum; + } + } + + str = best_sum; + // assert(best_off >= 0); + return off0 + best_off; +} + +int FT4::search_time_fine_known( + const std::vector> &bins, + int rate, + const std::vector &syms, + int off0, + int offN, + float hz, + int gran, + float &str +) +{ + if (off0 < 0) { + off0 = 0; + } + + // nearest FFT bin center. + float hz0 = round(hz / 6.25) * 6.25; + // move hz to hz0, so it is centered in a symbol-sized bin. + std::vector downsamples = fft_shift_f(bins, rate, hz - hz0); + int best_off = -1; + int block = blocksize(rate); + float best_sum = 0.0; + + for (int g = off0; g <= offN; g += gran) + { + if (g >= 0 && g + 79 * block <= (int)downsamples.size()) + { + float sum = one_strength_known(downsamples, rate, syms, hz0, g); + + if (sum > best_sum || best_off == -1) + { + best_off = g; + best_sum = sum; + } + } + } + + if (best_off < 0) { + return -1; + } + + str = best_sum; + return best_off; +} + +// +// 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 FT4::search_both( + const std::vector &samples200, + float hz0, + int hz_n, + float hz_win, + int off0, + int off_n, + int off_win +) +{ + // assert(hz0 >= 25 - 6.25 / 2 && hz0 <= 25 + 6.25 / 2); + + std::vector strengths; + + float hz_inc = 2 * hz_win / hz_n; + int off_inc = round(2 * off_win / (float)off_n); + + if (off_inc < 1) { + off_inc = 1; + } + + for (float hz = hz0 - hz_win; hz <= hz0 + hz_win + 0.01; hz += hz_inc) + { + float str = 0; + int off = search_time_fine( + samples200, + off0 - off_win, + off0 + off_win, hz, + off_inc, + str + ); + + if (off >= 0) + { + Strength st; + st.hz_ = hz; + st.off_ = off; + st.strength_ = str; + strengths.push_back(st); + } + } + + return strengths; +} + +void FT4::search_both_known( + const std::vector &samples, + int rate, + const std::vector &syms, + float hz0, + float off_secs0, // seconds + float &hz_out, + float &off_out +) +{ + // assert(hz0 >= 0 && hz0 + 50 < rate / 2); + + int off0 = round(off_secs0 * (float)rate); + int off_win = params.third_off_win * blocksize(rate_); + + if (off_win < 1) { + off_win = 1; + } + + int off_inc = trunc((2.0 * off_win) / (params.third_off_n - 1.0)); + + if (off_inc < 1) { + off_inc = 1; + } + + int got_best = 0; + float best_hz = 0; + int best_off = 0; + float best_strength = 0; + std::vector> bins = fftEngine_->one_fft(samples, 0, samples.size()); + float hz_start, hz_inc, hz_end; + + if (params.third_hz_n > 1) + { + hz_inc = (2.0 * params.third_hz_win) / (params.third_hz_n - 1.0); + hz_start = hz0 - params.third_hz_win; + hz_end = hz0 + params.third_hz_win; + } + else + { + hz_inc = 1; + hz_start = hz0; + hz_end = hz0; + } + + for (float hz = hz_start; hz <= hz_end + 0.0001; hz += hz_inc) + { + float strength = 0; + int off = search_time_fine_known( + bins, + rate, + syms, + off0 - off_win, + off0 + off_win, + hz, + off_inc, + strength + ); + + if (off >= 0 && (got_best == 0 || strength > best_strength)) + { + got_best = 1; + best_hz = hz; + best_off = off; + best_strength = strength; + } + } + + if (got_best) + { + hz_out = best_hz; + off_out = best_off / (float)rate; + } +} + +// +// 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 FT4::fft_shift( + const std::vector &samples, + int off, + int len, + int rate, + float hz +) +{ + std::vector> bins; + // horrible hack to avoid repeated FFTs on the same input. + hack_mu_.lock(); + + if ((int)samples.size() == hack_size_ && samples.data() == hack_data_ && + off == hack_off_ && len == hack_len_ && + samples[0] == hack_0_ && samples[1] == hack_1_) + { + bins = hack_bins_; + } + else + { + bins = fftEngine_->one_fft(samples, off, len); + hack_bins_ = bins; + hack_size_ = samples.size(); + hack_off_ = off; + hack_len_ = len; + hack_0_ = samples[0]; + hack_1_ = samples[1]; + hack_data_ = samples.data(); + } + + hack_mu_.unlock(); + return fft_shift_f(bins, rate, hz); +} + +// +// shift down by hz. +// +std::vector FT4::fft_shift_f( + const std::vector> &bins, + int rate, + float hz +) +{ + int nbins = bins.size(); + int len = (nbins - 1) * 2; + float bin_hz = rate / (float)len; + int down = round(hz / bin_hz); + std::vector> bins1(nbins); + + for (int i = 0; i < nbins; i++) + { + int j = i + down; + + if (j >= 0 && j < nbins) { + bins1[i] = bins[j]; + } else { + bins1[i] = 0; + } + } + + std::vector out = fftEngine_->one_ifft(bins1); + return out; +} + +// shift the frequency by a fraction of 6.25, +// to center hz on bin 4 (25 hz). +std::vector FT4::shift200( + const std::vector &samples200, + int off, + int len, + float hz +) +{ + if (std::abs(hz - 25) < 0.001 && off == 0 && len == (int)samples200.size()) { + return samples200; + } else { + return fft_shift(samples200, off, len, 200, hz - 25.0); + } + // return hilbert_shift(samples200, hz - 25.0, hz - 25.0, 200); +} + +// returns a mini-FFT of 79 8-tone symbols. +FFTEngine::ffts_t FT4::extract(const std::vector &samples200, float, int off) +{ + + FFTEngine::ffts_t bins3 = fftEngine_->ffts(samples200, off, 32); + FFTEngine::ffts_t m79(79); + + for (int si = 0; si < 79; si++) + { + m79[si].resize(8); + + if (si < (int)bins3.size()) + { + for (int bi = 0; bi < 8; bi++) + { + auto x = bins3[si][4 + bi]; + m79[si][bi] = x; + } + } + else + { + for (int bi = 0; bi < 8; bi++) { + m79[si][bi] = 0; + } + } + } + + return m79; +} + +// +// m79 is a 79x8 array of complex. +// +FFTEngine::ffts_t FT4::un_gray_code_c(const FFTEngine::ffts_t &m79) +{ + FFTEngine::ffts_t m79a(79); + int map[] = {0, 1, 3, 2, 6, 4, 5, 7}; + + for (int si = 0; si < 79; si++) + { + m79a[si].resize(8); + + for (int bi = 0; bi < 8; bi++) { + m79a[si][map[bi]] = m79[si][bi]; + } + } + + return m79a; +} + +// +// m79 is a 79x8 array of float. +// +std::vector> FT4::un_gray_code_r(const std::vector> &m79) +{ + std::vector> m79a(79); + int map[] = {0, 1, 3, 2, 6, 4, 5, 7}; + + for (int si = 0; si < 79; si++) + { + m79a[si].resize(8); + + for (int bi = 0; bi < 8; bi++) { + m79a[si][map[bi]] = m79[si][bi]; + } + } + + return m79a; +} + +// +// Generic Gray decoding for magnitudes (floats) +// +std::vector> FT4::un_gray_code_r_gen(const std::vector> &mags) +{ + if (mags.size() == 0) { + return mags; + } + + std::vector> magsa(mags.size()); + int nsyms = mags.front().size(); + + for (unsigned int si = 0; si < mags.size(); si++) + { + magsa[si].resize(nsyms); + + for (int bini = 0; bini < nsyms; bini++) + { + int grayi = bini ^ (bini >> 1); + magsa[si][bini] = mags[si][grayi]; + } + } + + return magsa; +} + +// +// normalize levels by windowed median. +// this helps, but why? +// +std::vector> FT4::convert_to_snr(const std::vector> &m79) +{ + if (params.snr_how < 0 || params.snr_win < 0) { + return m79; + } + + // + // for each symbol time, what's its "noise" level? + // + std::vector mm(79); + + for (int si = 0; si < 79; si++) + { + std::vector v(8); + float sum = 0.0; + + for (int bi = 0; bi < 8; bi++) + { + float x = m79[si][bi]; + v[bi] = x; + sum += x; + } + + if (params.snr_how != 1) { + std::sort(v.begin(), v.end()); + } + + if (params.snr_how == 0) { + // median + mm[si] = (v[3] + v[4]) / 2; + } else if (params.snr_how == 1) { + mm[si] = sum / 8; + } else if (params.snr_how == 2) { + // all but strongest tone. + mm[si] = (v[0] + v[1] + v[2] + v[3] + v[4] + v[5] + v[6]) / 7; + } else if (params.snr_how == 3) { + mm[si] = v[0]; // weakest tone + } else if (params.snr_how == 4) { + mm[si] = v[7]; // strongest tone + } else if (params.snr_how == 5) { + mm[si] = v[6]; // second-strongest tone + } else { + mm[si] = 1.0; + } + } + + // we're going to take a windowed average. + std::vector winwin; + + if (params.snr_win > 0) { + winwin = blackman(2 * params.snr_win + 1); + } else { + winwin.push_back(1.0); + } + + std::vector> n79(79); + + for (int si = 0; si < 79; si++) + { + float sum = 0; + + for (int dd = si - params.snr_win; dd <= si + params.snr_win; dd++) + { + int wi = dd - (si - params.snr_win); + + if (dd >= 0 && dd < 79) { + sum += mm[dd] * winwin[wi]; + } else if (dd < 0) { + sum += mm[0] * winwin[wi]; + } else { + sum += mm[78] * winwin[wi]; + } + } + + n79[si].resize(8); + + for (int bi = 0; bi < 8; bi++) { + n79[si][bi] = m79[si][bi] / sum; + } + } + + return n79; +} + +// +// normalize levels by windowed median. +// this helps, but why? +// +std::vector>> FT4::c_convert_to_snr( + const std::vector>> &m79 +) +{ + if (params.snr_how < 0 || params.snr_win < 0) { + return m79; + } + + // + // for each symbol time, what's its "noise" level? + // + std::vector mm(79); + + for (int si = 0; si < 79; si++) + { + std::vector v(8); + float sum = 0.0; + + for (int bi = 0; bi < 8; bi++) + { + float x = std::abs(m79[si][bi]); + v[bi] = x; + sum += x; + } + + if (params.snr_how != 1) { + std::sort(v.begin(), v.end()); + } + + if (params.snr_how == 0) { + // median + mm[si] = (v[3] + v[4]) / 2; + } else if (params.snr_how == 1) { + mm[si] = sum / 8; + } else if (params.snr_how == 2) { + // all but strongest tone. + mm[si] = (v[0] + v[1] + v[2] + v[3] + v[4] + v[5] + v[6]) / 7; + } else if (params.snr_how == 3) { + mm[si] = v[0]; // weakest tone + } else if (params.snr_how == 4) { + mm[si] = v[7]; // strongest tone + } else if (params.snr_how == 5) { + mm[si] = v[6]; // second-strongest tone + } else { + mm[si] = 1.0; + } + } + + // we're going to take a windowed average. + std::vector winwin; + + if (params.snr_win > 0) { + winwin = blackman(2 * params.snr_win + 1); + } else { + winwin.push_back(1.0); + } + + std::vector>> n79(79); + + for (int si = 0; si < 79; si++) + { + float sum = 0; + + for (int dd = si - params.snr_win; dd <= si + params.snr_win; dd++) + { + int wi = dd - (si - params.snr_win); + + if (dd >= 0 && dd < 79) { + sum += mm[dd] * winwin[wi]; + } else if (dd < 0) { + sum += mm[0] * winwin[wi]; + } else { + sum += mm[78] * winwin[wi]; + } + } + + n79[si].resize(8); + + for (int bi = 0; bi < 8; bi++) { + n79[si][bi] = m79[si][bi] / sum; + } + } + + return n79; +} + +std::vector> FT4::convert_to_snr_gen(const FT4Params& params, int nbSymbolBits, const std::vector> &mags) +{ + if (params.snr_how < 0 || params.snr_win < 0) { + return mags; + } + + // + // for each symbol time, what's its "noise" level? + // + std::vector mm(mags.size()); + int nbSymbols = 1< v(nbSymbols); + float sum = 0.0; + + for (int bini = 0; bini < nbSymbols; bini++) + { + float x = mags[si][bini]; + v[bini] = x; + sum += x; + } + + if (params.snr_how != 1) { + std::sort(v.begin(), v.end()); + } + + int mid = nbSymbols / 2; + + if (params.snr_how == 0) { + // median + mm[si] = (v[mid-1] + v[mid]) / 2; + } else if (params.snr_how == 1) { + mm[si] = sum / nbSymbols; + } else if (params.snr_how == 2) { + // all but strongest tone. + mm[si] = std::accumulate(v.begin(), v.end() - 1, 0.0f) / (v.size() - 1); + } else if (params.snr_how == 3) { + mm[si] = v.front(); // weakest tone + } else if (params.snr_how == 4) { + mm[si] = v.back(); // strongest tone + } else if (params.snr_how == 5) { + mm[si] = v[v.size()-2]; // second-strongest tone + } else { + mm[si] = 1.0; + } + } + + // we're going to take a windowed average. + std::vector winwin; + + if (params.snr_win > 0) { + winwin = blackman(2 * params.snr_win + 1); + } else { + winwin.push_back(1.0); + } + + std::vector> snr(mags.size()); + + for (int si = 0; si < (int) mags.size(); si++) + { + float sum = 0; + + for (int dd = si - params.snr_win; dd <= si + params.snr_win; dd++) + { + int wi = dd - (si - params.snr_win); + + if (dd >= 0 && dd < (int) mags.size()) { + sum += mm[dd] * winwin[wi]; + } else if (dd < 0) { + sum += mm[0] * winwin[wi]; + } else { + sum += mm[mags.size()-1] * winwin[wi]; + } + } + + snr[si].resize(nbSymbols); + + for (int bi = 0; bi < nbSymbols; bi++) { + snr[si][bi] = mags[si][bi] / sum; + } + } + + return snr; +} + +// +// statistics to decide soft probabilities, +// to drive LDPC decoder. +// distribution of strongest tones, and +// distribution of noise. +// +void FT4::make_stats( + const std::vector> &m79, + Stats &bests, + Stats &all +) +{ + int costas[] = {3, 1, 4, 0, 6, 5, 2}; + + for (int si = 0; si < 79; si++) + { + if (si < 7 || (si >= 36 && si < 36 + 7) || si >= 72) + { + // Costas. + int ci; + + if (si >= 72) { + ci = si - 72; + } else if (si >= 36) { + ci = si - 36; + } else { + ci = si; + } + + for (int bi = 0; bi < 8; bi++) + { + float x = m79[si][bi]; + all.add(x); + + if (bi == costas[ci]) { + bests.add(x); + } + } + } + else + { + float mx = 0; + + for (int bi = 0; bi < 8; bi++) + { + float x = m79[si][bi]; + + if (x > mx) { + mx = x; + } + + all.add(x); + } + + bests.add(mx); + } + } +} + +// +// generalized version of the above for any number of symbols and no Costas +// used by FT-chirp decoder +// +void FT4::make_stats_gen( + const std::vector> &mags, + int nbSymbolBits, + Stats &bests, + Stats &all +) +{ + int nbBins = 1< mx) { + mx = x; + } + + all.add(x); + } + + bests.add(mx); + } +} + +// +// 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> FT4::soft_c2m(const FFTEngine::ffts_t &c79) +{ + std::vector> m79(79); + std::vector raw_phases(79); // of strongest tone in each symbol time + + for (int si = 0; si < 79; si++) + { + m79[si].resize(8); + int mxi = -1; + float mx; + float mx_phase; + + for (int bi = 0; bi < 8; bi++) + { + float x = std::abs(c79[si][bi]); + m79[si][bi] = x; + + if (mxi < 0 || x > mx) + { + mxi = bi; + mx = x; + mx_phase = std::arg(c79[si][bi]); // -pi .. pi + } + } + + raw_phases[si] = mx_phase; + } + + if (params.soft_phase_win <= 0) { + return m79; + } + + // phase around each symbol. + std::vector phases(79); + + // for each symbol time, median of nearby phases + for (int si = 0; si < 79; si++) + { + std::vector v; + + for (int si1 = si - params.soft_phase_win; si1 <= si + params.soft_phase_win; si1++) + { + if (si1 >= 0 && si1 < 79) + { + float x = raw_phases[si1]; + v.push_back(x); + } + } + + // choose the phase that has the lowest total distance to other + // phases. like median but avoids -pi..pi wrap-around. + int n = v.size(); + int best = -1; + float best_score = 0; + for (int i = 0; i < n; i++) + + { + float score = 0; + + for (int j = 0; j < n; j++) + { + if (i == j) { + continue; + } + + float d = fabs(v[i] - v[j]); + + if (d > M_PI) { + d = 2 * M_PI - d; + } + + score += d; + } + + if (best == -1 || score < best_score) + { + best = i; + best_score = score; + } + } + + phases[si] = v[best]; + } + + // project each tone against the median phase around that symbol time. + for (int si = 0; si < 79; si++) + { + for (int bi = 0; bi < 8; bi++) + { + float mag = std::abs(c79[si][bi]); + float angle = std::arg(c79[si][bi]); + float d = angle - phases[si]; + float factor = 0.1; + + if (d < M_PI / 2 && d > -M_PI / 2) { + factor = cos(d); + } + + m79[si][bi] = factor * mag; + } + } + + return m79; +} + +// +// 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 FT4::bayes( + FT4Params& params, + float best_zero, + float best_one, + int lli, + Stats &bests, + Stats &all +) +{ + float maxlog = 4.97; + float ll = 0; + float pzero = 0.5; + float pone = 0.5; + + if (params.use_apriori) + { + pzero = 1.0 - apriori174[lli]; + pone = apriori174[lli]; + } + + // + // Bayes combining rule normalization from: + // http://cs.wellesley.edu/~anderson/writing/naive-bayes.pdf + // + // a = P(zero)P(e0|zero)P(e1|zero) + // b = P(one)P(e0|one)P(e1|one) + // p = a / (a + b) + // + // also see Mark Owen's book Practical Signal Processing, + // Chapter 6. + // + + // zero + float a = pzero * bests.problt(best_zero) * (1.0 - all.problt(best_one)); + // printf("FT8::bayes: a: %f bp: %f ap: %f \n", a, bests.problt(best_zero), all.problt(best_one)); + + if (params.bayes_how == 1) { + a *= all.problt(all.mean() + (best_zero - best_one)); + } + + // one + float b = pone * bests.problt(best_one) * (1.0 - all.problt(best_zero)); + // printf("FT8::bayes: b: %f bp: %f ap: %f \n", b, bests.problt(best_one), all.problt(best_zero)); + + if (params.bayes_how == 1) { + b *= all.problt(all.mean() + (best_one - best_zero)); + } + + float p; + + if (a + b == 0) { + p = 0.5; + } else { + p = a / (a + b); + } + + // printf("FT8::bayes: all.mean: %f a: %f b: %f p: %f\n", all.mean(), a, b, p); + + if (1 - p == 0.0) { + ll = maxlog; + } else { + ll = log(p / (1 - p)); + } + + if (ll > maxlog) { + ll = maxlog; + } + + if (ll < -maxlog) { + ll = -maxlog; + } + + return ll; +} + +// +// c79 is 79x8 complex tones, before un-gray-coding. +// +void FT4::soft_decode(const FFTEngine::ffts_t &c79, float ll174[]) +{ + std::vector> m79(79); + // m79 = absolute values of c79. + // still pre-un-gray-coding so we know which + // are the correct Costas tones. + m79 = soft_c2m(c79); + m79 = convert_to_snr(m79); + // statistics to decide soft probabilities. + // distribution of strongest tones, and + // distribution of noise. + Stats bests(params.problt_how_sig, params.log_tail, params.log_rate); + Stats all(params.problt_how_noise, params.log_tail, params.log_rate); + make_stats(m79, bests, all); + m79 = un_gray_code_r(m79); + int lli = 0; + + // tone numbers that make second index bit zero or one. + int zeroi[4][3]; + int onei[4][3]; + + for (int biti = 0; biti < 3; biti++) + { + if (biti == 0) + { + // high bit + zeroi[0][0] = 0; + zeroi[1][0] = 1; + zeroi[2][0] = 2; + zeroi[3][0] = 3; + onei[0][0] = 4; + onei[1][0] = 5; + onei[2][0] = 6; + onei[3][0] = 7; + } + + if (biti == 1) + { + // middle bit + zeroi[0][1] = 0; + zeroi[1][1] = 1; + zeroi[2][1] = 4; + zeroi[3][1] = 5; + onei[0][1] = 2; + onei[1][1] = 3; + onei[2][1] = 6; + onei[3][1] = 7; + } + + if (biti == 2) + { + // low bit + zeroi[0][2] = 0; + zeroi[1][2] = 2; + zeroi[2][2] = 4; + zeroi[3][2] = 6; + onei[0][2] = 1; + onei[1][2] = 3; + onei[2][2] = 5; + onei[3][2] = 7; + } + } + + for (int i79 = 0; i79 < 79; i79++) + { + if (i79 < 7 || (i79 >= 36 && i79 < 36 + 7) || i79 >= 72) { + // Costas, skip + continue; + } + + // for each of the three bits, look at the strongest tone + // that would make it a zero, and the strongest tone that + // would make it a one. use Bayes to decide which is more + // likely, comparing each against the distribution of noise + // and the distribution of strongest tones. + // most-significant-bit first. + + for (int biti = 0; biti < 3; biti++) + { + // strongest tone that would make this bit be zero. + int got_best_zero = 0; + float best_zero = 0; + + for (int i = 0; i < 4; i++) + { + float x = m79[i79][zeroi[i][biti]]; + + if (got_best_zero == 0 || x > best_zero) + { + got_best_zero = 1; + best_zero = x; + } + } + + // strongest tone that would make this bit be one. + int got_best_one = 0; + float best_one = 0; + + for (int i = 0; i < 4; i++) + { + float x = m79[i79][onei[i][biti]]; + if (got_best_one == 0 || x > best_one) + { + got_best_one = 1; + best_one = x; + } + } + + float ll = bayes(params, best_zero, best_one, lli, bests, all); + ll174[lli++] = ll; + } + } + // assert(lli == 174); +} + +// +// mags is the vector of 2^nbSymbolBits vector of magnitudes at each symbol time +// ll174 is the resulting 174 soft bits of payload +// used in FT-chirp modulation scheme - generalized to any number of symbol bits +// +void FT4::soft_decode_mags(FT4Params& params, const std::vector>& mags_, int nbSymbolBits, float ll174[]) +{ + if ((nbSymbolBits > 16) || (nbSymbolBits < 1)) { + return; + } + + std::vector> mags = convert_to_snr_gen(params, nbSymbolBits, mags_); + // statistics to decide soft probabilities. + // distribution of strongest tones, and + // distribution of noise. + Stats bests(params.problt_how_sig, params.log_tail, params.log_rate); + Stats all(params.problt_how_noise, params.log_tail, params.log_rate); + make_stats_gen(mags, nbSymbolBits, bests, all); + mags = un_gray_code_r_gen(mags); + int lli = 0; + int zoX = 1<<(nbSymbolBits-1); + int zoY = nbSymbolBits; + std::vector zeroi(zoX*zoY); + std::vector onei(zoX*zoY); + + for (int biti = 0; biti < nbSymbolBits; biti++) + { + int i = biti * zoX; + set_ones_zeroes(&onei[i], &zeroi[i], nbSymbolBits, biti); + } + + for (unsigned int si = 0; si < mags.size(); si++) + { + // for each of the symbol bits, look at the strongest tone + // that would make it a zero, and the strongest tone that + // would make it a one. use Bayes to decide which is more + // likely, comparing each against the distribution of noise + // and the distribution of strongest tones. + // most-significant-bit first. + for (int biti = nbSymbolBits - 1; biti >= 0; biti--) + { + // strongest tone that would make this bit be zero. + int got_best_zero = 0; + float best_zero = 0; + + for (int i = 0; i < 1<<(nbSymbolBits-1); i++) + { + float x = mags[si][zeroi[i+biti*zoX]]; + // printf("FT8::soft_decode_mags:: biti: %d i: %d zeroi: %d x: %f best_zero: %f\n", biti, i, zeroi[i+biti*zoX], x, best_zero); + + if (got_best_zero == 0 || x > best_zero) + { + got_best_zero = 1; + best_zero = x; + } + } + + // strongest tone that would make this bit be one. + int got_best_one = 0; + float best_one = 0; + + for (int i = 0; i < 1<<(nbSymbolBits-1); i++) + { + float x = mags[si][onei[i+biti*zoX]]; + // printf("FT8::soft_decode_mags:: biti: %d i: %d onei: %d x: %f best_one: %f\n", biti, i, onei[i+biti*zoX], x, best_one); + + if (got_best_one == 0 || x > best_one) + { + got_best_one = 1; + best_one = x; + } + } + + // printf("FT8::soft_decode_mags: biti: %d best_zero: %f best_one: %f\n", biti, best_zero, best_one); + + float ll = bayes(params, best_zero, best_one, lli, bests, all); + ll174[lli++] = ll; + } + } +} + +// +// c79 is 79x8 complex tones, before un-gray-coding. +// +void FT4::c_soft_decode(const FFTEngine::ffts_t &c79x, float ll174[]) +{ + FFTEngine::ffts_t c79 = c_convert_to_snr(c79x); + int costas[] = {3, 1, 4, 0, 6, 5, 2}; + std::complex maxes[79]; + + for (int i = 0; i < 79; i++) + { + std::complex m; + + if (i < 7) + { + // Costas. + m = c79[i][costas[i]]; + } + else if (i >= 36 && i < 36 + 7) + { + // Costas. + m = c79[i][costas[i - 36]]; + } + else if (i >= 72) + { + // Costas. + m = c79[i][costas[i - 72]]; + } + else + { + int got = 0; + + for (int j = 0; j < 8; j++) + { + if (got == 0 || std::abs(c79[i][j]) > std::abs(m)) + { + got = 1; + m = c79[i][j]; + } + } + } + + maxes[i] = m; + } + + std::vector> m79(79); + + for (int i = 0; i < 79; i++) + { + m79[i].resize(8); + + for (int j = 0; j < 8; j++) + { + std::complex c = c79[i][j]; + int n = 0; + float sum = 0; + + for (int k = i - params.c_soft_win; k <= i + params.c_soft_win; k++) + { + if (k < 0 || k >= 79) { + continue; + } + + if (k == i) + { + sum -= params.c_soft_weight * std::abs(c); + } + else + { + // we're expecting all genuine tones to have + // about the same phase and magnitude. + // so set m79[i][j] to the distance from the + // phase/magnitude predicted by surrounding + // genuine-looking tones. + std::complex c1 = maxes[k]; + std::complex d = c1 - c; + sum += std::abs(d); + } + + n += 1; + } + + m79[i][j] = 0 - (sum / n); + } + } + + // statistics to decide soft probabilities. + // distribution of strongest tones, and + // distribution of noise. + Stats bests(params.problt_how_sig, params.log_tail, params.log_rate); + Stats all(params.problt_how_noise, params.log_tail, params.log_rate); + make_stats(m79, bests, all); + m79 = un_gray_code_r(m79); + int lli = 0; + + // tone numbers that make second index bit zero or one. + int zeroi[4][3]; + int onei[4][3]; + + for (int biti = 0; biti < 3; biti++) + { + if (biti == 0) + { + // high bit + zeroi[0][0] = 0; + zeroi[1][0] = 1; + zeroi[2][0] = 2; + zeroi[3][0] = 3; + onei[0][0] = 4; + onei[1][0] = 5; + onei[2][0] = 6; + onei[3][0] = 7; + } + + if (biti == 1) + { + // middle bit + zeroi[0][1] = 0; + zeroi[1][1] = 1; + zeroi[2][1] = 4; + zeroi[3][1] = 5; + onei[0][1] = 2; + onei[1][1] = 3; + onei[2][1] = 6; + onei[3][1] = 7; + } + + if (biti == 2) + { + // low bit + zeroi[0][2] = 0; + zeroi[1][2] = 2; + zeroi[2][2] = 4; + zeroi[3][2] = 6; + onei[0][2] = 1; + onei[1][2] = 3; + onei[2][2] = 5; + onei[3][2] = 7; + } + } + + for (int i79 = 0; i79 < 79; i79++) + { + if (i79 < 7 || (i79 >= 36 && i79 < 36 + 7) || i79 >= 72) { + // Costas, skip + continue; + } + + // for each of the three bits, look at the strongest tone + // that would make it a zero, and the strongest tone that + // would make it a one. use Bayes to decide which is more + // likely, comparing each against the distribution of noise + // and the distribution of strongest tones. + // most-significant-bit first. + + for (int biti = 0; biti < 3; biti++) + { + // strongest tone that would make this bit be zero. + int got_best_zero = 0; + float best_zero = 0; + + for (int i = 0; i < 4; i++) + { + float x = m79[i79][zeroi[i][biti]]; + + if (got_best_zero == 0 || x > best_zero) + { + got_best_zero = 1; + best_zero = x; + } + } + + // strongest tone that would make this bit be one. + int got_best_one = 0; + float best_one = 0; + + for (int i = 0; i < 4; i++) + { + float x = m79[i79][onei[i][biti]]; + + if (got_best_one == 0 || x > best_one) + { + got_best_one = 1; + best_one = x; + } + } + + float ll = bayes(params, best_zero, best_one, lli, bests, all); + ll174[lli++] = ll; + } + } + // assert(lli == 174); +} + +// +// set ones and zero symbol indexes. Bit index is LSB +// +void FT4::set_ones_zeroes(int ones[], int zeroes[], int nbBits, int bitIndex) +{ + int nbIndexes = 1 << (nbBits - 1); + + if (bitIndex == 0) + { + for (int i = 0; i < nbIndexes; i++) + { + zeroes[i] = i<<1; + ones[i] = zeroes[i] | 1; + } + } + else if (bitIndex == nbBits - 1) + { + for (int i = 0; i < nbIndexes; i++) + { + zeroes[i] = i; + ones[i] = (1<<(nbBits-1)) | zeroes[i]; + } + } + else + { + int mask = (1< 0 for zero, +// scaled by str. +// +std::vector FT4::extract_bits(const std::vector &syms, const std::vector str) +{ + // assert(syms.size() == 79); + // assert(str.size() == 79); + + std::vector bits; + + for (int si = 0; si < 79; si++) + { + if (si < 7 || (si >= 36 && si < 36 + 7) || si >= 72) + { + // costas -- skip + } + else + { + bits.push_back((syms[si] & 4) == 0 ? str[si] : -str[si]); + bits.push_back((syms[si] & 2) == 0 ? str[si] : -str[si]); + bits.push_back((syms[si] & 1) == 0 ? str[si] : -str[si]); + } + } + + return bits; +} + +// decode successive pairs of symbols. exploits the likelihood +// that they have the same phase, by summing the complex +// correlations for each possible pair and using the max. +void FT4::soft_decode_pairs( + const FFTEngine::ffts_t &m79x, + float ll174[] +) +{ + FFTEngine::ffts_t m79 = c_convert_to_snr(m79x); + + struct BitInfo + { + float zero; // strongest correlation that makes it zero + float one; // and one + }; + + std::vector bitinfo(79 * 3); + + for (int i = 0; i < (int)bitinfo.size(); i++) + { + bitinfo[i].zero = 0; + bitinfo[i].one = 0; + } + + Stats all(params.problt_how_noise, params.log_tail, params.log_rate); + Stats bests(params.problt_how_sig, params.log_tail, params.log_rate); + int map[] = {0, 1, 3, 2, 6, 4, 5, 7}; // un-gray-code + + for (int si = 0; si < 79; si += 2) + { + float mx = 0; + float corrs[8 * 8]; + + for (int s1 = 0; s1 < 8; s1++) + { + for (int s2 = 0; s2 < 8; s2++) + { + // sum up the correlations. + std::complex csum = m79[si][s1]; + + if (si + 1 < 79) { + csum += m79[si + 1][s2]; + } + + float x = std::abs(csum); + corrs[s1 * 8 + s2] = x; + + if (x > mx) { + mx = x; + } + + all.add(x); + // first symbol + int i = map[s1]; + + for (int bit = 0; bit < 3; bit++) + { + int bitind = (si + 0) * 3 + (2 - bit); + + if ((i & (1 << bit))) + { + // symbol i would make this bit a one. + if (x > bitinfo[bitind].one) { + bitinfo[bitind].one = x; + } + } + else + { + // symbol i would make this bit a zero. + if (x > bitinfo[bitind].zero) { + bitinfo[bitind].zero = x; + } + } + } + + // second symbol + if (si + 1 < 79) + { + i = map[s2]; + + for (int bit = 0; bit < 3; bit++) + { + int bitind = (si + 1) * 3 + (2 - bit); + + if ((i & (1 << bit))) + { + // symbol i would make this bit a one. + if (x > bitinfo[bitind].one) { + bitinfo[bitind].one = x; + } + } + else + { + // symbol i would make this bit a zero. + if (x > bitinfo[bitind].zero) { + bitinfo[bitind].zero = x; + } + } + } + } + } + } + + if (si == 0 || si == 36 || si == 72) { + bests.add(corrs[3 * 8 + 1]); + } else if (si == 2 || si == 38 || si == 74) { + bests.add(corrs[4 * 8 + 0]); + } else if (si == 4 || si == 40 || si == 76) { + bests.add(corrs[6 * 8 + 5]); + } else { + bests.add(mx); + } + } + + int lli = 0; + + for (int si = 0; si < 79; si++) + { + if (si < 7 || (si >= 36 && si < 36 + 7) || si >= 72) { + // costas + continue; + } + + for (int i = 0; i < 3; i++) + { + float best_zero = bitinfo[si * 3 + i].zero; + float best_one = bitinfo[si * 3 + i].one; + // ll174[lli++] = best_zero > best_one ? 4.99 : -4.99; + float ll = bayes(params, best_zero, best_one, lli, bests, all); + ll174[lli++] = ll; + } + } + // assert(lli == 174); +} + +void FT4::soft_decode_triples( + const FFTEngine::ffts_t &m79x, + float ll174[] +) +{ + FFTEngine::ffts_t m79 = c_convert_to_snr(m79x); + + struct BitInfo + { + float zero; // strongest correlation that makes it zero + float one; // and one + }; + + std::vector bitinfo(79 * 3); + + for (int i = 0; i < (int)bitinfo.size(); i++) + { + bitinfo[i].zero = 0; + bitinfo[i].one = 0; + } + + Stats all(params.problt_how_noise, params.log_tail, params.log_rate); + Stats bests(params.problt_how_sig, params.log_tail, params.log_rate); + + int map[] = {0, 1, 3, 2, 6, 4, 5, 7}; // un-gray-code + + for (int si = 0; si < 79; si += 3) + { + float mx = 0; + float corrs[8 * 8 * 8]; + + for (int s1 = 0; s1 < 8; s1++) + { + for (int s2 = 0; s2 < 8; s2++) + { + for (int s3 = 0; s3 < 8; s3++) + { + std::complex csum = m79[si][s1]; + + if (si + 1 < 79) { + csum += m79[si + 1][s2]; + } + + if (si + 2 < 79) { + csum += m79[si + 2][s3]; + } + + float x = std::abs(csum); + corrs[s1 * 64 + s2 * 8 + s3] = x; + + if (x > mx) { + mx = x; + } + + all.add(x); + // first symbol + int i = map[s1]; + + for (int bit = 0; bit < 3; bit++) + { + int bitind = (si + 0) * 3 + (2 - bit); + + if ((i & (1 << bit))) + { + // symbol i would make this bit a one. + if (x > bitinfo[bitind].one) { + bitinfo[bitind].one = x; + } + } + else + { + // symbol i would make this bit a zero. + if (x > bitinfo[bitind].zero) { + bitinfo[bitind].zero = x; + } + } + } + + // second symbol + if (si + 1 < 79) + { + i = map[s2]; + + for (int bit = 0; bit < 3; bit++) + { + int bitind = (si + 1) * 3 + (2 - bit); + + if ((i & (1 << bit))) + { + // symbol i would make this bit a one. + if (x > bitinfo[bitind].one) { + bitinfo[bitind].one = x; + } + } + else + { + // symbol i would make this bit a zero. + if (x > bitinfo[bitind].zero) { + bitinfo[bitind].zero = x; + } + } + } + } + + // third symbol + if (si + 2 < 79) + { + i = map[s3]; + + for (int bit = 0; bit < 3; bit++) + { + int bitind = (si + 2) * 3 + (2 - bit); + + if ((i & (1 << bit))) + { + // symbol i would make this bit a one. + if (x > bitinfo[bitind].one) { + bitinfo[bitind].one = x; + } + } + else + { + // symbol i would make this bit a zero. + if (x > bitinfo[bitind].zero) { + bitinfo[bitind].zero = x; + } + } + } + } + } + } + } + + // costas: 3, 1, 4, 0, 6, 5, 2 + if (si == 0 || si == 36 || si == 72) { + bests.add(corrs[3 * 64 + 1 * 8 + 4]); + } else if (si == 3 || si == 39 || si == 75) { + bests.add(corrs[0 * 64 + 6 * 8 + 5]); + } else { + bests.add(mx); + } + } + + int lli = 0; + + for (int si = 0; si < 79; si++) + { + if (si < 7 || (si >= 36 && si < 36 + 7) || si >= 72) { + // costas + continue; + } + + for (int i = 0; i < 3; i++) + { + float best_zero = bitinfo[si * 3 + i].zero; + float best_one = bitinfo[si * 3 + i].one; + // ll174[lli++] = best_zero > best_one ? 4.99 : -4.99; + float ll = bayes(params, best_zero, best_one, lli, bests, all); + ll174[lli++] = ll; + } + } + // assert(lli == 174); +} + +// +// given log likelihood for each bit, try LDPC and OSD decoders. +// on success, puts corrected 174 bits into a174[]. +// +int FT4::decode(const float ll174[], int a174[], FT4Params& _params, int use_osd, std::string &comment) +{ + int plain[174]; // will be 0/1 bits. + int ldpc_ok = 0; // 83 will mean success. + LDPC::ldpc_decode((float *)ll174, _params.ldpc_iters, plain, &ldpc_ok); + int ok_thresh = 83; // 83 is perfect + + if (ldpc_ok >= ok_thresh) + { + // plain[] is 91 systematic data bits, 83 parity bits. + for (int i = 0; i < 174; i++) { + a174[i] = plain[i]; + } + + if (OSD::check_crc(a174)) { + // success! + return 1; + } else { + comment = "CRC fail"; + } + } + else + { + comment = "LDPC fail"; + } + + if (use_osd && _params.osd_depth >= 0 && ldpc_ok >= _params.osd_ldpc_thresh) + { + int oplain[91]; + int got_depth = -1; + int osd_ok = OSD::osd_decode((float *)ll174, _params.osd_depth, oplain, &got_depth); + + if (osd_ok) + { + // reconstruct all 174. + comment += "OSD-" + std::to_string(got_depth) + "-" + std::to_string(ldpc_ok); + OSD::ldpc_encode(oplain, a174); + return 1; + } + else + { + comment = "OSD fail"; + } + } + + return 0; +} + +// +// encode a 77 bit message into a 174 bit payload +// adds the 14 bit CRC to obtain 91 bits +// apply (174, 91) generator mastrix to obtain the 83 parity bits +// append the 83 bits to the 91 bits message + crc to obbain the 174 bit payload +// +void FT4::encode(int a174[], const int s77[]) +{ + int a91[91]; // msg + CRC + std::fill(a91, a91 + 91, 0); + std::copy(s77, s77+77, a91); // copy msg + LDPC::ft8_crc(a91, 82, &a91[77]); // append CRC - to match OSD::check_crc + std::copy(a91, a91+91, a174); // copy msg + CRC + int sum, n, ni, b; + + for (int i=0; i<83; i++) + { + sum = 0; + + for (int j=0; j<91; j++) + { + n = j/8; // byte index in the generator matrix + ni = j%8; // bit index in the generator matrix byte LSB first + b = (Arrays::Gm[i][n] >> (7-ni)) & 1; // bit in the generator matrix + sum += a91[j] * b; + } + + a174[91+i] = sum % 2; // sum modulo 2 + } +} + +// +// 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> FT4::fbandpass( + const std::vector> &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 +) +{ + // assert(low_outer <= low_inner); + // assert(low_inner <= high_inner); + // assert(high_inner <= high_outer); + + int nbins = bins0.size(); + std::vector> bins1(nbins); + + for (int i = 0; i < nbins; i++) + { + float ihz = i * bin_hz; + // cos(x)+flat+cos(x) taper + float factor; + if (ihz <= low_outer || ihz >= high_outer) + { + factor = 0; + } + else if (ihz >= low_outer && ihz < low_inner) + { + // rising shoulder +#if 1 + factor = (ihz - low_outer) / (low_inner - low_outer); // 0 .. 1 +#else + float theta = (ihz - low_outer) / (low_inner - low_outer); // 0 .. 1 + theta -= 1; // -1 .. 0 + theta *= 3.14159; // -pi .. 0 + factor = cos(theta); // -1 .. 1 + factor = (factor + 1) / 2; // 0 .. 1 +#endif + } + else if (ihz > high_inner && ihz <= high_outer) + { + // falling shoulder +#if 1 + factor = (high_outer - ihz) / (high_outer - high_inner); // 1 .. 0 +#else + float theta = (high_outer - ihz) / (high_outer - high_inner); // 1 .. 0 + theta = 1.0 - theta; // 0 .. 1 + theta *= 3.14159; // 0 .. pi + factor = cos(theta); // 1 .. -1 + factor = (factor + 1) / 2; // 1 .. 0 +#endif + } + else + { + factor = 1.0; + } + + bins1[i] = bins0[i] * factor; + } + + return bins1; +} + +// +// 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 FT4::down_v7(const std::vector &samples, float hz) +{ + int len = samples.size(); + std::vector> bins = fftEngine_->one_fft(samples, 0, len); + + return down_v7_f(bins, len, hz); +} + +std::vector FT4::down_v7_f(const std::vector> &bins, int len, float hz) +{ + int nbins = bins.size(); + float bin_hz = rate_ / (float)len; + int down = round((hz - 25) / bin_hz); + std::vector> bins1(nbins); + + for (int i = 0; i < nbins; i++) + { + int j = i + down; + + if (j >= 0 && j < nbins) { + bins1[i] = bins[j]; + } else { + bins1[i] = 0; + } + } + + // now filter to fit in 200 samples/second. + + float low_inner = 25.0 - params.shoulder200_extra; + float low_outer = low_inner - params.shoulder200; + + if (low_outer < 0) { + low_outer = 0; + } + + float high_inner = 75 - 6.25 + params.shoulder200_extra; + float high_outer = high_inner + params.shoulder200; + + if (high_outer > 100) { + high_outer = 100; + } + + bins1 = fbandpass( + bins1, + bin_hz, + low_outer, + low_inner, + high_inner, + high_outer + ); + + // convert back to time domain and down-sample to 200 samples/second. + int blen = round(len * (200.0 / rate_)); + std::vector> bbins(blen / 2 + 1); + + for (int i = 0; i < (int)bbins.size(); i++) { + bbins[i] = bins1[i]; + } + + std::vector out = fftEngine_->one_ifft(bbins); + return out; +} + +// + // 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 FT4::one_merge(const std::vector> &bins, int len, float hz, int off) +{ + // + // set up to search for best frequency and time offset. + // + + // + // move down to 25 hz and re-sample to 200 samples/second, + // i.e. 32 samples/symbol. + // + std::vector samples200 = down_v7_f(bins, len, hz); + int off200 = round((off / (float)rate_) * 200.0); + int ret = one_iter(samples200, off200, hz); + return ret; +} + +// 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 FT4::one_iter(const std::vector &samples200, int best_off, float hz_for_cb) +{ + if (params.do_second) + { + std::vector strengths = search_both( + samples200, + 25, + params.second_hz_n, + params.second_hz_win, + best_off, + params.second_off_n, + params.second_off_win * 32 + ); + // + // sort strongest-first. + // + std::sort( + strengths.begin(), + strengths.end(), + [](const Strength &a, const Strength &b) -> bool { + return a.strength_ > b.strength_; + } + ); + + for (int i = 0; i < (int)strengths.size() && i < params.second_count; i++) + { + float hz = strengths[i].hz_; + int off = strengths[i].off_; + int ret = one_iter1(samples200, off, hz, hz_for_cb, hz_for_cb); + + if (ret > 0) { + return ret; + } + } + } + else + { + int ret = one_iter1(samples200, best_off, 25, hz_for_cb, hz_for_cb); + return ret; + } + + return 0; +} + +// +// estimate SNR, yielding numbers vaguely similar to WSJT-X. +// m79 is a 79x8 complex FFT output. +// +float FT4::guess_snr(const FFTEngine::ffts_t &m79) +{ + int costas[] = {3, 1, 4, 0, 6, 5, 2}; + float pnoises = 0; + float psignals = 0; + + for (int i = 0; i < 7; i++) + { + psignals += std::abs(m79[i][costas[i]]); + psignals += std::abs(m79[36 + i][costas[i]]); + psignals += std::abs(m79[72 + i][costas[i]]); + pnoises += std::abs(m79[i][(costas[i] + 4) % 8]); + pnoises += std::abs(m79[36 + i][(costas[i] + 4) % 8]); + pnoises += std::abs(m79[72 + i][(costas[i] + 4) % 8]); + } + + for (int i = 0; i < 79; i++) + { + if (i < 7 || (i >= 36 && i < 36 + 7) || (i >= 72 && i < 72 + 7)) { + continue; + } + + std::vector v(8); + + for (int j = 0; j < 8; j++) { + v[j] = std::abs(m79[i][j]); + } + + std::sort(v.begin(), v.end()); + psignals += v[7]; // strongest tone, probably the signal + pnoises += (v[2] + v[3] + v[4]) / 3; + } + + pnoises /= 79; + psignals /= 79; + pnoises *= pnoises; // square yields power + psignals *= psignals; + float raw = psignals / pnoises; + raw -= 1; // turn (s+n)/n into s/n + + if (raw < 0.1) { + raw = 0.1; + } + + raw /= (2500.0 / 2.7); // 2.7 hz noise b/w -> 2500 hz b/w + float snr = 10 * log10(raw); + snr += 5; + snr *= 1.4; + return snr; +} + +// +// 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 FT4::fine(const FFTEngine::ffts_t &m79, int, float &adj_hz, float &adj_off) +{ + adj_hz = 0.0; + adj_off = 0.0; + // tone number for each of the 79 symbols. + int sym[79]; + float symval[79]; + float symphase[79]; + int costas[] = {3, 1, 4, 0, 6, 5, 2}; + + for (int i = 0; i < 79; i++) + { + if (i < 7) + { + sym[i] = costas[i]; + } + else if (i >= 36 && i < 36 + 7) + { + sym[i] = costas[i - 36]; + } + else if (i >= 72) + { + sym[i] = costas[i - 72]; + } + else + { + int mxj = -1; + float mx = 0; + + for (int j = 0; j < 8; j++) + { + float x = std::abs(m79[i][j]); + + if (mxj < 0 || x > mx) + { + mx = x; + mxj = j; + } + } + + sym[i] = mxj; + } + + symphase[i] = std::arg(m79[i][sym[i]]); + symval[i] = std::abs(m79[i][sym[i]]); + } + + float sum = 0; + float weight_sum = 0; + + for (int i = 0; i < 79 - 1; i++) + { + float d = symphase[i + 1] - symphase[i]; + + while (d > M_PI) { + d -= 2 * M_PI; + } + + while (d < -M_PI) { + d += 2 * M_PI; + } + + float w = symval[i]; + sum += d * w; + weight_sum += w; + } + + float mean = sum / weight_sum; + float err_rad = mean; // radians per symbol time + float err_hz = (err_rad / (2 * M_PI)) / 0.16; // cycles per symbol time + + // if each symbol's phase is a bit more than we expect, + // that means the real frequency is a bit higher + // than we thought, so increase our estimate. + adj_hz = err_hz; + + // + // now think about offset error. + // + // the higher tones have many cycles per + // symbol -- e.g. tone 7 has 11 cycles + // in each symbol. a one- or two-sample + // offset error at such a high tone will + // change the phase by pi or more, + // which makes the phase-to-samples + // conversion ambiguous. so only try + // to distinguish early-ontime-late, + // not the amount. + // + int nearly = 0; + int nlate = 0; + float early = 0.0; + float late = 0.0; + + for (int i = 1; i < 79; i++) + { + float ph0 = std::arg(m79[i - 1][sym[i - 1]]); + float ph = std::arg(m79[i][sym[i]]); + float d = ph - ph0; + d -= err_rad; // correct for hz error. + + while (d > M_PI) { + d -= 2 * M_PI; + } + + while (d < -M_PI) { + d += 2 * M_PI; + } + + // if off is correct, each symbol will have the same phase (modulo + // the above hz correction), since each FFT bin holds an integer + // number of cycles. + + // if off is too small, the phase is altered by the trailing part + // of the previous symbol. if the previous tone was lower, + // the phase won't have advanced as much as expected, and + // this symbol's phase will be lower than the previous phase. + // if the previous tone was higher, the phase will be more + // advanced than expected. thus off too small leads to + // a phase difference that's the reverse of the tone difference. + + // if off is too high, then the FFT started a little way into + // this symbol, which causes the phase to be advanced a bit. + // of course the previous symbol's phase was also advanced + // too much. if this tone is higher than the previous symbol, + // its phase will be more advanced than the previous. if + // less, less. + + // the point: if successive phases and tone differences + // are positively correlated, off is too high. if negatively, + // too low. + + // fine_max_tone: + // if late, ignore if a high tone, since ambiguous. + // if early, ignore if prev is a high tone. + + if (sym[i] > sym[i - 1]) + { + if (d > 0 && sym[i] <= params.fine_max_tone) + { + nlate++; + late += d / std::abs(sym[i] - sym[i - 1]); + } + + if (d < 0 && sym[i - 1] <= params.fine_max_tone) + { + nearly++; + early += fabs(d) / std::abs(sym[i] - sym[i - 1]); + } + } + else if (sym[i] < sym[i - 1]) + { + if (d > 0 && sym[i - 1] <= params.fine_max_tone) + { + nearly++; + early += d / std::abs(sym[i] - sym[i - 1]); + } + + if (d < 0 && sym[i] <= params.fine_max_tone) + { + nlate++; + late += fabs(d) / std::abs(sym[i] - sym[i - 1]); + } + } + } + + if (nearly > 0) { + early /= nearly; + } + + if (nlate > 0) { + late /= nlate; + } + + // qDebug("early %d %.1f, late %d %.1f", nearly, early, nlate, late); + + // assumes 32 samples/symbol. + if (nearly > 2 * nlate) + { + adj_off = round(32 * early / params.fine_thresh); + + if (adj_off > params.fine_max_off) { + adj_off = params.fine_max_off; + } + } + else if (nlate > 2 * nearly) + { + adj_off = 0 - round(32 * late / params.fine_thresh); + + if (fabs(adj_off) > params.fine_max_off) { + adj_off -= params.fine_max_off; + } + } +} + +// +// 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 FT4::one_iter1( + const std::vector &samples200x, + int best_off, + float best_hz, + float hz0_for_cb, + float hz1_for_cb +) +{ + // put best_hz in the middle of bin 4, at 25.0. + std::vector samples200 = shift200( + samples200x, + 0, + samples200x.size(), + best_hz + ); + + // mini 79x8 FFT. + FFTEngine::ffts_t m79 = extract(samples200, 25, best_off); + + // look at symbol-to-symbol phase change to try + // to improve best_hz and best_off. + if (params.do_fine_hz || params.do_fine_off) + { + float adj_hz = 0; + float adj_off = 0; + fine(m79, 4, adj_hz, adj_off); + + if (params.do_fine_hz == 0) { + adj_hz = 0; + } + + if (params.do_fine_off == 0) { + adj_off = 0; + } + + if (fabs(adj_hz) < 6.25 / 4 && fabs(adj_off) < 4) + { + best_hz += adj_hz; + best_off += round(adj_off); + + if (best_off < 0) { + best_off = 0; + } + + samples200 = shift200(samples200x, 0, samples200x.size(), best_hz); + m79 = extract(samples200, 25, best_off); + } + } + + float ll174[174]; + + if (params.soft_ones) + { + if (params.soft_ones == 1) { + soft_decode(m79, ll174); + } else { + c_soft_decode(m79, ll174); + } + + int ret = try_decode( + samples200, + ll174, + best_hz, + best_off, + hz0_for_cb, + hz1_for_cb, + params.use_osd, + "", + m79 + ); + + if (ret) { + return ret; + } + } + + if (params.soft_pairs) + { + float p174[174]; + soft_decode_pairs(m79, p174); + int ret = try_decode( + samples200, + p174, + best_hz, + best_off, + hz0_for_cb, + hz1_for_cb, + params.use_osd, + "", + m79 + ); + + if (ret) { + return ret; + } + + if (params.soft_ones == 0) { + std::copy(p174, p174 + 174, ll174); + } + } + + if (params.soft_triples) + { + float p174[174]; + soft_decode_triples(m79, p174); + int ret = try_decode( + samples200, + p174, + best_hz, + best_off, + hz0_for_cb, + hz1_for_cb, + params.use_osd, + "", + m79 + ); + + if (ret) { + return ret; + } + } + + if (params.use_hints) + { + for (int hi = 0; hi < (int)hints1_.size(); hi++) + { + int h = hints1_[hi]; // 28-bit number, goes in ll174 0..28 + + if (params.use_hints == 2 && h != 2) + { + // just CQ + continue; + } + + float n174[174]; + + for (int i = 0; i < 174; i++) + { + if (i < 28) + { + int bit = h & (1 << 27); + + if (bit) { + n174[i] = -4.97; + } else { + n174[i] = 4.97; + } + + h <<= 1; + } + else + { + n174[i] = ll174[i]; + } + } + + int ret = try_decode( + samples200, + n174, + best_hz, + best_off, + hz0_for_cb, + hz1_for_cb, + 0, + "hint1", + m79 + ); + + if (ret) { + return ret; + } + } + } + + if (params.use_hints == 1) + { + for (int hi = 0; hi < (int)hints2_.size(); hi++) + { + int h = hints2_[hi]; // 28-bit number, goes in ll174 29:29+28 + float n174[174]; + + for (int i = 0; i < 174; i++) + { + if (i >= 29 && i < 29 + 28) + { + int bit = h & (1 << 27); + + if (bit) { + n174[i] = -4.97; + } else { + n174[i] = 4.97; + } + + h <<= 1; + } + else + { + n174[i] = ll174[i]; + } + } + + int ret = try_decode( + samples200, + n174, + best_hz, + best_off, + hz0_for_cb, + hz1_for_cb, + 0, + "hint2", + m79 + ); + + if (ret) { + return ret; + } + } + } + + return 0; +} + +// +// 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 FT4::subtract( + const std::vector re79, + float hz0, + float hz1, + float off_sec +) +{ + int block = blocksize(rate_); + float bin_hz = rate_ / (float)block; + int off0 = off_sec * rate_; + float mhz = (hz0 + hz1) / 2.0; + int bin0 = round(mhz / bin_hz); + // move nsamples so that signal is centered in bin0. + float diff0 = (bin0 * bin_hz) - hz0; + float diff1 = (bin0 * bin_hz) - hz1; + std::vector moved = fftEngine_->hilbert_shift(nsamples_, diff0, diff1, rate_); + FFTEngine::ffts_t bins = fftEngine_->ffts(moved, off0, block); + + if (bin0 + 8 > (int)bins[0].size()) { + return; + } + + if ((int)bins.size() < 79) { + return; + + } + + std::vector phases(79); + std::vector amps(79); + + for (int i = 0; i < 79; i++) + { + int sym = bin0 + re79[i]; + std::complex c = bins[i][sym]; + phases[i] = std::arg(c); + // FFT multiplies magnitudes by number of bins, + // or half the number of samples. + amps[i] = std::abs(c) / (block / 2.0); + } + + int ramp = round(block * params.subtract_ramp); + + if (ramp < 1) { + ramp = 1; + } + + // initial ramp part of first symbol. + { + int sym = bin0 + re79[0]; + float phase = phases[0]; + float amp = amps[0]; + float hz = 6.25 * sym; + float dtheta = 2 * M_PI / (rate_ / hz); // advance per sample + + for (int jj = 0; jj < ramp; jj++) + { + float theta = phase + jj * dtheta; + float x = amp * cos(theta); + x *= jj / (float)ramp; + int iii = off0 + block * 0 + jj; + moved[iii] -= x; + } + } + + for (int si = 0; si < 79; si++) + { + int sym = bin0 + re79[si]; + float phase = phases[si]; + float amp = amps[si]; + float hz = 6.25 * sym; + float dtheta = 2 * M_PI / (rate_ / hz); // advance per sample + + // we've already done the first ramp for this symbol. + // now for the steady part between ramps. + for (int jj = ramp; jj < block - ramp; jj++) + { + float theta = phase + jj * dtheta; + float x = amp * cos(theta); + int iii = off0 + block * si + jj; + moved[iii] -= x; + } + + // now the two ramps, from us to the next symbol. + // we need to smoothly change the frequency, + // approximating wsjt-x's gaussian frequency shift, + // and also end up matching the next symbol's phase, + // which is often different from this symbol due + // to inaccuracies in hz or offset. + + // at start of this symbol's off-ramp. + float theta = phase + (block - ramp) * dtheta; + float hz1; + float phase1; + + if (si + 1 >= 79) + { + hz1 = hz; + phase1 = phase; + } + else + { + int sym1 = bin0 + re79[si + 1]; + hz1 = 6.25 * sym1; + phase1 = phases[si + 1]; + } + + float dtheta1 = 2 * M_PI / (rate_ / hz1); + // add this to dtheta for each sample, to gradually + // change the frequency. + float inc = (dtheta1 - dtheta) / (2.0 * ramp); + // after we've applied all those inc's, what will the + // phase be at the end of the next symbol's initial ramp, + // if we don't do anything to correct it? + float actual = theta + dtheta * 2.0 * ramp + inc * 4.0 * ramp * ramp / 2.0; + // what phase does the next symbol want to be at when + // its on-ramp finishes? + float target = phase1 + dtheta1 * ramp; + + // ??? + while (fabs(target - actual) > M_PI) + { + if (target < actual) { + target += (2 * M_PI) - 1e-3; // plus epsilonn to break possible infinite loop + } else { + target -= (2 * M_PI) + 1e-3; // plus epsilonn to break possible infinite loop + } + } + + // adj is to be spread evenly over the off-ramp and on-ramp samples. + float adj = target - actual; + int end = block + ramp; + + if (si == 79 - 1) { + end = block; + } + + for (int jj = block - ramp; jj < end; jj++) + { + int iii = off0 + block * si + jj; + float x = amp * cos(theta); + + // trail off to zero at the very end. + if (si == 79 - 1) { + x *= 1.0 - ((jj - (block - ramp)) / (float)ramp); + } + + moved[iii] -= x; + theta += dtheta; + dtheta += inc; + theta += adj / (2.0 * ramp); + } + } + + nsamples_ = fftEngine_->hilbert_shift(moved, -diff0, -diff1, rate_); +} + +// +// 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 FT4::try_decode( + const std::vector &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 +) +{ + int a174[174]; + std::string comment(comment1); + + if (decode(ll174, a174, params, use_osd, comment)) + { + // a174 is corrected 91 bits of plain message plus 83 bits of LDPC parity. + // how many of the corrected 174 bits match the received signal in ll174? + int correct_bits = 0; + + for (int i = 0; i < 174; i++) + { + if (ll174[i] < 0 && a174[i] == 1) { + correct_bits += 1; + } else if (ll174[i] > 0 && a174[i] == 0) { + correct_bits += 1; + } + } + + // reconstruct correct 79 symbols from LDPC output. + std::vector re79 = recode(a174); + + if (params.do_third == 1) + { + // fine-tune offset and hz for better subtraction. + float best_off = best_off_samples / 200.0; + search_both_known( + samples200, + 200, + re79, + best_hz, + best_off, + best_hz, + best_off + ); + best_off_samples = round(best_off * 200.0); + } + + // convert starting sample # from 200 samples/second back to rate_. + // also hz. + float best_off = best_off_samples / 200.0; // convert to seconds + best_hz = hz0_for_cb + (best_hz - 25.0); + + if (params.do_third == 2) + { + // fine-tune offset and hz for better subtraction. + search_both_known( + samples_, + rate_, + re79, + best_hz, + best_off, + best_hz, + best_off + ); + } + + float snr = guess_snr(m79); + + if (cb_) + { + cb_mu_.lock(); + int ret = cb_->hcb( + a174, + best_hz + down_hz_, + best_off, + comment.c_str(), + snr, + pass_, + correct_bits + ); + + cb_mu_.unlock(); + + if (ret == 2) + { + // a new decode. subtract it from nsamples_. + subtract(re79, best_hz, best_hz, best_off); + } + + return ret; + } + + return 1; + } + else + { + return 0; + } +} + +// +// 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 FT4::recode(int a174[]) +{ + int i174 = 0; + int costas[] = {3, 1, 4, 0, 6, 5, 2}; + std::vector out79; + + for (int i79 = 0; i79 < 79; i79++) + { + if (i79 < 7) + { + out79.push_back(costas[i79]); + } + else if (i79 >= 36 && i79 < 36 + 7) + { + out79.push_back(costas[i79 - 36]); + } + else if (i79 >= 72) + { + out79.push_back(costas[i79 - 72]); + } + else + { + int sym = (a174[i174 + 0] << 2) | (a174[i174 + 1] << 1) | (a174[i174 + 2] << 0); + i174 += 3; + // gray code + int map[] = {0, 1, 3, 2, 5, 6, 4, 7}; + sym = map[sym]; + out79.push_back(sym); + } + } + // assert(out79.size() == 79); + // assert(i174 == 174); + return out79; +} + +FT4DecoderLight::~FT4DecoderLight() { forceQuit(); } -void FT4Decoder::entry( +void FT4DecoderLight::entry( float xsamples[], int nsamples, int start, @@ -82,7 +3657,7 @@ void FT4Decoder::entry( } } -void FT4Decoder::wait(double time_left) +void FT4DecoderLight::wait(double time_left) { unsigned long thread_timeout = time_left * 1000; @@ -92,7 +3667,7 @@ void FT4Decoder::wait(double time_left) if (!success) { - qDebug("FT8::FT4Decoder::wait: thread timed out"); + qDebug("FT8::FT4DecoderLight::wait: thread timed out"); thread_timeout = 50; } @@ -100,7 +3675,7 @@ void FT4Decoder::wait(double time_left) } } -void FT4Decoder::forceQuit() +void FT4DecoderLight::forceQuit() { while (!threads.empty()) { @@ -110,4 +3685,139 @@ void FT4Decoder::forceQuit() } } +FT4Decoder::~FT4Decoder() +{ + forceQuit(); // stop all remaining running threads if any + + for (auto& fftEngine : fftEngines) { + delete fftEngine; + } +} + +// +// Launch decoding +// +void FT4Decoder::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 nprevdecs, + struct cdecode *xprevdecs +) +{ + double t0 = now(); + double deadline = t0 + time_left; + double final_deadline = t0 + total_time_left; + // decodes from previous runs, for subtraction. + std::vector prevdecs; + + for (int i = 0; i < nprevdecs; i++) { + prevdecs.push_back(xprevdecs[i]); + } + + std::vector samples(nsamples); + + for (int i = 0; i < nsamples; i++) { + samples[i] = xsamples[i]; + } + + if (min_hz < 0) { + min_hz = 0; + } + + if (max_hz > rate / 2) { + max_hz = rate / 2; + } + + float per = (max_hz - min_hz) / params.nthreads; + + for (int i = 0; i < params.nthreads; i++) + { + float hz0 = min_hz + i * per; + + if (i > 0 || params.overlap_edges) { + hz0 -= params.overlap; + } + + float hz1 = min_hz + (i + 1) * per; + + if (i != params.nthreads - 1 || params.overlap_edges) { + hz1 += params.overlap; + } + + hz0 = std::max(hz0, 0.0f); + hz1 = std::min(hz1, (rate / 2.0f) - 50); + + if (i == (int) fftEngines.size()) { + fftEngines.push_back(new FFTEngine()); + } + + FT4 *ft4 = new FT4( + samples, + hz0, + hz1, + start, + rate, + hints1, + hints2, + deadline, + final_deadline, + cb, + prevdecs, + fftEngines[i] + ); + + ft4->getParams() = getParams(); // transfer parameters + int npasses = nprevdecs > 0 ? params.npasses_two : params.npasses_one; + ft4->set_npasses(npasses); + QThread *th = new QThread(); + th->setObjectName(tr("ft4:%1:%2").arg(cb->get_name()).arg(i)); + threads.push_back(th); + // std::thread *th = new std::thread([ft4, npasses] () { ft4->go(npasses); }); + // thv.push_back(std::pair(ft4, th)); + ft4->moveToThread(th); + QObject::connect(th, &QThread::started, ft4, &FT4::start_work); + QObject::connect(ft4, &FT4::finished, th, &QThread::quit, Qt::DirectConnection); + QObject::connect(th, &QThread::finished, ft4, &QObject::deleteLater); + QObject::connect(th, &QThread::finished, th, &QThread::deleteLater); + th->start(); + } +} + +void FT4Decoder::wait(double time_left) +{ + unsigned long thread_timeout = time_left * 1000; + + while (threads.size() != 0) + { + bool success = threads.front()->wait(thread_timeout); + + if (!success) + { + qDebug("FT4::FT4Decoder::wait: thread timed out"); + thread_timeout = 50; // only 50ms for the rest + } + + threads.erase(threads.begin()); + } +} + +void FT4Decoder::forceQuit() +{ + while (threads.size() != 0) + { + threads.front()->quit(); + threads.front()->wait(); + threads.erase(threads.begin()); + } +} + } // namespace FT8 diff --git a/ft8/ft4.h b/ft8/ft4.h index 57f6724cb..fb934a4ed 100644 --- a/ft8/ft4.h +++ b/ft8/ft4.h @@ -17,8 +17,163 @@ class QThread; namespace FT8 { +// 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 class FT8_API FT4Params { +public: + 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 use_osd; + 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; + + FT4Params() + { + 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; // Gaussian + problt_how_sig = 0; // Gaussian + use_apriori = 1; + use_hints = 2; // 1 means use all hints, 2 means just CQ hints + win_type = 1; + use_osd = 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 FT4Params + +class FT8_API FT4ParamsLight +{ public: int nthreads; int ldpc_iters; @@ -27,7 +182,7 @@ public: int osd_ldpc_thresh; int max_candidates; - FT4Params() : + FT4ParamsLight() : nthreads(8), ldpc_iters(25), use_osd(1), @@ -37,9 +192,473 @@ public: {} }; -class FT8_API FT4Decoder : public QObject +// The FT8 worker +class FT8_API FT4 : public QObject { Q_OBJECT +public: + FT4( + const std::vector &samples, + float min_hz, + float max_hz, + int start, + int rate, + int hints1[], + int hints2[], + double deadline, + double final_deadline, + CallbackInterface *cb, + std::vector prevdecs, + FFTEngine *fftEngine + ); + ~FT4(); + // Number of passes + void set_npasses(int npasses) { npasses_ = npasses; } + // Start the worker + void start_work(); + // 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 coarse(const FFTEngine::ffts_t &bins, int si0, int si1); + + FT4Params& getParams() { return params; } + // + // given log likelihood for each bit, try LDPC and OSD decoders. + // on success, puts corrected 174 bits into a174[]. + // + static int decode(const float ll174[], int a174[], FT4Params& params, int use_osd, std::string &comment); + // encode a 77 bit message into a 174 bit payload + // adds the 14 bit CRC to obtain 91 bits + // apply (174, 91) generator mastrix to obtain the 83 parity bits + // append the 83 bits to the 91 bits message e+ crc to obtain the 174 bit payload + static void encode(int a174[], const int s77[]); + + // + // set ones and zero symbol indexes + // + static void set_ones_zeroes(int ones[], int zeroes[], int nbBits, int bitIndex); + + // + // mags is the vector of 2^nbSymbolBits vector of magnitudes at each symbol time + // ll174 is the resulting 174 soft bits of payload + // used in FT-chirp modulation scheme - generalized to any number of symbol bits + // + static void soft_decode_mags(FT4Params& params, const std::vector>& mags, int nbSymbolBits, float ll174[]); + + // + // Generic Gray decoding for magnitudes (floats) + // + static std::vector> un_gray_code_r_gen(const std::vector> &mags); + +private: + // + // 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 reduce_rate( + const std::vector &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 &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 &samples, + int rate, + const std::vector &syms, + float hz, + int off + ); + int search_time_fine( + const std::vector &samples200, + int off0, + int offN, + float hz, + int gran, + float &str + ); + int search_time_fine_known( + const std::vector> &bins, + int rate, + const std::vector &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 search_both( + const std::vector &samples200, + float hz0, + int hz_n, + float hz_win, + int off0, + int off_n, + int off_win + ); + void search_both_known( + const std::vector &samples, + int rate, + const std::vector &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 fft_shift( + const std::vector &samples, + int off, + int len, + int rate, + float hz + ); + // + // shift down by hz. + // + std::vector fft_shift_f( + const std::vector> &bins, + int rate, + float hz + ); + // shift the frequency by a fraction of 6.25, + // to center hz on bin 4 (25 hz). + std::vector shift200( + const std::vector &samples200, + int off, + int len, + float hz + ); + // returns a mini-FFT of 79 8-tone symbols. + FFTEngine::ffts_t extract(const std::vector &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> un_gray_code_r(const std::vector> &m79); + // + // normalize levels by windowed median. + // this helps, but why? + // + std::vector> convert_to_snr(const std::vector> &m79); + // + // normalize levels by windowed median. + // this helps, but why? + // + static std::vector> convert_to_snr_gen(const FT4Params& params, int nbSymbolBits, const std::vector> &mags); + // + // normalize levels by windowed median. + // this helps, but why? + // + std::vector>> c_convert_to_snr( + const std::vector>> &m79 + ); + // + // statistics to decide soft probabilities, + // to drive LDPC decoder. + // distribution of strongest tones, and + // distribution of noise. + // + static void make_stats( + const std::vector> &m79, + Stats &bests, + Stats &all + ); + // + // generalized version of the above for any number of symbols and no Costas + // used by FT-chirp decoder + // + static void make_stats_gen( + const std::vector> &mags, + int nbSymbolBits, + 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> 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. + // + static float bayes( + FT4Params& params, + 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 extract_bits(const std::vector &syms, const std::vector str); + // decode successive pairs of symbols. exploits the likelihood + // 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[] + ); + // + // 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> fbandpass( + const std::vector> &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 down_v7(const std::vector &samples, float hz); + std::vector down_v7_f(const std::vector> &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_merge(const std::vector> &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 &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 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 &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 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 &samples200x, + int best_off, + float best_hz, + float hz0_for_cb, + float hz1_for_cb + ); + +signals: + void finished(); +private: + FT4Params params; + FFTEngine *fftEngine_; + int npasses_; + static const double apriori174[]; + + float min_hz_; + float max_hz_; + std::vector samples_; // input to each pass + std::vector 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 hints1_; + std::vector hints2_; + int pass_; + float down_hz_; + + QMutex cb_mu_; + CallbackInterface *cb_; // call-back interface + + QMutex hack_mu_; + int hack_size_; + int hack_off_; + int hack_len_; + float hack_0_; + float hack_1_; + const float *hack_data_; + std::vector> hack_bins_; + std::vector prevdecs_; +}; // class FT4 + + +class FT8_API FT4DecoderLight : public QObject +{ + Q_OBJECT +public: + ~FT4DecoderLight(); + 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 * + ); + void wait(double time_left); + void forceQuit(); + FT4ParamsLight& getParams() { return params; } + +private: + FT4ParamsLight params; + std::vector threads; +}; + +class FT8_API FT4Decoder : public QObject { + Q_OBJECT public: ~FT4Decoder(); void entry( @@ -57,14 +676,15 @@ public: int, struct cdecode * ); - void wait(double time_left); - void forceQuit(); + void wait(double time_left); //!< wait for all threads to finish + void forceQuit(); //!< force quit all threads FT4Params& getParams() { return params; } - private: FT4Params params; std::vector threads; -}; + std::vector fftEngines; +}; // FT4Decoder + } // namespace FT8 diff --git a/ft8/ft4worker.cpp b/ft8/ft4worker.cpp index 8ebf414dd..a27f556a5 100644 --- a/ft8/ft4worker.cpp +++ b/ft8/ft4worker.cpp @@ -62,7 +62,7 @@ FT4Worker::FT4Worker( int start, int rate, CallbackInterface *cb, - const FT4Params& params + const FT4ParamsLight& params ) : m_samples(samples), m_minHz(minHz), @@ -369,6 +369,7 @@ void FT4Worker::decode() int plain[174]; int ldpcOk = 0; + int correct_bits = 0; LDPC::ldpc_decode(llr.data(), m_params.ldpc_iters, plain, &ldpcOk); if (ldpcOk < ldpcThreshold) @@ -378,8 +379,15 @@ void FT4Worker::decode() int a174[174]; - for (int i = 0; i < 174; i++) { + for (int i = 0; i < 174; i++) + { a174[i] = plain[i]; + + if (llr.data()[i] < 0 && a174[i] == 1) { + correct_bits += 1; + } else if (llr.data()[i] > 0 && a174[i] == 0) { + correct_bits += 1; + } } bool crcOk = OSD::check_crc(a174); @@ -420,7 +428,7 @@ void FT4Worker::decode() } const float off = static_cast(m_start) / m_rate + static_cast(candidate.start) / m_rate; - m_cb->hcb(a91, tone0, off, "FT4-EXP", snr, 0, ldpcOk); + m_cb->hcb(a91, tone0, off, "FT4-EXP", snr, 0, correct_bits); } } diff --git a/ft8/ft4worker.h b/ft8/ft4worker.h index dc74b9df5..82fbb36a4 100644 --- a/ft8/ft4worker.h +++ b/ft8/ft4worker.h @@ -40,7 +40,7 @@ public: int start, int rate, CallbackInterface *cb, - const FT4Params& params + const FT4ParamsLight& params ); void start_work(); @@ -59,7 +59,7 @@ private: int m_start; int m_rate; CallbackInterface *m_cb; - FT4Params m_params; + FT4ParamsLight m_params; FFTEngine m_fftEngine; }; diff --git a/ft8/ft8.cpp b/ft8/ft8.cpp index 85894e3a8..fa2bf6636 100644 --- a/ft8/ft8.cpp +++ b/ft8/ft8.cpp @@ -47,78 +47,6 @@ namespace FT8 { -// -// return a Hamming window of length n. -// -std::vector hamming(int n) -{ - std::vector h(n); - - for (int k = 0; k < n; k++) { - h[k] = 0.54 - 0.46 * cos(2 * M_PI * k / (n - 1.0)); - } - - return h; -} - -// -// blackman window -// -std::vector blackman(int n) -{ - std::vector h(n); - - for (int k = 0; k < n; k++) { - h[k] = 0.42 - 0.5 * cos(2 * M_PI * k / n) + 0.08 * cos(4 * M_PI * k / n); - } - - return h; -} - -// -// symmetric blackman window -// -std::vector sym_blackman(int n) -{ - std::vector h(n); - - for (int k = 0; k < (n / 2) + 1; k++) { - h[k] = 0.42 - 0.5 * cos(2 * M_PI * k / n) + 0.08 * cos(4 * M_PI * k / n); - } - - for (int k = n - 1; k >= (n / 2) + 1; --k) { - h[k] = h[(n - 1) - k]; - } - - return h; -} - -// -// blackman-harris window -// -std::vector blackmanharris(int n) -{ - float a0 = 0.35875; - float a1 = 0.48829; - float a2 = 0.14128; - float a3 = 0.01168; - std::vector h(n); - - for (int k = 0; k < n; k++) - { - // symmetric - h[k] = a0 - a1 * cos(2 * M_PI * k / (n - 1)) + a2 * cos(4 * M_PI * k / (n - 1)) - a3 * cos(6 * M_PI * k / (n - 1)); - // periodic - // h[k] = - // a0 - // - a1 * cos(2 * M_PI * k / n) - // + a2 * cos(4 * M_PI * k / n) - // - a3 * cos(6 * M_PI * k / n); - } - - return h; -} - // a-priori probability of each of the 174 LDPC codeword // bits being one. measured from reconstructed correct // codewords, into ft8bits, then python bprob.py. diff --git a/ft8/util.cpp b/ft8/util.cpp index 8fa0f6f0b..1f823a023 100644 --- a/ft8/util.cpp +++ b/ft8/util.cpp @@ -338,4 +338,76 @@ std::string trim(const std::string &s) { return rtrim(ltrim(s)); } +// +// return a Hamming window of length n. +// +std::vector hamming(int n) +{ + std::vector h(n); + + for (int k = 0; k < n; k++) { + h[k] = 0.54 - 0.46 * cos(2 * M_PI * k / (n - 1.0)); + } + + return h; +} + +// +// blackman window +// +std::vector blackman(int n) +{ + std::vector h(n); + + for (int k = 0; k < n; k++) { + h[k] = 0.42 - 0.5 * cos(2 * M_PI * k / n) + 0.08 * cos(4 * M_PI * k / n); + } + + return h; +} + +// +// symmetric blackman window +// +std::vector sym_blackman(int n) +{ + std::vector h(n); + + for (int k = 0; k < (n / 2) + 1; k++) { + h[k] = 0.42 - 0.5 * cos(2 * M_PI * k / n) + 0.08 * cos(4 * M_PI * k / n); + } + + for (int k = n - 1; k >= (n / 2) + 1; --k) { + h[k] = h[(n - 1) - k]; + } + + return h; +} + +// +// blackman-harris window +// +std::vector blackmanharris(int n) +{ + float a0 = 0.35875; + float a1 = 0.48829; + float a2 = 0.14128; + float a3 = 0.01168; + std::vector h(n); + + for (int k = 0; k < n; k++) + { + // symmetric + h[k] = a0 - a1 * cos(2 * M_PI * k / (n - 1)) + a2 * cos(4 * M_PI * k / (n - 1)) - a3 * cos(6 * M_PI * k / (n - 1)); + // periodic + // h[k] = + // a0 + // - a1 * cos(2 * M_PI * k / n) + // + a2 * cos(4 * M_PI * k / n) + // - a3 * cos(6 * M_PI * k / n); + } + + return h; +} + } // namespace FT8 diff --git a/ft8/util.h b/ft8/util.h index aad571560..371e348b7 100644 --- a/ft8/util.h +++ b/ft8/util.h @@ -50,6 +50,11 @@ std::vector gfsk_r( std::vector gfsk_window(int samples_per_symbol, float b); std::string trim(const std::string &s); +std::vector hamming(int n); // return a Hamming window of length n. +std::vector blackman(int n); // blackman window +std::vector sym_blackman(int n); // symmetric blackman window +std::vector blackmanharris(int n); // blackman-harris window + typedef unsigned long ulong; typedef unsigned int uint; diff --git a/plugins/channelrx/demodft8/ft8demodworker.h b/plugins/channelrx/demodft8/ft8demodworker.h index af749827d..d11f72eb9 100644 --- a/plugins/channelrx/demodft8/ft8demodworker.h +++ b/plugins/channelrx/demodft8/ft8demodworker.h @@ -112,7 +112,7 @@ private: bool m_invalidSequence; qint64 m_baseFrequency; FT8::FT8Decoder m_ft8Decoder; - FT8::FT4Decoder m_ft4Decoder; + FT8::FT4DecoderLight m_ft4Decoder; FT8::Packing m_packing; MessageQueue *m_guiReportingMessageQueue; MessageQueue *m_pskReportingMessageQueue; diff --git a/sdrbench/test_ft4.cpp b/sdrbench/test_ft4.cpp index 90ba74922..ee1e59ec7 100644 --- a/sdrbench/test_ft4.cpp +++ b/sdrbench/test_ft4.cpp @@ -216,7 +216,7 @@ void MainBench::testFT4(const QString& wavFile, const QString& argsStr) wfile.close(); - FT8::FT4Decoder decoder; + FT8::FT4DecoderLight decoder; decoder.getParams().nthreads = nthreads; decoder.getParams().use_osd = 1; decoder.getParams().osd_depth = 6;