From 8e77ad01adad67f9129ee9494d28c37b4b89ed1b Mon Sep 17 00:00:00 2001 From: f4exb Date: Wed, 3 Apr 2024 12:54:03 +0200 Subject: [PATCH] FT8: added soft decode on magnitudes and genaralize to any number of bits per symbol --- ft8/CMakeLists.txt | 2 + ft8/ft8.cpp | 418 ++++++++++-------- ft8/ft8.h | 84 ++-- ft8/ft8stats.cpp | 200 +++++++++ ft8/ft8stats.h | 79 ++++ .../demodchirpchat/chirpchatdemodgui.ui | 5 + .../demodchirpchat/chirpchatdemodmsg.h | 6 + .../demodchirpchat/chirpchatdemodsettings.h | 3 +- .../demodchirpchat/chirpchatdemodsink.cpp | 153 ++++++- .../demodchirpchat/chirpchatdemodsink.h | 10 + sdrbench/test_ft8protocols.cpp | 151 ++++++- 11 files changed, 853 insertions(+), 258 deletions(-) create mode 100644 ft8/ft8stats.cpp create mode 100644 ft8/ft8stats.h diff --git a/ft8/CMakeLists.txt b/ft8/CMakeLists.txt index 75d62d0f7..b0060df7f 100644 --- a/ft8/CMakeLists.txt +++ b/ft8/CMakeLists.txt @@ -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 diff --git a/ft8/ft8.cpp b/ft8/ft8.cpp index cb9046a6e..617f97bf3 100644 --- a/ft8/ft8.cpp +++ b/ft8/ft8.cpp @@ -30,7 +30,6 @@ #include // #include #include -#include #include #include #include @@ -120,178 +119,6 @@ std::vector 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>> FT8::c_convert_to_snr( return n79; } +std::vector> FT8::convert_to_snr_gen(const FT8Params& params, int nbSymbolBits, const std::vector> &mags) +{ + if (params.snr_how < 0 || params.snr_win < 0) { + return mags; + } + + // + // for each symbol time, what's its "noise" level? + // + std::vector mm(mags.size()); + int nbSymbols = 1< v(nbSymbols); + float sum = 0.0; + + for (int bini = 0; bini < nbSymbols; bini++) + { + float x = mags[si][bini]; + v[bini] = x; + sum += x; + } + + if (params.snr_how != 1) { + std::sort(v.begin(), v.end()); + } + + int mid = nbSymbols / 2; + + if (params.snr_how == 0) { + // median + mm[si] = (v[mid-1] + v[mid]) / 2; + } else if (params.snr_how == 1) { + mm[si] = sum / nbSymbols; + } else if (params.snr_how == 2) { + // all but strongest tone. + mm[si] = std::accumulate(v.begin(), v.end() - 1, 0.0f) / (v.size() - 1); + } else if (params.snr_how == 3) { + mm[si] = v.front(); // weakest tone + } else if (params.snr_how == 4) { + mm[si] = v.back(); // strongest tone + } else if (params.snr_how == 5) { + mm[si] = v[v.size()-2]; // second-strongest tone + } else { + mm[si] = 1.0; + } + } + + // we're going to take a windowed average. + std::vector winwin; + + if (params.snr_win > 0) { + winwin = blackman(2 * params.snr_win + 1); + } else { + winwin.push_back(1.0); + } + + std::vector> snr(mags.size()); + + for (int si = 0; si < (int) mags.size(); si++) + { + float sum = 0; + + for (int dd = si - params.snr_win; dd <= si + params.snr_win; dd++) + { + int wi = dd - (si - params.snr_win); + + if (dd >= 0 && dd < (int) mags.size()) { + sum += mm[dd] * winwin[wi]; + } else if (dd < 0) { + sum += mm[0] * winwin[wi]; + } else { + sum += mm[mags.size()-1] * winwin[wi]; + } + } + + snr[si].resize(nbSymbols); + + for (int bi = 0; bi < nbSymbols; bi++) { + snr[si][bi] = mags[si][bi] / sum; + } + } + + return snr; +} + // // statistics to decide soft probabilities, // to drive LDPC decoder. @@ -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> &mags, + int nbSymbolBits, + Stats &bests, + Stats &all +) +{ + int nbBins = 1< mx) { + mx = x; + } + + all.add(x); + } + + bests.add(mx); + } +} + // // convert 79x8 complex FFT bins to magnitudes. // @@ -1767,6 +1713,7 @@ std::vector> 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>& mags_, int nbSymbolBits, float ll174[]) +{ + std::vector> mags = convert_to_snr_gen(params, nbSymbolBits, mags_); + // statistics to decide soft probabilities. + // distribution of strongest tones, and + // distribution of noise. + Stats bests(params.problt_how_sig, params.log_tail, params.log_rate); + Stats all(params.problt_how_noise, params.log_tail, params.log_rate); + make_stats_gen(mags, nbSymbolBits, bests, all); + 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< 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; } } diff --git a/ft8/ft8.h b/ft8/ft8.h index 41b3cdf7c..004e411ea 100644 --- a/ft8/ft8.h +++ b/ft8/ft8.h @@ -30,6 +30,7 @@ #include #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 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>& 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> convert_to_snr_gen(const FT8Params& params, int nbSymbolBits, const std::vector> &mags); + // + // normalize levels by windowed median. + // this helps, but why? + // std::vector>> c_convert_to_snr( const std::vector>> &m79 ); @@ -453,12 +424,22 @@ private: // distribution of strongest tones, and // distribution of noise. // - void make_stats( + static void make_stats( const std::vector> &m79, Stats &bests, Stats &all ); // + // generalized version of the above for any number of symbols and no Costas + // used by FT-chirp decoder + // + static void make_stats_gen( + const std::vector> &mags, + int nbSymbolBits, + Stats &bests, + Stats &all + ); + // // convert 79x8 complex FFT bins to magnitudes. // // exploits local phase coherence by decreasing magnitudes of bins @@ -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, diff --git a/ft8/ft8stats.cpp b/ft8/ft8stats.cpp new file mode 100644 index 000000000..8fc17181d --- /dev/null +++ b/ft8/ft8stats.cpp @@ -0,0 +1,200 @@ +/////////////////////////////////////////////////////////////////////////////////// +// Copyright (C) 2024 Edouard Griffiths, F4EXB // +// // +// 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 . // +/////////////////////////////////////////////////////////////////////////////////// + +#include +#include + +#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 diff --git a/ft8/ft8stats.h b/ft8/ft8stats.h new file mode 100644 index 000000000..dff8dc2f9 --- /dev/null +++ b/ft8/ft8stats.h @@ -0,0 +1,79 @@ +/////////////////////////////////////////////////////////////////////////////////// +// Copyright (C) 2024 Edouard Griffiths, F4EXB // +// // +// 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 . // +/////////////////////////////////////////////////////////////////////////////////// + +#ifndef _ft8stats_h_ +#define _ft8stats_h_ + +#include + +#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 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_ diff --git a/plugins/channelrx/demodchirpchat/chirpchatdemodgui.ui b/plugins/channelrx/demodchirpchat/chirpchatdemodgui.ui index d3e938fcc..39d5ba3d3 100644 --- a/plugins/channelrx/demodchirpchat/chirpchatdemodgui.ui +++ b/plugins/channelrx/demodchirpchat/chirpchatdemodgui.ui @@ -585,6 +585,11 @@ TTY + + + FT + + diff --git a/plugins/channelrx/demodchirpchat/chirpchatdemodmsg.h b/plugins/channelrx/demodchirpchat/chirpchatdemodmsg.h index 801daa848..58ee3da20 100644 --- a/plugins/channelrx/demodchirpchat/chirpchatdemodmsg.h +++ b/plugins/channelrx/demodchirpchat/chirpchatdemodmsg.h @@ -28,6 +28,7 @@ namespace ChirpChatDemodMsg public: const std::vector& getSymbols() const { return m_symbols; } + const std::vector>& 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& magnitudes) { + m_magnitudes.push_back(magnitudes); + } + static MsgDecodeSymbols* create() { return new MsgDecodeSymbols(); } @@ -57,6 +62,7 @@ namespace ChirpChatDemodMsg private: std::vector m_symbols; + std::vector> m_magnitudes; unsigned int m_syncWord; float m_signalDb; float m_noiseDb; diff --git a/plugins/channelrx/demodchirpchat/chirpchatdemodsettings.h b/plugins/channelrx/demodchirpchat/chirpchatdemodsettings.h index aef215fbc..051312d9c 100644 --- a/plugins/channelrx/demodchirpchat/chirpchatdemodsettings.h +++ b/plugins/channelrx/demodchirpchat/chirpchatdemodsettings.h @@ -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; diff --git a/plugins/channelrx/demodchirpchat/chirpchatdemodsink.cpp b/plugins/channelrx/demodchirpchat/chirpchatdemodsink.cpp index e7de652fd..3fe3f8fd2 100644 --- a/plugins/channelrx/demodchirpchat/chirpchatdemodsink.cpp +++ b/plugins/channelrx/demodchirpchat/chirpchatdemodsink.cpp @@ -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 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& 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& 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< 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< 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) diff --git a/plugins/channelrx/demodchirpchat/chirpchatdemodsink.h b/plugins/channelrx/demodchirpchat/chirpchatdemodsink.h index e670d3436..c8313f802 100644 --- a/plugins/channelrx/demodchirpchat/chirpchatdemodsink.h +++ b/plugins/channelrx/demodchirpchat/chirpchatdemodsink.h @@ -138,6 +138,16 @@ private: Complex *specBuffer, unsigned int specDecim ); + unsigned int extractMagnitudes( + std::vector& 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); diff --git a/sdrbench/test_ft8protocols.cpp b/sdrbench/test_ft8protocols.cpp index 1754deec0..b5cdc4d41 100644 --- a/sdrbench/test_ft8protocols.cpp +++ b/sdrbench/test_ft8protocols.cpp @@ -18,6 +18,7 @@ #include #include #include +#include #include @@ -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< 12)) + { + qWarning("TestFT8Protocols::testSoftDecode: bits peer symbols invalid: %d", nbBits); + return; + } + + int symbolSize = 1< magSymbols(symbolSize); + std::vector> mags; + std::random_device rd; + std::mt19937 gen(rd()); + std::uniform_real_distribution 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