mirror of
https://github.com/f4exb/sdrangel.git
synced 2024-11-21 23:55:13 -05:00
3827 lines
99 KiB
C++
3827 lines
99 KiB
C++
///////////////////////////////////////////////////////////////////////////////////
|
|
// Copyright (C) 2023 Edouard Griffiths, F4EXB <f4exb06@gmail.com> //
|
|
// //
|
|
// This is the code from ft8mon: https://github.com/rtmrtmrtmrtm/ft8mon //
|
|
// reformatted and adapted to Qt and SDRangel context //
|
|
// //
|
|
// This program is free software; you can redistribute it and/or modify //
|
|
// it under the terms of the GNU General Public License as published by //
|
|
// the Free Software Foundation as version 3 of the License, or //
|
|
// (at your option) any later version. //
|
|
// //
|
|
// This program is distributed in the hope that it will be useful, //
|
|
// but WITHOUT ANY WARRANTY; without even the implied warranty of //
|
|
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the //
|
|
// GNU General Public License V3 for more details. //
|
|
// //
|
|
// You should have received a copy of the GNU General Public License //
|
|
// along with this program. If not, see <http://www.gnu.org/licenses/>. //
|
|
///////////////////////////////////////////////////////////////////////////////////
|
|
|
|
//
|
|
// An FT8 receiver in C++.
|
|
//
|
|
// Many ideas and protocol details borrowed from Franke
|
|
// and Taylor's WSJT-X code.
|
|
//
|
|
// Robert Morris, AB1HL
|
|
//
|
|
|
|
#include <stdio.h>
|
|
// #include <assert.h>
|
|
#include <math.h>
|
|
#include <fftw3.h>
|
|
#include <algorithm>
|
|
#include <complex>
|
|
#include <random>
|
|
#include <functional>
|
|
#include <map>
|
|
|
|
#include <QThread>
|
|
|
|
#include "util.h"
|
|
#include "ft8.h"
|
|
#include "libldpc.h"
|
|
#include "osd.h"
|
|
#include "arrays.h"
|
|
|
|
namespace FT8 {
|
|
|
|
//
|
|
// return a Hamming window of length n.
|
|
//
|
|
std::vector<float> hamming(int n)
|
|
{
|
|
std::vector<float> 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<float> blackman(int n)
|
|
{
|
|
std::vector<float> 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<float> sym_blackman(int n)
|
|
{
|
|
std::vector<float> 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<float> blackmanharris(int n)
|
|
{
|
|
float a0 = 0.35875;
|
|
float a1 = 0.48829;
|
|
float a2 = 0.14128;
|
|
float a3 = 0.01168;
|
|
std::vector<float> 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.
|
|
// from ft8-n4
|
|
const double FT8::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,
|
|
};
|
|
|
|
FT8::FT8(
|
|
const std::vector<float> &samples,
|
|
float min_hz,
|
|
float max_hz,
|
|
int start,
|
|
int rate,
|
|
int hints1[],
|
|
int hints2[],
|
|
double deadline,
|
|
double final_deadline,
|
|
CallbackInterface *cb,
|
|
std::vector<cdecode> prevdecs,
|
|
FFTEngine *fftEngine
|
|
)
|
|
{
|
|
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;
|
|
}
|
|
|
|
FT8::~FT8()
|
|
{
|
|
}
|
|
|
|
void FT8::start_work()
|
|
{
|
|
go(npasses_);
|
|
emit finished();
|
|
}
|
|
|
|
// strength of costas block of signal with tone 0 at bi0,
|
|
// and symbol zero at si0.
|
|
float FT8::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 FT8::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<Strength> FT8::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<Strength> strengths;
|
|
|
|
for (int bi = min_bin; bi < max_bin && bi + 8 <= nbins; bi++)
|
|
{
|
|
std::vector<Strength> 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<float> FT8::reduce_rate(
|
|
const std::vector<float> &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<std::complex<float>> 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<std::complex<float>> 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<float> vvv = fftEngine_->one_ifft(bbins);
|
|
delta_hz = delta * bin_hz;
|
|
|
|
return vvv;
|
|
}
|
|
|
|
void FT8::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<int> distribution(0, samples_.size() - 1);
|
|
auto rnd = std::bind(distribution, generator);
|
|
std::vector<float> 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<int> 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<Strength> 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<std::complex<float>> 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<float> 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<Strength> 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 FT8::one_strength(const std::vector<float> &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 FT8::one_strength_known(
|
|
const std::vector<float> &samples,
|
|
int rate,
|
|
const std::vector<int> &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<float> 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<float> 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 FT8::search_time_fine(
|
|
const std::vector<float> &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<float> 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 FT8::search_time_fine_known(
|
|
const std::vector<std::complex<float>> &bins,
|
|
int rate,
|
|
const std::vector<int> &syms,
|
|
int off0,
|
|
int offN,
|
|
float hz,
|
|
int gran,
|
|
float &str
|
|
)
|
|
{
|
|
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<float> 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<Strength> FT8::search_both(
|
|
const std::vector<float> &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<Strength> 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 FT8::search_both_known(
|
|
const std::vector<float> &samples,
|
|
int rate,
|
|
const std::vector<int> &syms,
|
|
float hz0,
|
|
float off_secs0, // seconds
|
|
float &hz_out,
|
|
float &off_out
|
|
)
|
|
{
|
|
// 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<std::complex<float>> 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<float> FT8::fft_shift(
|
|
const std::vector<float> &samples,
|
|
int off,
|
|
int len,
|
|
int rate,
|
|
float hz
|
|
)
|
|
{
|
|
std::vector<std::complex<float>> 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<float> FT8::fft_shift_f(
|
|
const std::vector<std::complex<float>> &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<std::complex<float>> 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<float> 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<float> FT8::shift200(
|
|
const std::vector<float> &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 FT8::extract(const std::vector<float> &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 FT8::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<std::vector<float>> FT8::un_gray_code_r(const std::vector<std::vector<float>> &m79)
|
|
{
|
|
std::vector<std::vector<float>> 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<std::vector<float>> FT8::un_gray_code_r_gen(const std::vector<std::vector<float>> &mags)
|
|
{
|
|
if (mags.size() == 0) {
|
|
return mags;
|
|
}
|
|
|
|
std::vector<std::vector<float>> 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<std::vector<float>> FT8::convert_to_snr(const std::vector<std::vector<float>> &m79)
|
|
{
|
|
if (params.snr_how < 0 || params.snr_win < 0) {
|
|
return m79;
|
|
}
|
|
|
|
//
|
|
// for each symbol time, what's its "noise" level?
|
|
//
|
|
std::vector<float> mm(79);
|
|
|
|
for (int si = 0; si < 79; si++)
|
|
{
|
|
std::vector<float> 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<float> winwin;
|
|
|
|
if (params.snr_win > 0) {
|
|
winwin = blackman(2 * params.snr_win + 1);
|
|
} else {
|
|
winwin.push_back(1.0);
|
|
}
|
|
|
|
std::vector<std::vector<float>> 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<std::vector<std::complex<float>>> FT8::c_convert_to_snr(
|
|
const std::vector<std::vector<std::complex<float>>> &m79
|
|
)
|
|
{
|
|
if (params.snr_how < 0 || params.snr_win < 0) {
|
|
return m79;
|
|
}
|
|
|
|
//
|
|
// for each symbol time, what's its "noise" level?
|
|
//
|
|
std::vector<float> mm(79);
|
|
|
|
for (int si = 0; si < 79; si++)
|
|
{
|
|
std::vector<float> 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<float> winwin;
|
|
|
|
if (params.snr_win > 0) {
|
|
winwin = blackman(2 * params.snr_win + 1);
|
|
} else {
|
|
winwin.push_back(1.0);
|
|
}
|
|
|
|
std::vector<std::vector<std::complex<float>>> 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<std::vector<float>> FT8::convert_to_snr_gen(const FT8Params& params, int nbSymbolBits, const std::vector<std::vector<float>> &mags)
|
|
{
|
|
if (params.snr_how < 0 || params.snr_win < 0) {
|
|
return mags;
|
|
}
|
|
|
|
//
|
|
// for each symbol time, what's its "noise" level?
|
|
//
|
|
std::vector<float> mm(mags.size());
|
|
int nbSymbols = 1<<nbSymbolBits;
|
|
|
|
for (int si = 0; si < (int) mags.size(); si++)
|
|
{
|
|
std::vector<float> 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<float> winwin;
|
|
|
|
if (params.snr_win > 0) {
|
|
winwin = blackman(2 * params.snr_win + 1);
|
|
} else {
|
|
winwin.push_back(1.0);
|
|
}
|
|
|
|
std::vector<std::vector<float>> 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 FT8::make_stats(
|
|
const std::vector<std::vector<float>> &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 FT8::make_stats_gen(
|
|
const std::vector<std::vector<float>> &mags,
|
|
int nbSymbolBits,
|
|
Stats &bests,
|
|
Stats &all
|
|
)
|
|
{
|
|
int nbBins = 1<<nbSymbolBits;
|
|
|
|
for (unsigned int si = 0; si < mags.size(); si++)
|
|
{
|
|
float mx = 0;
|
|
|
|
for (int bi = 0; bi < nbBins; bi++)
|
|
{
|
|
float x = mags[si][bi];
|
|
|
|
if (x > 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<std::vector<float>> FT8::soft_c2m(const FFTEngine::ffts_t &c79)
|
|
{
|
|
std::vector<std::vector<float>> m79(79);
|
|
std::vector<float> 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<float> phases(79);
|
|
|
|
// for each symbol time, median of nearby phases
|
|
for (int si = 0; si < 79; si++)
|
|
{
|
|
std::vector<float> 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 FT8::bayes(
|
|
FT8Params& 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 FT8::soft_decode(const FFTEngine::ffts_t &c79, float ll174[])
|
|
{
|
|
std::vector<std::vector<float>> 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 FT8::soft_decode_mags(FT8Params& params, const std::vector<std::vector<float>>& mags_, int nbSymbolBits, float ll174[])
|
|
{
|
|
if ((nbSymbolBits > 16) || (nbSymbolBits < 1)) {
|
|
return;
|
|
}
|
|
|
|
std::vector<std::vector<float>> 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<int> zeroi(zoX*zoY);
|
|
std::vector<int> 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 FT8::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<float> maxes[79];
|
|
|
|
for (int i = 0; i < 79; i++)
|
|
{
|
|
std::complex<float> 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<std::vector<float>> m79(79);
|
|
|
|
for (int i = 0; i < 79; i++)
|
|
{
|
|
m79[i].resize(8);
|
|
|
|
for (int j = 0; j < 8; j++)
|
|
{
|
|
std::complex<float> 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<float> c1 = maxes[k];
|
|
std::complex<float> 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 FT8::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<<nbBits) - 1;
|
|
int maskLow = (1<<bitIndex) - 1;
|
|
int maskHigh = mask ^ maskLow;
|
|
|
|
for (int i = 0; i < nbIndexes; i++)
|
|
{
|
|
zeroes[i] = (i & maskLow) + ((i & maskHigh)<<1);
|
|
ones[i] = zeroes[i] + (1<<bitIndex);
|
|
}
|
|
}
|
|
}
|
|
|
|
//
|
|
// turn 79 symbol numbers into 174 bits.
|
|
// strip out the three Costas sync blocks,
|
|
// leaving 58 symbol numbers.
|
|
// each represents three bits.
|
|
// (all post-un-gray-code).
|
|
// str is per-symbol strength; must be positive.
|
|
// each returned element is < 0 for 1, > 0 for zero,
|
|
// scaled by str.
|
|
//
|
|
std::vector<float> FT8::extract_bits(const std::vector<int> &syms, const std::vector<float> str)
|
|
{
|
|
// assert(syms.size() == 79);
|
|
// assert(str.size() == 79);
|
|
|
|
std::vector<float> 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 FT8::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> 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<float> 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 FT8::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> 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<float> 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 FT8::decode(const float ll174[], int a174[], FT8Params& _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 FT8::encode(int a174[], 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<std::complex<float>> FT8::fbandpass(
|
|
const std::vector<std::complex<float>> &bins0,
|
|
float bin_hz,
|
|
float low_outer, // start of transition
|
|
float low_inner, // start of flat area
|
|
float high_inner, // end of flat area
|
|
float high_outer // end of transition
|
|
)
|
|
{
|
|
// assert(low_outer <= low_inner);
|
|
// assert(low_inner <= high_inner);
|
|
// assert(high_inner <= high_outer);
|
|
|
|
int nbins = bins0.size();
|
|
std::vector<std::complex<float>> 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<float> FT8::down_v7(const std::vector<float> &samples, float hz)
|
|
{
|
|
int len = samples.size();
|
|
std::vector<std::complex<float>> bins = fftEngine_->one_fft(samples, 0, len);
|
|
|
|
return down_v7_f(bins, len, hz);
|
|
}
|
|
|
|
std::vector<float> FT8::down_v7_f(const std::vector<std::complex<float>> &bins, int len, float hz)
|
|
{
|
|
int nbins = bins.size();
|
|
float bin_hz = rate_ / (float)len;
|
|
int down = round((hz - 25) / bin_hz);
|
|
std::vector<std::complex<float>> 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<std::complex<float>> bbins(blen / 2 + 1);
|
|
|
|
for (int i = 0; i < (int)bbins.size(); i++) {
|
|
bbins[i] = bins1[i];
|
|
}
|
|
|
|
std::vector<float> 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 FT8::one_merge(const std::vector<std::complex<float>> &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<float> 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 FT8::one_iter(const std::vector<float> &samples200, int best_off, float hz_for_cb)
|
|
{
|
|
if (params.do_second)
|
|
{
|
|
std::vector<Strength> 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 FT8::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<float> 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 FT8::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 FT8::one_iter1(
|
|
const std::vector<float> &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<float> 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 FT8::subtract(
|
|
const std::vector<int> 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<float> 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<float> phases(79);
|
|
std::vector<float> amps(79);
|
|
|
|
for (int i = 0; i < 79; i++)
|
|
{
|
|
int sym = bin0 + re79[i];
|
|
std::complex<float> 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 FT8::try_decode(
|
|
const std::vector<float> &samples200,
|
|
float ll174[174],
|
|
float best_hz,
|
|
int best_off_samples,
|
|
float hz0_for_cb,
|
|
float,
|
|
int use_osd,
|
|
const char *comment1,
|
|
const FFTEngine::ffts_t &m79
|
|
)
|
|
{
|
|
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<int> 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<int> FT8::recode(int a174[])
|
|
{
|
|
int i174 = 0;
|
|
int costas[] = {3, 1, 4, 0, 6, 5, 2};
|
|
std::vector<int> 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;
|
|
}
|
|
|
|
FT8Decoder::~FT8Decoder()
|
|
{
|
|
forceQuit(); // stop all remaining running threads if any
|
|
|
|
for (auto& fftEngine : fftEngines) {
|
|
delete fftEngine;
|
|
}
|
|
}
|
|
|
|
//
|
|
// Launch decoding
|
|
//
|
|
void FT8Decoder::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<cdecode> prevdecs;
|
|
|
|
for (int i = 0; i < nprevdecs; i++) {
|
|
prevdecs.push_back(xprevdecs[i]);
|
|
}
|
|
|
|
std::vector<float> 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());
|
|
}
|
|
|
|
FT8 *ft8 = new FT8(
|
|
samples,
|
|
hz0,
|
|
hz1,
|
|
start,
|
|
rate,
|
|
hints1,
|
|
hints2,
|
|
deadline,
|
|
final_deadline,
|
|
cb,
|
|
prevdecs,
|
|
fftEngines[i]
|
|
);
|
|
|
|
ft8->getParams() = getParams(); // transfer parameters
|
|
int npasses = nprevdecs > 0 ? params.npasses_two : params.npasses_one;
|
|
ft8->set_npasses(npasses);
|
|
QThread *th = new QThread();
|
|
th->setObjectName(tr("ft8:%1:%2").arg(cb->get_name()).arg(i));
|
|
threads.push_back(th);
|
|
// std::thread *th = new std::thread([ft8, npasses] () { ft8->go(npasses); });
|
|
// thv.push_back(std::pair<FT8*, std::thread*>(ft8, th));
|
|
ft8->moveToThread(th);
|
|
QObject::connect(th, &QThread::started, ft8, &FT8::start_work);
|
|
QObject::connect(ft8, &FT8::finished, th, &QThread::quit, Qt::DirectConnection);
|
|
QObject::connect(th, &QThread::finished, ft8, &QObject::deleteLater);
|
|
QObject::connect(th, &QThread::finished, th, &QThread::deleteLater);
|
|
th->start();
|
|
}
|
|
}
|
|
|
|
void FT8Decoder::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("FT8::FT8Decoder::wait: thread timed out");
|
|
thread_timeout = 50; // only 50ms for the rest
|
|
}
|
|
|
|
threads.erase(threads.begin());
|
|
}
|
|
}
|
|
|
|
void FT8Decoder::forceQuit()
|
|
{
|
|
while (threads.size() != 0)
|
|
{
|
|
threads.front()->quit();
|
|
threads.front()->wait();
|
|
threads.erase(threads.begin());
|
|
}
|
|
}
|
|
|
|
} // namespace FT8
|