Merge branch 'f4exb:master' into freq_scanner

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

ft8/ft8.cpp

@@ -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;
         }
     }

ft8/ft8.h

@@ -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_noise = 0;    // Gaussian
-        problt_how_sig = 0;
+        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,

ft8/ft8stats.cpp

@@ -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

ft8/ft8stats.h

@@ -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_

plugins/channelrx/demodchirpchat/chirpchatdemodgui.ui

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

plugins/channelrx/demodchirpchat/chirpchatdemodmsg.h

@@ -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;

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;

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(
+            if (m_settings.m_codingScheme == ChirpChatDemodSettings::CodingFT)
-                argmax(
+            {
-                    m_fft->out(),
+                std::vector<float> magnitudes;
-                    m_fftInterpolation,
+                symbol = evalSymbol(
-                    m_fftLength,
+                    extractMagnitudes(
-                    magsq,
+                        magnitudes,
-                    magsqTotal,
+                        m_fft->out(),
-                    m_spectrumBuffer,
+                        m_fftInterpolation,
-                    m_fftInterpolation
+                        m_fftLength,
-                )
+                        magsq,
-            ) % m_nbSymbolsEff;
+                        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;
+    magsqNoise /= (nbsymbols - 1);
-    magSqTotal /= fftMult*fftLength;
+    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)

plugins/channelrx/demodchirpchat/chirpchatdemodsink.h

@@ -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);

sdrbench/test_ft8protocols.cpp

@@ -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

sdrgui/mainspectrum/readme.md

@@ -22,7 +22,7 @@ The tooltip shows the device type, sequence number and serial number of the devi
 
 You may click on this area and drag the window with the mouse.
 
-<h3>3: Title</h3>
+<h3>2: Title</h3>
 
 The window title shows the device type and a sequence number for which the spectrum is displayed. It is the same as in the corresponding device window.