VE7UWU/sdrangel
1
0
Fork 0
mirror of https://github.com/f4exb/sdrangel.git synced 2024-11-25 01:18:38 -05:00
Code Issues Projects Releases Wiki Activity

FT8: added soft decode on magnitudes and genaralize to any number of bits per symbol

Browse Source
This commit is contained in:
f4exb 2024-04-03 12:54:03 +02:00
parent 0b780b158b
commit 8e77ad01ad
11 changed files with 853 additions and 258 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
ft8/CMakeLists.txt
View File

@ -7,6 +7,7 @@ set(ft8_SOURCES
ft8.cpp
ft8plan.cpp
ft8plans.cpp
ft8stats.cpp
libldpc.cpp
osd.cpp
packing.cpp
@ -22,6 +23,7 @@ set(ft8_HEADERS
ft8.h
ft8plan.h
ft8plans.h
ft8stats.h
libldpc.h
osd.h
packing.h

418
ft8/ft8.cpp
View File

@ -30,7 +30,6 @@
#include <stdio.h>
// #include <assert.h>
#include <math.h>
#include <complex>
#include <fftw3.h>
#include <algorithm>
#include <complex>
@ -120,178 +119,6 @@ std::vector<float> blackmanharris(int n)
return h;
}
Stats::Stats(int how, float log_tail, float log_rate) :
sum_(0),
finalized_(false),
how_(how),
log_tail_(log_tail),
log_rate_(log_rate)
{}
void Stats::add(float x)
{
a_.push_back(x);
sum_ += x;
finalized_ = false;
}
void Stats::finalize()
{
finalized_ = true;
int n = a_.size();
mean_ = sum_ / n;
float var = 0;
float bsum = 0;
for (int i = 0; i < n; i++)
{
float y = a_[i] - mean_;
var += y * y;
bsum += fabs(y);
}
var /= n;
stddev_ = sqrt(var);
b_ = bsum / n;
// prepare for binary search to find where values lie
// in the distribution.
if (how_ != 0 && how_ != 5) {
std::sort(a_.begin(), a_.end());
}
}
float Stats::mean()
{
if (!finalized_) {
finalize();
}
return mean_;
}
float Stats::stddev()
{
if (!finalized_) {
finalize();
}
return stddev_;
}
// fraction of distribution that's less than x.
// assumes normal distribution.
// this is PHI(x), or the CDF at x,
// or the integral from -infinity
// to x of the PDF.
float Stats::gaussian_problt(float x)
{
float SDs = (x - mean()) / stddev();
float frac = 0.5 * (1.0 + erf(SDs / sqrt(2.0)));
return frac;
}
// https://en.wikipedia.org/wiki/Laplace_distribution
// m and b from page 116 of Mark Owen's Practical Signal Processing.
float Stats::laplace_problt(float x)
{
float m = mean();
float cdf;
if (x < m) {
cdf = 0.5 * exp((x - m) / b_);
} else {
cdf = 1.0 - 0.5 * exp(-(x - m) / b_);
}
return cdf;
}
// look into the actual distribution.
float Stats::problt(float x)
{
if (!finalized_) {
finalize();
}
if (how_ == 0) {
return gaussian_problt(x);
}
if (how_ == 5) {
return laplace_problt(x);
}
// binary search.
auto it = std::lower_bound(a_.begin(), a_.end(), x);
int i = it - a_.begin();
int n = a_.size();
if (how_ == 1)
{
// index into the distribution.
// works poorly for values that are off the ends
// of the distribution, since those are all
// mapped to 0.0 or 1.0, regardless of magnitude.
return i / (float)n;
}
if (how_ == 2)
{
// use a kind of logistic regression for
// values near the edges of the distribution.
if (i < log_tail_ * n)
{
float x0 = a_[(int)(log_tail_ * n)];
float y = 1.0 / (1.0 + exp(-log_rate_ * (x - x0)));
// y is 0..0.5
y /= 5;
return y;
}
else if (i > (1 - log_tail_) * n)
{
float x0 = a_[(int)((1 - log_tail_) * n)];
float y = 1.0 / (1.0 + exp(-log_rate_ * (x - x0)));
// y is 0.5..1
// we want (1-log_tail)..1
y -= 0.5;
y *= 2;
y *= log_tail_;
y += (1 - log_tail_);
return y;
}
else
{
return i / (float)n;
}
}
if (how_ == 3)
{
// gaussian for values near the edge of the distribution.
if (i < log_tail_ * n) {
return gaussian_problt(x);
} else if (i > (1 - log_tail_) * n) {
return gaussian_problt(x);
} else {
return i / (float)n;
}
}
if (how_ == 4)
{
// gaussian for values outside the distribution.
if (x < a_[0] || x > a_.back()) {
return gaussian_problt(x);
} else {
return i / (float)n;
}
}
return 0;
}
// a-priori probability of each of the 174 LDPC codeword
// bits being one. measured from reconstructed correct
// codewords, into ft8bits, then python bprob.py.
@ -1584,6 +1411,93 @@ std::vector<std::vector<std::complex<float>>> FT8::c_convert_to_snr(
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.
@ -1643,6 +1557,38 @@ void FT8::make_stats(
}
}
//
// 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.
//
@ -1767,6 +1713,7 @@ std::vector<std::vector<float>> FT8::soft_c2m(const FFTEngine::ffts_t &c79)
// returns log-likelihood, zero is positive, one is negative.
//
float FT8::bayes(
FT8Params& params,
float best_zero,
float best_one,
int lli,
@ -1799,6 +1746,7 @@ float FT8::bayes(
// 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));
@ -1806,6 +1754,7 @@ float FT8::bayes(
// 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));
@ -1819,6 +1768,8 @@ float FT8::bayes(
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 {
@ -1944,13 +1895,89 @@ void FT8::soft_decode(const FFTEngine::ffts_t &c79, float ll174[])
}
}
float ll = bayes(best_zero, best_one, lli, bests, all);
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[])
{
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);
int lli = 0;
int zoX = 1<<(nbSymbolBits-1);
int zoY = nbSymbolBits;
int *zeroi = new int[zoX*zoY];
int *onei = new int[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.
//
@ -2135,13 +2162,50 @@ void FT8::c_soft_decode(const FFTEngine::ffts_t &c79x, float ll174[])
}
}
float ll = bayes(best_zero, best_one, lli, bests, all);
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,
@ -2304,7 +2368,7 @@ void FT8::soft_decode_pairs(
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(best_zero, best_one, lli, bests, all);
float ll = bayes(params, best_zero, best_one, lli, bests, all);
ll174[lli++] = ll;
}
}
@ -2468,7 +2532,7 @@ void FT8::soft_decode_triples(
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(best_zero, best_one, lli, bests, all);
float ll = bayes(params, best_zero, best_one, lli, bests, all);
ll174[lli++] = ll;
}
}

84
ft8/ft8.h
View File

@ -30,6 +30,7 @@
#include <QString>
#include "fft.h"
#include "ft8stats.h"
#include "export.h"
class QThread;
@ -51,53 +52,6 @@ public:
virtual QString get_name() = 0;
};
//
// manage statistics for soft decoding, to help
// decide how likely each symbol is to be correct,
// to drive LDPC decoding.
//
// meaning of the how (problt_how) parameter:
// 0: gaussian
// 1: index into the actual distribution
// 2: do something complex for the tails.
// 3: index into the actual distribution plus gaussian for tails.
// 4: similar to 3.
// 5: laplace
//
class FT8_API Stats
{
public:
std::vector<float> a_;
float sum_;
bool finalized_;
float mean_; // cached
float stddev_; // cached
float b_; // cached
int how_;
public:
Stats(int how, float log_tail, float log_rate);
void add(float x);
void finalize();
float mean();
float stddev();
// fraction of distribution that's less than x.
// assumes normal distribution.
// this is PHI(x), or the CDF at x,
// or the integral from -infinity
// to x of the PDF.
float gaussian_problt(float x);
// https://en.wikipedia.org/wiki/Laplace_distribution
// m and b from page 116 of Mark Owen's Practical Signal Processing.
float laplace_problt(float x);
// look into the actual distribution.
float problt(float x);
private:
float log_tail_;
float log_rate_;
};
class FT8_API Strength
{
@ -220,8 +174,8 @@ struct FT8_API FT8Params
third_off_win = 0.075;
log_tail = 0.1;
log_rate = 8.0;
problt_how_noise = 0;
problt_how_sig = 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;
@ -320,6 +274,18 @@ public:
// append the 83 bits to the 91 bits messag e+ crc to obbain the 174 bit payload
static void encode(int a174[], 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(FT8Params& params, const std::vector<std::vector<float>>& mags, int nbSymbolBits, float ll174[]);
private:
//
// reduce the sample rate from arate to brate.
@ -444,6 +410,11 @@ private:
// normalize levels by windowed median.
// this helps, but why?
//
static std::vector<std::vector<float>> convert_to_snr_gen(const FT8Params& params, int nbSymbolBits, const std::vector<std::vector<float>> &mags);
//
// normalize levels by windowed median.
// this helps, but why?
//
std::vector<std::vector<std::complex<float>>> c_convert_to_snr(
const std::vector<std::vector<std::complex<float>>> &m79
);
@ -453,12 +424,22 @@ private:
// distribution of strongest tones, and
// distribution of noise.
//
void make_stats(
static void make_stats(
const std::vector<std::vector<float>> &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<std::vector<float>> &mags,
int nbSymbolBits,
Stats &bests,
Stats &all
);
//
// convert 79x8 complex FFT bins to magnitudes.
//
// exploits local phase coherence by decreasing magnitudes of bins
@ -477,7 +458,8 @@ private:
//
// returns log-likelihood, zero is positive, one is negative.
//
float bayes(
static float bayes(
FT8Params& params,
float best_zero,
float best_one,
int lli,

200
ft8/ft8stats.cpp Normal file
View File

@ -0,0 +1,200 @@
///////////////////////////////////////////////////////////////////////////////////
// Copyright (C) 2024 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/>. //
///////////////////////////////////////////////////////////////////////////////////
#include <math.h>
#include <algorithm>
#include "ft8stats.h"
namespace FT8 {
Stats::Stats(int how, float log_tail, float log_rate) :
sum_(0),
finalized_(false),
how_(how),
log_tail_(log_tail),
log_rate_(log_rate)
{}
void Stats::add(float x)
{
a_.push_back(x);
sum_ += x;
finalized_ = false;
}
void Stats::finalize()
{
finalized_ = true;
int n = a_.size();
mean_ = sum_ / n;
float var = 0;
float bsum = 0;
for (int i = 0; i < n; i++)
{
float y = a_[i] - mean_;
var += y * y;
bsum += fabs(y);
}
var /= n;
stddev_ = sqrt(var);
b_ = bsum / n;
// prepare for binary search to find where values lie
// in the distribution.
if (how_ != 0 && how_ != 5) {
std::sort(a_.begin(), a_.end());
}
}
float Stats::mean()
{
if (!finalized_) {
finalize();
}
return mean_;
}
float Stats::stddev()
{
if (!finalized_) {
finalize();
}
return stddev_;
}
// fraction of distribution that's less than x.
// assumes normal distribution.
// this is PHI(x), or the CDF at x,
// or the integral from -infinity
// to x of the PDF.
float Stats::gaussian_problt(float x)
{
float SDs = (x - mean()) / stddev();
float frac = 0.5 * (1.0 + erf(SDs / sqrt(2.0)));
return frac;
}
// https://en.wikipedia.org/wiki/Laplace_distribution
// m and b from page 116 of Mark Owen's Practical Signal Processing.
float Stats::laplace_problt(float x)
{
float m = mean();
float cdf;
if (x < m) {
cdf = 0.5 * exp((x - m) / b_);
} else {
cdf = 1.0 - 0.5 * exp(-(x - m) / b_);
}
return cdf;
}
// look into the actual distribution.
float Stats::problt(float x)
{
if (!finalized_) {
finalize();
}
if (how_ == 0) {
return gaussian_problt(x);
}
if (how_ == 5) {
return laplace_problt(x);
}
// binary search.
auto it = std::lower_bound(a_.begin(), a_.end(), x);
int i = it - a_.begin();
int n = a_.size();
if (how_ == 1)
{
// index into the distribution.
// works poorly for values that are off the ends
// of the distribution, since those are all
// mapped to 0.0 or 1.0, regardless of magnitude.
return i / (float)n;
}
if (how_ == 2)
{
// use a kind of logistic regression for
// values near the edges of the distribution.
if (i < log_tail_ * n)
{
float x0 = a_[(int)(log_tail_ * n)];
float y = 1.0 / (1.0 + exp(-log_rate_ * (x - x0)));
// y is 0..0.5
y /= 5;
return y;
}
else if (i > (1 - log_tail_) * n)
{
float x0 = a_[(int)((1 - log_tail_) * n)];
float y = 1.0 / (1.0 + exp(-log_rate_ * (x - x0)));
// y is 0.5..1
// we want (1-log_tail)..1
y -= 0.5;
y *= 2;
y *= log_tail_;
y += (1 - log_tail_);
return y;
}
else
{
return i / (float)n;
}
}
if (how_ == 3)
{
// gaussian for values near the edge of the distribution.
if (i < log_tail_ * n) {
return gaussian_problt(x);
} else if (i > (1 - log_tail_) * n) {
return gaussian_problt(x);
} else {
return i / (float)n;
}
}
if (how_ == 4)
{
// gaussian for values outside the distribution.
if (x < a_[0] || x > a_.back()) {
return gaussian_problt(x);
} else {
return i / (float)n;
}
}
return 0;
}
} // namespace FT8

79
ft8/ft8stats.h Normal file
View File

@ -0,0 +1,79 @@
///////////////////////////////////////////////////////////////////////////////////
// Copyright (C) 2024 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/>. //
///////////////////////////////////////////////////////////////////////////////////
#ifndef _ft8stats_h_
#define _ft8stats_h_
#include <vector>
#include "export.h"
namespace FT8 {
//
// manage statistics for soft decoding, to help
// decide how likely each symbol is to be correct,
// to drive LDPC decoding.
//
// meaning of the how (problt_how) parameter:
// 0: gaussian
// 1: index into the actual distribution
// 2: do something complex for the tails.
// 3: index into the actual distribution plus gaussian for tails.
// 4: similar to 3.
// 5: laplace
//
class FT8_API Stats
{
public:
std::vector<float> a_;
float sum_;
bool finalized_;
float mean_; // cached
float stddev_; // cached
float b_; // cached
int how_;
public:
Stats(int how, float log_tail, float log_rate);
void add(float x);
void finalize();
float mean();
float stddev();
// fraction of distribution that's less than x.
// assumes normal distribution.
// this is PHI(x), or the CDF at x,
// or the integral from -infinity
// to x of the PDF.
float gaussian_problt(float x);
// https://en.wikipedia.org/wiki/Laplace_distribution
// m and b from page 116 of Mark Owen's Practical Signal Processing.
float laplace_problt(float x);
// look into the actual distribution.
float problt(float x);
private:
float log_tail_;
float log_rate_;
};
} //namespace FT8
#endif // _ft8stats_h_

5
plugins/channelrx/demodchirpchat/chirpchatdemodgui.ui
View File

@ -585,6 +585,11 @@
<string>TTY</string>
</property>
</item>
<item>
<property name="text">
<string>FT</string>
</property>
</item>
</widget>
</item>
<item>

6
plugins/channelrx/demodchirpchat/chirpchatdemodmsg.h
View File

@ -28,6 +28,7 @@ namespace ChirpChatDemodMsg
public:
const std::vector<unsigned short>& getSymbols() const { return m_symbols; }
const std::vector<std::vector<float>>& getMagnitudes() const { return m_magnitudes; }
unsigned int getSyncWord() const { return m_syncWord; }
float getSingalDb() const { return m_signalDb; }
float getNoiseDb() const { return m_noiseDb; }
@ -48,6 +49,10 @@ namespace ChirpChatDemodMsg
m_noiseDb = db;
}
void pushBackMagnitudes(const std::vector<float>& magnitudes) {
m_magnitudes.push_back(magnitudes);
}
static MsgDecodeSymbols* create() {
return new MsgDecodeSymbols();
}
@ -57,6 +62,7 @@ namespace ChirpChatDemodMsg
private:
std::vector<unsigned short> m_symbols;
std::vector<std::vector<float>> m_magnitudes;
unsigned int m_syncWord;
float m_signalDb;
float m_noiseDb;

3
plugins/channelrx/demodchirpchat/chirpchatdemodsettings.h
View File

@ -37,7 +37,8 @@ struct ChirpChatDemodSettings
{
CodingLoRa, //!< Standard LoRa
CodingASCII, //!< plain ASCII (7 bits)
CodingTTY //!< plain TTY (5 bits)
CodingTTY, //!< plain TTY (5 bits)
CodingFT //!< FT8/4 scheme (payload 174 bits LDPC)
};
int m_inputFrequencyOffset;

153
plugins/channelrx/demodchirpchat/chirpchatdemodsink.cpp
View File

@ -382,18 +382,63 @@ void ChirpChatDemodSink::processSample(const Complex& ci)
m_fft->transform();
m_fftCounter = 0;
double magsq, magsqTotal;
unsigned short symbol;
unsigned short symbol = evalSymbol(
argmax(
m_fft->out(),
m_fftInterpolation,
m_fftLength,
magsq,
magsqTotal,
m_spectrumBuffer,
m_fftInterpolation
)
) % m_nbSymbolsEff;
if (m_settings.m_codingScheme == ChirpChatDemodSettings::CodingFT)
{
std::vector<float> magnitudes;
symbol = evalSymbol(
extractMagnitudes(
magnitudes,
m_fft->out(),
m_fftInterpolation,
m_fftLength,
magsq,
magsqTotal,
m_spectrumBuffer,
m_fftInterpolation
)
) % m_nbSymbolsEff;
m_decodeMsg->pushBackSymbol(symbol);
m_decodeMsg->pushBackMagnitudes(magnitudes);
}
else
{
int imax;
if (m_settings.m_deBits > 0)
{
double magSqNoise;
imax = argmaxSpreaded(
m_fft->out(),
m_fftInterpolation,
m_fftLength,
magsq,
magSqNoise,
magsqTotal,
m_spectrumBuffer,
m_fftInterpolation
);
// double dbS = CalcDb::dbPower(magsq);
// double dbN = CalcDb::dbPower(magSqNoise);
// qDebug("ChirpChatDemodSink::processSample: S: %5.2f N: %5.2f S/N: %5.2f", dbS, dbN, dbS - dbN);
}
else
{
imax = argmax(
m_fft->out(),
m_fftInterpolation,
m_fftLength,
magsq,
magsqTotal,
m_spectrumBuffer,
m_fftInterpolation
);
}
symbol = evalSymbol(imax) % m_nbSymbolsEff;
m_decodeMsg->pushBackSymbol(symbol);
}
if (m_spectrumSink) {
m_spectrumSink->feed(m_spectrumBuffer, m_nbSymbols);
@ -405,13 +450,18 @@ void ChirpChatDemodSink::processSample(const Complex& ci)
m_magsqTotalAvg(magsq);
m_decodeMsg->pushBackSymbol(symbol);
if ((m_chirpCount == 0)
|| (m_settings.m_eomSquelchTenths == 121) // max - disable squelch
|| ((m_settings.m_eomSquelchTenths*magsq)/10.0 > m_magsqMax))
{
qDebug("ChirpChatDemodSink::processSample: symbol %02u: %4u|%11.6f", m_chirpCount, symbol, magsq);
// const std::vector<float>& magnitudes = m_decodeMsg->getMagnitudes().back();
// int i = 0;
// for (auto magnitude : magnitudes)
// {
// qDebug("ChirpChatDemodSink::processSample: mag[%02d] = %11.6f", i, magnitude);
// i++;
// }
m_magsqOnAvg(magsq);
m_chirpCount++;
@ -508,38 +558,93 @@ unsigned int ChirpChatDemodSink::argmax(
return imax;
}
unsigned int ChirpChatDemodSink::extractMagnitudes(
std::vector<float>& magnitudes,
const Complex *fftBins,
unsigned int fftMult,
unsigned int fftLength,
double& magsqMax,
double& magsqTotal,
Complex *specBuffer,
unsigned int specDecim)
{
magsqMax = 0.0;
magsqTotal = 0.0;
unsigned int imax = 0;
double magSum = 0.0;
unsigned int spread = fftMult * (1<<m_settings.m_deBits);
unsigned int istart = fftMult*fftLength - spread/2 + 1;
float magnitude = 0.0;
for (unsigned int i2 = istart; i2 < istart + fftMult*fftLength; i2++)
{
int i = i2 % (fftMult*fftLength);
double magsq = std::norm(fftBins[i]);
magsqTotal += magsq;
magnitude += magsq;
if (i % spread == (spread/2)-1) // boundary (inclusive)
{
if (magnitude > magsqMax)
{
imax = (i/spread)*spread;
magsqMax = magnitude;
}
magnitudes.push_back(magnitude);
magnitude = 0.0;
}
if (specBuffer)
{
magSum += magsq;
if (i % specDecim == specDecim - 1)
{
specBuffer[i/specDecim] = Complex(std::polar(magSum, 0.0));
magSum = 0.0;
}
}
}
magsqTotal /= fftMult*fftLength;
return imax;
}
unsigned int ChirpChatDemodSink::argmaxSpreaded(
const Complex *fftBins,
unsigned int fftMult,
unsigned int fftLength,
double& magsqMax,
double& magsqNoise,
double& magSqTotal,
double& magsqTotal,
Complex *specBuffer,
unsigned int specDecim)
{
magsqMax = 0.0;
magsqNoise = 0.0;
magSqTotal = 0.0;
magsqTotal = 0.0;
unsigned int imax = 0;
double magSum = 0.0;
double magSymbol = 0.0;
unsigned int nbsymbols = 1<<(m_settings.m_spreadFactor - m_settings.m_deBits);
unsigned int spread = fftMult * (1<<m_settings.m_deBits);
unsigned int istart = fftMult*fftLength - spread/2 + 1;
double magSymbol = 0.0;
for (unsigned int i2 = istart; i2 < istart + fftMult*fftLength; i2++)
{
unsigned int i = i2 % (fftMult*fftLength);
int i = i2 % (fftMult*fftLength);
double magsq = std::norm(fftBins[i]);
magsqTotal += magsq;
magSymbol += magsq;
magSqTotal += magsq;
if (i % spread == spread/2) // boundary (inclusive)
if (i % spread == (spread/2)-1) // boundary (inclusive)
{
if (magSymbol > magsqMax)
{
imax = (i/spread)*spread;
magsqMax = magSymbol;
imax = i;
}
magsqNoise += magSymbol;
@ -559,10 +664,12 @@ unsigned int ChirpChatDemodSink::argmaxSpreaded(
}
magsqNoise -= magsqMax;
magsqNoise /= fftLength;
magSqTotal /= fftMult*fftLength;
magsqNoise /= (nbsymbols - 1);
magsqTotal /= nbsymbols;
// magsqNoise /= fftLength;
// magsqTotal /= fftMult*fftLength;
return imax / spread;
return imax;
}
void ChirpChatDemodSink::decimateSpectrum(Complex *in, Complex *out, unsigned int size, unsigned int decimation)

10
plugins/channelrx/demodchirpchat/chirpchatdemodsink.h
View File

@ -138,6 +138,16 @@ private:
Complex *specBuffer,
unsigned int specDecim
);
unsigned int extractMagnitudes(
std::vector<float>& magnitudes,
const Complex *fftBins,
unsigned int fftMult,
unsigned int fftLength,
double& magsqMax,
double& magSqTotal,
Complex *specBuffer,
unsigned int specDecim
);
void decimateSpectrum(Complex *in, Complex *out, unsigned int size, unsigned int decimation);
int toSigned(int u, int intSize);
unsigned int evalSymbol(unsigned int rawSymbol);

151
sdrbench/test_ft8protocols.cpp
View File

@ -18,6 +18,7 @@
#include <iostream>
#include <fstream>
#include <regex>
#include <random>
#include <QTextStream>
@ -40,10 +41,12 @@ class TestFT8Protocols
public:
static void testMsg1(const QStringList& argElements, bool runLDPC = false);
static void testMsg00(const QStringList& argElements, bool runLDPC = false);
static void testOnesZeroes(const QStringList& argElements);
static void testSoftDecode(const QStringList& argElements);
private:
static bool testLDPC(int a77[]);
static bool compare174(int a174[], int r174[]);
static bool compareBits(int a[], int r[], int nbBits = 174);
static void debugIntArray(int a[], int length);
};
@ -68,6 +71,10 @@ void MainBench::testFT8Protocols(const QString& argsStr)
TestFT8Protocols::testMsg1(argElements, true); // type 1 message test with LDPC encoding/decoding test
} else if (testType == "msg00L") {
TestFT8Protocols::testMsg00(argElements, true); // type 0.0 message test with LDPC encoding/decoding test
} else if (testType == "zeroones") {
TestFT8Protocols::testOnesZeroes(argElements);
} else if (testType == "softdec") {
TestFT8Protocols::testSoftDecode(argElements);
} else {
qWarning("MainBench::testFT8Protocols: unrecognized test type");
}
@ -192,17 +199,17 @@ bool TestFT8Protocols::testLDPC(int a77[])
}
else
{
return compare174(a174, r174);
return compareBits(a174, r174);
}
}
bool TestFT8Protocols::compare174(int a174[], int r174[])
bool TestFT8Protocols::compareBits(int a[], int r[], int nbBits)
{
for (int i=0; i < 174; i++)
for (int i=0; i < nbBits; i++)
{
if (a174[i] != r174[i])
if (a[i] != r[i])
{
qDebug("TestFT8Protocols::compare174: failed at index %d: %d != %d", i, a174[i], r174[i]);
qDebug("TestFT8Protocols::compareBits: failed at index %d: %d != %d", i, a[i], r[i]);
return false;
}
}
@ -222,4 +229,136 @@ void TestFT8Protocols::debugIntArray(int a[], int length)
qDebug("TestFT8Protocols::debugIntArray: %s", qPrintable(s));
}
void TestFT8Protocols::testOnesZeroes(const QStringList& argElements)
{
if (argElements.size() < 3)
{
qWarning("TestFT8Protocols::testOnesZeroes: not enough elements");
return;
}
int nbBits, bitIndex;
bool intOK;
nbBits = argElements[1].toInt(&intOK);
if (!intOK)
{
qWarning("TestFT8Protocols::testOnesZeroes: first argument is not numeric: %s", qPrintable(argElements[1]));
return;
}
bitIndex = argElements[2].toInt(&intOK);
if (!intOK)
{
qWarning("TestFT8Protocols::testOnesZeroes: second argument is not numeric: %s", qPrintable(argElements[2]));
return;
}
if (nbBits < 2)
{
qWarning("TestFT8Protocols::testOnesZeroes: nbBits too small: %d", nbBits);
return;
}
bitIndex = bitIndex > nbBits - 1 ? nbBits - 1 : bitIndex;
int *ones = new int[1<<nbBits];
int *zeroes = new int[1<<nbBits];
FT8::FT8::set_ones_zeroes(ones, zeroes, nbBits, bitIndex);
QString s;
QTextStream os(&s);
for (int i = 0; i < (1<<(nbBits-1)); i++) {
os << i << ": " << zeroes[i] << ", " << ones[i] << "\n";
}
qInfo("TestFT8Protocols::testOnesZeroes: (%d,%d) index: zeroes, ones:\n%s", nbBits, bitIndex, qPrintable(s));
}
void TestFT8Protocols::testSoftDecode(const QStringList& argElements)
{
if (argElements.size() < 3)
{
qWarning("TestFT8Protocols::testSoftDecode: not enough elements");
return;
}
bool intOK;
int nbBits = argElements[1].toInt(&intOK);
if (!intOK)
{
qWarning("TestFT8Protocols::testSoftDecode: first argument is not numeric: %s", qPrintable(argElements[1]));
return;
}
if ((nbBits < 2) || (nbBits > 12))
{
qWarning("TestFT8Protocols::testSoftDecode: bits peer symbols invalid: %d", nbBits);
return;
}
int symbolSize = 1<<nbBits;
std::vector<float> magSymbols(symbolSize);
std::vector<std::vector<float>> mags;
std::random_device rd;
std::mt19937 gen(rd());
std::uniform_real_distribution<float> dist(0.0, 0.01);
for (int i = 2; i < argElements.size(); i++)
{
int symbol = argElements[i].toInt(&intOK);
if (!intOK)
{
qWarning("TestFT8Protocols::testSoftDecode: symbol is not numeric: %s", qPrintable(argElements[i]));
return;
}
for (auto& m : magSymbols) {
m = 0.01 + dist(gen);
}
symbol = symbol % symbolSize;
magSymbols[symbol] += 0.01;
mags.push_back(magSymbols);
}
QString s;
QTextStream os(&s);
qDebug("TestFT8Protocols::testSoftDecode: mags:");
for (const auto& magrow : mags)
{
for (const auto& mag : magrow) {
os << mag << " ";
}
qDebug("TestFT8Protocols::testSoftDecode: %s", qPrintable(s));
s.clear();
}
float *lls = new float[mags.size()*nbBits];
std::fill(lls, lls+mags.size()*nbBits, 0.0);
FT8::FT8Params params;
FT8::FT8::soft_decode_mags(params, mags, nbBits, lls);
for (unsigned int si = 0; si < mags.size(); si++)
{
for (int biti = 0; biti < nbBits; biti++) {
os << " " << lls[nbBits*si + biti];
}
os << " ";
}
// for (unsigned int i = 0; i < mags.size()*nbBits; i++) {
// os << " " << lls[i];
// }
qInfo("TestFT8Protocols::testSoftDecode: lls: %s", qPrintable(s));
delete[] lls;
}
#endif