ham-ref/js8call
1
0
Fork 0
mirror of https://github.com/JS8Call-improved/JS8Call-improved synced 2026-08-13 17:47:36 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions
js8call/JS8.cpp
Allan Bazinet fddae761ab CMake cleanup and modernization 
Allan:
A gutting and modernization of the cmake configuration. This first pass is only expected to work on OSX, and will require additional work for Windows and Linux.
Fixes, formatting, and commentary.
Kill off doc_dir and data_dir; unused. These properties are legacy, and no longer used in current code.
Eliminate unused accessors

Chris:
Fix MACOSX_BUNDLE_ICON_FILE variable for Info.plist
Fix MACOSX_BUNDLE_SHORT_VERSION_STRING for Info.plist
Remove MACOSX_BUNDLE_LONG_VERSION_STRING which is no longer used
Clean up legacy Darwin files for building .dmg

Allan:
Add Qt Network library to find_library() call
Stub in (untested) windeployqt phase
Fix windeployqt pathing
Update Windows section; make OmniRig optional
Use std::numbers::pi instead of M_PI
Use std::numbers::pi in lieu of M_PI
Properly guard optional OmnRig include
Update CRC.h to address compiler pragma warnings
Additional OmniRig guarding
Simplify CMAKE_MODULE_PATH addition
Use list functionality, and append to the list, rather than prepend.

Chris:
Add Entitlements.plist for code signing with hardened runtime on MacOS
Add package directory for windeployqt to package the plugins in the proper directory structure for linking
add Qt Creator CMakeLists.txt.user  and build folder to .gitignore
move minimum required Qt version to 6.5 - qt_add_resources() only supported in Qt6.5 or newer
2025-10-13 18:18:55 -05:00

2800 lines
111 KiB
C++
Raw Blame History

This file contains ambiguous Unicode characters

This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.

/**
* (C) 2025 Allan Bazinet <w6baz@arrl.net> - All Rights Reserved
**/
#include "JS8.hpp"
#include <algorithm>
#include <atomic>
#include <cmath>
#include <complex>
#include <concepts>
#include <cstddef>
#include <cstdint>
#include <initializer_list>
#include <limits>
#include <memory>
#include <mutex>
#include <numbers>
#include <numeric>
#include <stdexcept>
#include <string_view>
#include <unordered_map>
#include <utility>
#include <vector>
#include <boost/crc.hpp>
#include <boost/math/ccmath/round.hpp>
#include <boost/multi_index_container.hpp>
#include <boost/multi_index/key.hpp>
#include <boost/multi_index/ordered_index.hpp>
#include <boost/multi_index/ranked_index.hpp>
#include <fftw3.h>
#include <vendor/Eigen/Dense>
#include <QDebug>
#include "commons.h"
// A C++ conversion of the Fortran JS8 encoding and decoder function.
// Some notes on the conversion:
//
// 1. Names of variables and functions as much as possible match those
// of the Fortran routines, for ease in cross-referencing during the
// debug comparison phase of testing. You don't have to like them; I
// don't like them either, frankly, but it's the reasonable approach
// to the problem as of this writing; we can make 'em pretty later.
//
// 2. The BP decoder should be a faithful reproduction of the Fortran
// version, albeit modified for the column-major vs. row-major
// differences between the two languages.
//
// 3. The OSD decoder is no longer used, and the depth is now fixed at
// 2, instead of being variable 1 to 4.
//
// 4. The Fortran version didn't compute the 40% rank consistently in
// syncjs8(); this version does. It wasn't typically off by much, but
// it was reliably not going to be at 40%. Hopefully, this change will
// result in more predictable first-pass candidate selection.
//
// 5. The Fortran version was very subject to Runge's phenomenon when
// computing the baseline in baselinejs8(), and was using a ton of
// data points below the 10% threshold for the polynomial determination.
// Neither of these seemed to be helpful, so in contrast we're using
// a number of Chebyshev nodes proportional to the desired polynomial
// terms.
//
// 6. The Fortran version normalized `s1` by dividing by the median in
// js8dec(), but did so in a naive manner, not checking for a median
// of zero. Testing indicates that the normalization does not appear
// to contribute to decoder yield, so it's been removed.
//
// 7. Translating array indices from the world of Fortran to that of C++
// is no one's fun task. If you see things that aren't behaving as
// expected, look at the Fortran code and compare the array indexing;
// would not be surprised in the least to have off-by-one errors here.
/******************************************************************************/
// Compilation Utilities
/******************************************************************************/
namespace
{
// Full-range cosine function using symmetries of cos(x). std::cos
// isn't constexpr until C++20, and we're targeting C++17 at the
// moment. We only use this function during compilation; std::cos
// is the better choice at runtime. Once we move to requiring a
// C++20 compiler, we can just use std::cos.
constexpr auto
cos(double x)
{
constexpr auto RAD_360 = std::numbers::pi * 2;
constexpr auto RAD_180 = std::numbers::pi;
constexpr auto RAD_90 = std::numbers::pi / 2;
// Polynomial approximation of cos(x) for x in [0, RAD_90],
// Accuracy here in theory is 1e-18, but double precision
// itself is only 1-e16, so within the domain of doubles,
// this should be extremely accurate.
constexpr auto cos = [](double const x)
{
constexpr std::array coefficients =
{
1.0, // Coefficient for x^0
-0.49999999999999994, // Coefficient for x^2
0.041666666666666664, // Coefficient for x^4
-0.001388888888888889, // Coefficient for x^6
0.000024801587301587, // Coefficient for x^8
-0.00000027557319223986, // Coefficient for x^10
0.00000000208767569878681, // Coefficient for x^12
-0.00000000001147074513875176, // Coefficient for x^14
0.0000000000000477947733238733 // Coefficient for x^16
};
auto const x2 = x * x;
auto const x4 = x2 * x2;
auto const x6 = x4 * x2;
auto const x8 = x4 * x4;
auto const x10 = x8 * x2;
auto const x12 = x8 * x4;
auto const x14 = x12 * x2;
auto const x16 = x8 * x8;
return coefficients[0]
+ coefficients[1] * x2
+ coefficients[2] * x4
+ coefficients[3] * x6
+ coefficients[4] * x8
+ coefficients[5] * x10
+ coefficients[6] * x12
+ coefficients[7] * x14
+ coefficients[8] * x16;
};
// Reduce x to [0, RAD_360)
x -= static_cast<long long>(x / RAD_360) * RAD_360;
// Map x to [0, RAD_180]
if (x > RAD_180) x = RAD_360 - x;
// Map x to [0, RAD_90] and evaluate the polynomial;
// flip the sign for angles in the second quadrant.
return x > RAD_90 ? -cos(RAD_180 - x) : cos(x);
};
}
/******************************************************************************/
// Constants
/******************************************************************************/
namespace
{
/* COMMON PARAMETERS */
// !Common
//
// parameter (KK=87) !Information bits (75 + CRC12)
// parameter (ND=58) !Data symbols
// parameter (NS=21) !Sync symbols (3 @ Costas 7x7)
// parameter (NN=NS+ND) !Total channel symbols (79)
// parameter (ASYNCMIN=1.5) !Minimum Sync
// parameter (NFSRCH=5) !Search frequency range in Hz (i.e., +/- 2.5 Hz)
// parameter (NMAXCAND=300) !Maximum number of candidate signals
// Parameter Value Description
// KK 87 Number of information bits (75 message bits + 12 CRC bits).
// ND 58 Number of data symbols in the JS8 transmission.
// NS 21 Number of synchronization symbols (3 Costas arrays of size 7).
// NN 79 Total number of channel symbols (NN = NS + ND).
// ASYNCMIN 1.5 Minimum sync value for successful decoding.
// NFSRCH 5 Search frequency range in Hz (±2.5 Hz).
// NMAXCAND 300 Maximum number of candidate signals.
constexpr int N = 174; // Total bits
constexpr int K = 87; // Message bits
constexpr int M = N - K; // Check bits
constexpr int KK = 87; // Information bits (75 + CRC12)
constexpr int ND = 58; // Data symbols
constexpr int NS = 21; // Sync symbols (3 @ Costas 7x7)
constexpr int NN = NS + ND; // Total channel symbols (79)
constexpr float ASYNCMIN = 1.5f; // Minimum sync
constexpr int NFSRCH = 5; // Search frequency range in Hz (i.e., +/- 2.5 Hz)
constexpr std::size_t NMAXCAND = 300; // Maxiumum number of candidate signals
constexpr int NFILT = 1400; // Filter length
constexpr int NROWS = 8;
constexpr int NFOS = 2;
constexpr int NSSY = 4;
constexpr int NP = 3200;
constexpr int NP2 = 2812;
constexpr float TAU = 2.0f * std::numbers::pi_v<float>;
constexpr auto ZERO = std::complex<float>{0.0f, 0.0f};
// Key for the constants that follow:
//
// NSUBMODE - ID of the submode
// NCOSTAS - Which JS8 Costas Arrays to use
// NSPS - Number of samples per second
// NTXDUR - Duration of the transmission in seconds.
// NDOWNSPS - Number of samples per symbol after downsampling.
// NDD - Parameter used in waveform tapering and related calculations. XXX
// JZ - Range of symbol offsets considered during decoding.
// ASTART - Start delay in seconds for decoding.
// BASESUB - XXX
// NMAX - Samples in input wave
// NSTEP - Rough time-sync step size
// NHSYM - Number of symbol spectra (1/4-sym steps)
// NDOW - Downsample factor to 32 samples per symbol
// NQSYMBOL - Downsample factor of a quarter symbol
/* A MODE DECODER */
struct ModeA
{
// Static constants
inline static constexpr int NSUBMODE = 0;
inline static constexpr auto NCOSTAS = JS8::Costas::Type::ORIGINAL;
inline static constexpr int NSPS = JS8A_SYMBOL_SAMPLES;
inline static constexpr int NTXDUR = JS8A_TX_SECONDS;
inline static constexpr int NDOWNSPS = 32;
inline static constexpr int NDD = 100;
inline static constexpr int JZ = 62;
inline static constexpr float ASTART = 0.5f;
inline static constexpr float BASESUB = 40.0f;
// Derived parameters
inline static constexpr float AZ = (12000.0f / NSPS) * 0.64f;
inline static constexpr int NMAX = NTXDUR * JS8_RX_SAMPLE_RATE;
inline static constexpr int NFFT1 = NSPS * NFOS;
inline static constexpr int NSTEP = NSPS / NSSY;
inline static constexpr int NHSYM = NMAX / NSTEP - 3;
inline static constexpr int NDOWN = NSPS / NDOWNSPS;
inline static constexpr int NQSYMBOL = NDOWNSPS / 4;
inline static constexpr int NDFFT1 = NSPS * NDD;
inline static constexpr int NDFFT2 = NDFFT1 / NDOWN;
inline static constexpr int NP2 = NN * NDOWNSPS;
inline static constexpr float TSTEP = NSTEP / 12000.0f;
inline static constexpr int JSTRT = ASTART / TSTEP;
inline static constexpr float DF = 12000.0f / NFFT1;
};
/* B MODE DECODER */
struct ModeB
{
// Static constants
inline static constexpr int NSUBMODE = 1;
inline static constexpr auto NCOSTAS = JS8::Costas::Type::MODIFIED;
inline static constexpr int NSPS = JS8B_SYMBOL_SAMPLES;
inline static constexpr int NTXDUR = JS8B_TX_SECONDS;
inline static constexpr int NDOWNSPS = 20;
inline static constexpr int NDD = 100;
inline static constexpr int JZ = 144;
inline static constexpr float ASTART = 0.2f;
inline static constexpr float BASESUB = 39.0f;
// Derived parameters
inline static constexpr float AZ = (12000.0f / NSPS) * 0.8f;
inline static constexpr int NMAX = NTXDUR * JS8_RX_SAMPLE_RATE;
inline static constexpr int NFFT1 = NSPS * NFOS;
inline static constexpr int NSTEP = NSPS / NSSY;
inline static constexpr int NHSYM = NMAX / NSTEP - 3;
inline static constexpr int NDOWN = NSPS / NDOWNSPS;
inline static constexpr int NQSYMBOL = NDOWNSPS / 4;
inline static constexpr int NDFFT1 = NSPS * NDD;
inline static constexpr int NDFFT2 = NDFFT1 / NDOWN;
inline static constexpr int NP2 = NN * NDOWNSPS;
inline static constexpr float TSTEP = NSTEP / 12000.0f;
inline static constexpr int JSTRT = ASTART / TSTEP;
inline static constexpr float DF = 12000.0f / NFFT1;
};
/* C MODE DECODER */
struct ModeC
{
// Static constants
inline static constexpr int NSUBMODE = 2;
inline static constexpr auto NCOSTAS = JS8::Costas::Type::MODIFIED;
inline static constexpr int NSPS = JS8C_SYMBOL_SAMPLES;
inline static constexpr int NTXDUR = JS8C_TX_SECONDS;
inline static constexpr int NDOWNSPS = 12;
inline static constexpr int NDD = 120;
inline static constexpr int JZ = 172;
inline static constexpr float ASTART = 0.1f;
inline static constexpr float BASESUB = 38.0f;
// Derived parameters
inline static constexpr float AZ = (12000.0f / NSPS) * 0.6f;
inline static constexpr int NMAX = NTXDUR * JS8_RX_SAMPLE_RATE;
inline static constexpr int NFFT1 = NSPS * NFOS;
inline static constexpr int NSTEP = NSPS / NSSY;
inline static constexpr int NHSYM = NMAX / NSTEP - 3;
inline static constexpr int NDOWN = NSPS / NDOWNSPS;
inline static constexpr int NQSYMBOL = NDOWNSPS / 4;
inline static constexpr int NDFFT1 = NSPS * NDD;
inline static constexpr int NDFFT2 = NDFFT1 / NDOWN;
inline static constexpr int NP2 = NN * NDOWNSPS;
inline static constexpr float TSTEP = NSTEP / 12000.0f;
inline static constexpr int JSTRT = ASTART / TSTEP;
inline static constexpr float DF = 12000.0f / NFFT1;
};
/* E MODE DECODER */
// Note that the original used 28 for NTXDUR and 90 for NDD, but the
// corresponding C++ mainline side used 30 for NTXDUR, so for the
// moment, we're matching that here, which seems logical at present.
struct ModeE
{
// Static constants
inline static constexpr int NSUBMODE = 4;
inline static constexpr auto NCOSTAS = JS8::Costas::Type::MODIFIED;
inline static constexpr int NSPS = JS8E_SYMBOL_SAMPLES;
inline static constexpr int NTXDUR = JS8E_TX_SECONDS; // XXX was 28 in Fortran
inline static constexpr int NDOWNSPS = 32;
inline static constexpr int NDD = 94; // XXX was 90 in Fortran
inline static constexpr int JZ = 32;
inline static constexpr float ASTART = 0.5f;
inline static constexpr float BASESUB = 42.0f;
// Derived parameters
inline static constexpr float AZ = (12000.0f / NSPS) * 0.64f;
inline static constexpr int NMAX = NTXDUR * JS8_RX_SAMPLE_RATE;
inline static constexpr int NFFT1 = NSPS * NFOS;
inline static constexpr int NSTEP = NSPS / NSSY;
inline static constexpr int NHSYM = NMAX / NSTEP - 3;
inline static constexpr int NDOWN = NSPS / NDOWNSPS;
inline static constexpr int NQSYMBOL = NDOWNSPS / 4;
inline static constexpr int NDFFT1 = NSPS * NDD;
inline static constexpr int NDFFT2 = NDFFT1 / NDOWN;
inline static constexpr int NP2 = NN * NDOWNSPS;
inline static constexpr float TSTEP = NSTEP / 12000.0f;
inline static constexpr int JSTRT = ASTART / TSTEP;
inline static constexpr float DF = 12000.0f / NFFT1;
};
/* I MODE DECODER */
struct ModeI
{
// Static constants
inline static constexpr int NSUBMODE = 8;
inline static constexpr auto NCOSTAS = JS8::Costas::Type::MODIFIED;
inline static constexpr int NSPS = JS8I_SYMBOL_SAMPLES;
inline static constexpr int NTXDUR = JS8I_TX_SECONDS;
inline static constexpr int NDOWNSPS = 12;
inline static constexpr int NDD = 125;
inline static constexpr int JZ = 250;
inline static constexpr float ASTART = 0.1f;
inline static constexpr float BASESUB = 36.0f;
// Derived parameters
inline static constexpr float AZ = (12000.0f / NSPS) * 0.64f;
inline static constexpr int NMAX = NTXDUR * JS8_RX_SAMPLE_RATE;
inline static constexpr int NFFT1 = NSPS * NFOS;
inline static constexpr int NSTEP = NSPS / NSSY;
inline static constexpr int NHSYM = NMAX / NSTEP - 3;
inline static constexpr int NDOWN = NSPS / NDOWNSPS;
inline static constexpr int NQSYMBOL = NDOWNSPS / 4;
inline static constexpr int NDFFT1 = NSPS * NDD;
inline static constexpr int NDFFT2 = NDFFT1 / NDOWN;
inline static constexpr int NP2 = NN * NDOWNSPS;
inline static constexpr float TSTEP = NSTEP / 12000.0f;
inline static constexpr int JSTRT = ASTART / TSTEP;
inline static constexpr float DF = 12000.0f / NFFT1;
};
// Tunable settings; degree of the polynomial used for the baseline
// curve fit, and the percentile of the span at which to sample. In
// general, a 5th degree polynomial and the 10th percentile should
// be optimal.
constexpr auto BASELINE_DEGREE = 5;
constexpr auto BASELINE_SAMPLE = 10;
// Define the closed range in Hz that we'll consider to be the window
// for baseline determination.
constexpr auto BASELINE_MIN = 500;
constexpr auto BASELINE_MAX = 2500;
// We're going to do a pairwise Estrin's evaluation of the polynomial
// coefficients, so it's critical that the degree of the polynomial is
// odd, resulting in an even number of coefficients.
static_assert(BASELINE_DEGREE & 1, "Degree must be odd");
static_assert(BASELINE_SAMPLE >= 0 &&
BASELINE_SAMPLE <= 100, "Sample must be a percentage");
// Since we know the degree of the polynomial, and thus the number of
// nodes that we're going to use, we can do all the trigonometry work
// required to calculate the Chebyshev nodes in advance, by computing
// them over the range [0, 1]; we can then scale these at runtime to
// a span of any size by simple multiplication.
//
// Downside to this with C++17 is that std::cos() is not yet constexpr,
// as it is in C++20, so we must provide our own implementation until
// then.
constexpr auto BASELINE_NODES = []()
{
// Down to the actual business of generating Chebyshev nodes
// suitable for scaling; once we move to C++20 as the minimum
// compiler, we can remove the cos() function above and instead
// call std::cos() here, as it's required to be constexpr in
// C++20 and above, and presumably it'll be of high quality.
auto nodes = std::array<double, BASELINE_DEGREE + 1>{};
constexpr auto slice = std::numbers::pi / (2.0 * nodes.size());
for (std::size_t i = 0; i < nodes.size(); ++i)
{
nodes[i] = 0.5 * (1.0 - cos(slice * (2.0 * i + 1)));
}
return nodes;
}();
}
/******************************************************************************/
// Local Types
/******************************************************************************/
namespace
{
// Accumulation of rounding errors in IEEE 754 values can be a problem
// when summing large numbers of small values; a Kahan summation class
// by which to compensate for them.
//
// Fortran, or at least, gfortran, will use this technique under the
// covers in various scenarios. While it'd be reasonable to expect it
// to be used in sum(), that's typically not the case.
//
// However, for example, it'll use it here for the value that goes into
// win(i), and naive summation in C++ will as a result not produce the
// same values without using compensation.
//
// subroutine nuttal_window(win,n)
// real win(n)
// pi=4.0*atan(1.0)
// a0=0.3635819
// a1=-0.4891775;
// a2=0.1365995;
// a3=-0.0106411;
// do i=1,n
// win(i)=a0+a1*cos(2*pi*(i-1)/(n))+ &
// a2*cos(4*pi*(i-1)/(n))+ &
// a3*cos(6*pi*(i-1)/(n))
// enddo
// return
// end subroutine nuttal_window
template <std::floating_point T>
class KahanSum
{
T m_sum; // Accumulated sum
T m_compensation; // Compensation for lost low-order bits
public:
KahanSum(T sum = T(0))
: m_sum(sum)
, m_compensation(T(0))
{}
KahanSum &
operator=(T const sum)
{
m_sum = sum;
m_compensation = T(0);
return *this;
}
KahanSum &
operator+=(T const value)
{
T const y = value - m_compensation; // Correct the value
T const t = m_sum + y; // Perform the sum
m_compensation = (t - m_sum) - y; // Update compensation
m_sum = t; // Update the sum
return *this;
}
operator T() const { return m_sum; }
};
// Management of dynamic FFTW plan storage.
class FFTWPlanManager
{
public:
enum class Type
{
DS,
BB,
CF,
CB,
SD,
CS,
count
};
// Disallow copying and moving
FFTWPlanManager (FFTWPlanManager const &) = delete;
FFTWPlanManager & operator=(FFTWPlanManager const &) = delete;
FFTWPlanManager (FFTWPlanManager &&) = delete;
FFTWPlanManager & operator=(FFTWPlanManager &&) = delete;
// Constructor
FFTWPlanManager()
{
m_plans.fill(nullptr);
}
// Destructor
~FFTWPlanManager()
{
std::lock_guard<std::mutex> lock(fftw_mutex);
for (auto & plan : m_plans)
{
if (plan) fftwf_destroy_plan(plan);
}
}
// Accessor
fftwf_plan const &
operator[](Type const type) const noexcept
{
return m_plans[static_cast<std::size_t>(type)];
}
// Manipulator
fftwf_plan &
operator[](Type const type) noexcept
{
return m_plans[static_cast<std::size_t>(type)];
}
// Iteration support
auto begin() noexcept { return m_plans.begin(); }
auto end() noexcept { return m_plans.end(); }
auto begin() const noexcept { return m_plans.begin(); }
auto end() const noexcept { return m_plans.end(); }
private:
// Data members
std::array<fftwf_plan, static_cast<std::size_t>(Type::count)> m_plans;
};
// Encapsulates the first-order search results provided by syncjs8().
struct Sync
{
float freq;
float step;
float sync;
// Constructor for convenience.
Sync(float const freq,
float const step,
float const sync)
: freq(freq)
, step(step)
, sync(sync)
{}
};
// Tag structs so that we can refer to multi index container indices
// by a descriptive tag instead of by the index of the index. These
// don't need to be anything but a name.
namespace Tag
{
struct Freq {};
struct Rank {};
struct Sync {};
}
// Container indexing Sync objects in useful ways, used by syncjs8().
namespace MI = boost::multi_index;
using SyncIndex = MI::multi_index_container
<
Sync,
MI::indexed_by
<
MI::ordered_non_unique<
MI::tag<Tag::Freq>,
MI::key<&Sync::freq>
>,
MI::ranked_non_unique<
MI::tag<Tag::Rank>,
MI::key<&Sync::sync>
>,
MI::ordered_non_unique<
MI::tag<Tag::Sync>,
MI::key<&Sync::sync>,
std::greater<>
>
>
>;
// Represents a decoded message, i.e., the 3-bit message type
// and the 12 bytes that result from decoding a message.
class Decode
{
public:
int type;
std::string data;
Decode(int type,
std::string data)
: type(type)
, data(std::move(data))
{}
bool operator==(Decode const &) const noexcept = default;
struct Hash
{
std::size_t
operator()(Decode const & decode) const noexcept
{
std::size_t const h1 = std::hash<int>{}(decode.type);
std::size_t const h2 = std::hash<std::string>{}(decode.data);
return h1 ^ (h2 + 0x9e3779b9 + (h1 << 6) + (h1 >> 2));
}
};
using Map = std::unordered_map<Decode, int, Hash>;
};
}
/******************************************************************************/
// Belief Propagation Decoder
/******************************************************************************/
namespace
{
constexpr int BP_MAX_ROWS = 7; // Max rows per column in Nm
constexpr int BP_MAX_CHECKS = 3; // Max checks per bit in Mn
constexpr int BP_MAX_ITERATIONS = 30; // Max iterations in BP decoder
constexpr std::array<std::array<int, BP_MAX_CHECKS>, N> Mn =
{{
{ 0, 24, 68}, { 1, 4, 72}, { 2, 31, 67}, { 3, 50, 60}, { 5, 62, 69}, { 6, 32, 78},
{ 7, 49, 85}, { 8, 36, 42}, { 9, 40, 64}, {10, 13, 63}, {11, 74, 76}, {12, 22, 80},
{14, 15, 81}, {16, 55, 65}, {17, 52, 59}, {18, 30, 51}, {19, 66, 83}, {20, 28, 71},
{21, 23, 43}, {25, 34, 75}, {26, 35, 37}, {27, 39, 41}, {29, 53, 54}, {33, 48, 86},
{38, 56, 57}, {44, 73, 82}, {45, 61, 79}, {46, 47, 84}, {58, 70, 77}, { 0, 49, 52},
{ 1, 46, 83}, { 2, 24, 78}, { 3, 5, 13}, { 4, 6, 79}, { 7, 33, 54}, { 8, 35, 68},
{ 9, 42, 82}, {10, 22, 73}, {11, 16, 43}, {12, 56, 75}, {14, 26, 55}, {15, 27, 28},
{17, 18, 58}, {19, 39, 62}, {20, 34, 51}, {21, 53, 63}, {23, 61, 77}, {25, 31, 76},
{29, 71, 84}, {30, 64, 86}, {32, 38, 50}, {36, 47, 74}, {37, 69, 70}, {40, 41, 67},
{44, 66, 85}, {45, 80, 81}, {48, 65, 72}, {57, 59, 65}, {60, 64, 84}, { 0, 13, 20},
{ 1, 12, 58}, { 2, 66, 81}, { 3, 31, 72}, { 4, 35, 53}, { 5, 42, 45}, { 6, 27, 74},
{ 7, 32, 70}, { 8, 48, 75}, { 9, 57, 63}, {10, 47, 67}, {11, 18, 44}, {14, 49, 60},
{15, 21, 25}, {16, 71, 79}, {17, 39, 54}, {19, 34, 50}, {22, 24, 33}, {23, 62, 86},
{26, 38, 73}, {28, 77, 82}, {29, 69, 76}, {30, 68, 83}, {21, 36, 85}, {37, 40, 80},
{41, 43, 56}, {46, 52, 61}, {51, 55, 78}, {59, 74, 80}, { 0, 38, 76}, { 1, 15, 40},
{ 2, 30, 53}, { 3, 35, 77}, { 4, 44, 64}, { 5, 56, 84}, { 6, 13, 48}, { 7, 20, 45},
{ 8, 14, 71}, { 9, 19, 61}, {10, 16, 70}, {11, 33, 46}, {12, 67, 85}, {17, 22, 42},
{18, 63, 72}, {23, 47, 78}, {24, 69, 82}, {25, 79, 86}, {26, 31, 39}, {27, 55, 68},
{28, 62, 65}, {29, 41, 49}, {32, 36, 81}, {34, 59, 73}, {37, 54, 83}, {43, 51, 60},
{50, 52, 71}, {57, 58, 66}, {46, 55, 75}, { 0, 18, 36}, { 1, 60, 74}, { 2, 7, 65},
{ 3, 59, 83}, { 4, 33, 38}, { 5, 25, 52}, { 6, 31, 56}, { 8, 51, 66}, { 9, 11, 14},
{10, 50, 68}, {12, 13, 64}, {15, 30, 42}, {16, 19, 35}, {17, 79, 85}, {20, 47, 58},
{21, 39, 45}, {22, 32, 61}, {23, 29, 73}, {24, 41, 63}, {26, 48, 84}, {27, 37, 72},
{28, 43, 80}, {34, 67, 69}, {40, 62, 75}, {44, 48, 70}, {49, 57, 86}, {47, 53, 82},
{12, 54, 78}, {76, 77, 81}, { 0, 1, 23}, { 2, 5, 74}, { 3, 55, 86}, { 4, 43, 52},
{ 6, 49, 82}, { 7, 9, 27}, { 8, 54, 61}, {10, 28, 66}, {11, 32, 39}, {13, 15, 19},
{14, 34, 72}, {16, 30, 38}, {17, 35, 56}, {18, 45, 75}, {20, 41, 83}, {21, 33, 58},
{22, 25, 60}, {24, 59, 64}, {26, 63, 79}, {29, 36, 65}, {31, 44, 71}, {37, 50, 85},
{40, 76, 78}, {42, 55, 67}, {46, 73, 81}, {39, 51, 77}, {53, 60, 70}, {45, 57, 68}
}};
struct CheckNode
{
int valid_neighbors;
std::array<int, BP_MAX_ROWS> neighbors;
};
constexpr std::array<CheckNode, M> Nm =
{{
{6, { 0, 29, 59, 88, 117, 146, 0}}, {6, { 1, 30, 60, 89, 118, 146, 0}}, {6, { 2, 31, 61, 90, 119, 147, 0}},
{6, { 3, 32, 62, 91, 120, 148, 0}}, {6, { 1, 33, 63, 92, 121, 149, 0}}, {6, { 4, 32, 64, 93, 122, 147, 0}},
{6, { 5, 33, 65, 94, 123, 150, 0}}, {6, { 6, 34, 66, 95, 119, 151, 0}}, {6, { 7, 35, 67, 96, 124, 152, 0}},
{6, { 8, 36, 68, 97, 125, 151, 0}}, {6, { 9, 37, 69, 98, 126, 153, 0}}, {6, {10, 38, 70, 99, 125, 154, 0}},
{6, {11, 39, 60, 100, 127, 144, 0}}, {6, { 9, 32, 59, 94, 127, 155, 0}}, {6, {12, 40, 71, 96, 125, 156, 0}},
{6, {12, 41, 72, 89, 128, 155, 0}}, {6, {13, 38, 73, 98, 129, 157, 0}}, {6, {14, 42, 74, 101, 130, 158, 0}},
{6, {15, 42, 70, 102, 117, 159, 0}}, {6, {16, 43, 75, 97, 129, 155, 0}}, {6, {17, 44, 59, 95, 131, 160, 0}},
{6, {18, 45, 72, 82, 132, 161, 0}}, {6, {11, 37, 76, 101, 133, 162, 0}}, {6, {18, 46, 77, 103, 134, 146, 0}},
{6, { 0, 31, 76, 104, 135, 163, 0}}, {6, {19, 47, 72, 105, 122, 162, 0}}, {6, {20, 40, 78, 106, 136, 164, 0}},
{6, {21, 41, 65, 107, 137, 151, 0}}, {6, {17, 41, 79, 108, 138, 153, 0}}, {6, {22, 48, 80, 109, 134, 165, 0}},
{6, {15, 49, 81, 90, 128, 157, 0}}, {6, {2, 47, 62, 106, 123, 166, 0}}, {6, { 5, 50, 66, 110, 133, 154, 0}},
{6, {23, 34, 76, 99, 121, 161, 0}}, {6, {19, 44, 75, 111, 139, 156, 0}}, {6, {20, 35, 63, 91, 129, 158, 0}},
{6, { 7, 51, 82, 110, 117, 165, 0}}, {6, {20, 52, 83, 112, 137, 167, 0}}, {6, {24, 50, 78, 88, 121, 157, 0}},
{7, {21, 43, 74, 106, 132, 154, 171}}, {6, { 8, 53, 83, 89, 140, 168, 0}}, {6, {21, 53, 84, 109, 135, 160, 0}},
{6, { 7, 36, 64, 101, 128, 169, 0}}, {6, {18, 38, 84, 113, 138, 149, 0}}, {6, {25, 54, 70, 92, 141, 166, 0}},
{7, {26, 55, 64, 95, 132, 159, 173}}, {6, {27, 30, 85, 99, 116, 170, 0}}, {6, {27, 51, 69, 103, 131, 143, 0}},
{6, {23, 56, 67, 94, 136, 141, 0}}, {6, {6, 29, 71, 109, 142, 150, 0}}, {6, { 3, 50, 75, 114, 126, 167, 0}},
{6, {15, 44, 86, 113, 124, 171, 0}}, {6, {14, 29, 85, 114, 122, 149, 0}}, {6, {22, 45, 63, 90, 143, 172, 0}},
{6, {22, 34, 74, 112, 144, 152, 0}}, {7, {13, 40, 86, 107, 116, 148, 169}}, {6, {24, 39, 84, 93, 123, 158, 0}},
{6, {24, 57, 68, 115, 142, 173, 0}}, {6, {28, 42, 60, 115, 131, 161, 0}}, {6, {14, 57, 87, 111, 120, 163, 0}},
{7, { 3, 58, 71, 113, 118, 162, 172}}, {6, {26, 46, 85, 97, 133, 152, 0}}, {5, { 4, 43, 77, 108, 140, 0, 0}},
{6, { 9, 45, 68, 102, 135, 164, 0}}, {6, { 8, 49, 58, 92, 127, 163, 0}}, {6, {13, 56, 57, 108, 119, 165, 0}},
{6, {16, 54, 61, 115, 124, 153, 0}}, {6, { 2, 53, 69, 100, 139, 169, 0}}, {6, { 0, 35, 81, 107, 126, 173, 0}},
{5, { 4, 52, 80, 104, 139, 0, 0}}, {6, {28, 52, 66, 98, 141, 172, 0}}, {6, {17, 48, 73, 96, 114, 166, 0}},
{6, { 1, 56, 62, 102, 137, 156, 0}}, {6, {25, 37, 78, 111, 134, 170, 0}}, {6, {10, 51, 65, 87, 118, 147, 0}},
{6, {19, 39, 67, 116, 140, 159, 0}}, {6, {10, 47, 80, 88, 145, 168, 0}}, {6, {28, 46, 79, 91, 145, 171, 0}},
{6, { 5, 31, 86, 103, 144, 168, 0}}, {6, {26, 33, 73, 105, 130, 164, 0}}, {5, {11, 55, 83, 87, 138, 0, 0}},
{6, {12, 55, 61, 110, 145, 170, 0}}, {6, {25, 36, 79, 104, 143, 150, 0}}, {6, {16, 30, 81, 112, 120, 160, 0}},
{5, {27, 48, 58, 93, 136, 0, 0}}, {6, { 6, 54, 82, 100, 130, 167, 0}}, {6, {23, 49, 77, 105, 142, 148, 0}}
}};
// Belief Propagation Decoder
int
bpdecode174(std::array<float, N> const & llr,
std::array<int8_t, K> & decoded,
std::array<int8_t, N> & cw)
{
// Initialize messages and variables
std::array<std::array<float, BP_MAX_CHECKS>, N> tov = {}; // Messages to variable nodes
std::array<std::array<float, BP_MAX_ROWS>, M> toc = {}; // Messages to check nodes
std::array<std::array<float, BP_MAX_ROWS> , M> tanhtoc = {}; // Tanh of messages
std::array<float, N> zn = {}; // Bit log likelihood ratios
std::array<int, M> synd = {}; // Syndrome for checks
int ncnt = 0;
int nclast = 0;
// Initialize toc (messages from bits to checks)
for (int i = 0; i < M; ++i) {
for (int j = 0; j < Nm[i].valid_neighbors; ++j) {
toc[i][j] = llr[Nm[i].neighbors[j]];
}
}
// Iterative decoding
for (int iter = 0; iter <= BP_MAX_ITERATIONS; ++iter) {
// Update bit log likelihood ratios
for (int i = 0; i < N; ++i) {
zn[i] = llr[i] + std::accumulate(tov[i].begin(), tov[i].begin() + BP_MAX_CHECKS, 0.0f);
}
// Check if we have a valid codeword
for (int i = 0; i < N; ++i) cw[i] = zn[i] > 0 ? 1 : 0;
int ncheck = 0;
for (int i = 0; i < M; ++i) {
synd[i] = 0;
for (int j = 0; j < Nm[i].valid_neighbors; ++j) {
synd[i] += cw[Nm[i].neighbors[j]];
}
if (synd[i] % 2 != 0) ++ncheck;
}
if (ncheck == 0)
{
// Extract decoded bits (last N-M bits of codeword)
std::copy(cw.begin() + M, cw.end(), decoded.begin());
// Count errors
int nerr = 0;
for (int i = 0; i < N; ++i) {
if ((2 * cw[i] - 1) * llr[i] < 0.0f) {
++nerr;
}
}
return nerr;
}
// Early stopping criterion
if (iter > 0) {
int nd = ncheck - nclast;
ncnt = (nd < 0) ? 0 : ncnt + 1;
if (ncnt >= 5 && iter >= 10 && ncheck > 15) {
return -1;
}
}
nclast = ncheck;
// Send messages from bits to check nodes
for (int i = 0; i < M; ++i) {
for (int j = 0; j < Nm[i].valid_neighbors; ++j) {
int ibj = Nm[i].neighbors[j];
toc[i][j] = zn[ibj];
for (int k = 0; k < BP_MAX_CHECKS; ++k) {
if (Mn[ibj][k] == i) {
toc[i][j] -= tov[ibj][k];
}
}
}
}
// Send messages from check nodes to variable nodes
for (int i = 0; i < M; ++i) {
for (int j = 0; j < 7; ++j) { // Fixed range [0, 7) to match Fortran's 1:7, could be nrw[j], or 7 logically
tanhtoc[i][j] = std::tanh(-toc[i][j] / 2.0f);
}
}
for (int i = 0; i < N; ++i) {
for (int j = 0; j < BP_MAX_CHECKS; ++j) {
int ichk = Mn[i][j];
if (ichk >= 0) {
float Tmn = 1.0f;
for (int k = 0; k < Nm[ichk].valid_neighbors; ++k) {
if (Nm[ichk].neighbors[k] != i) {
Tmn *= tanhtoc[ichk][k];
}
}
tov[i][j] = 2.0f * std::atanh(-Tmn);
}
}
}
}
return -1; // Decoding failed
}
}
/******************************************************************************/
// Local Routines
/******************************************************************************/
namespace
{
constexpr std::string_view alphabet = "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz-+";
static_assert(alphabet.size() == 64);
// Function that either translates valid JS8 message characters to their
// corresponding 6-bit word value, or throws. This will end up doing a
// direct index operation into a 256-byte table, the creation of which
// must be constexpr under C++17.
constexpr auto alphabetWord = []()
{
constexpr std::uint8_t invalid = 0xff;
constexpr auto words = []()
{
std::array<std::uint8_t, 256> words{};
for (auto & word : words) word = invalid;
for (std::size_t i = 0; i < alphabet.size(); ++i)
{
words[static_cast<std::uint8_t>(alphabet[i])] = static_cast<std::uint8_t>(i);
}
return words;
}();
return [words](char const value)
{
if (auto const word = words[value];
word != invalid)
{
return word;
}
throw std::runtime_error("Invalid character in message");
};
}();
// Sanity check key bounds of the 6-bit encoding table.
static_assert(alphabetWord('0') == 0);
static_assert(alphabetWord('A') == 10);
static_assert(alphabetWord('a') == 36);
static_assert(alphabetWord('-') == 62);
static_assert(alphabetWord('+') == 63);
template <typename T>
std::uint16_t
CRC12(T const & range)
{
return boost::augmented_crc<12, 0xc06>(range.data(),
range.size()) ^ 42;
}
bool
checkCRC12(std::array<std::int8_t, KK> const & decoded)
{
std::array<uint8_t, 11> bits = {};
for (std::size_t i = 0; i < decoded.size(); ++i)
{
if (decoded[i]) bits[i / 8] |= (1 << (7 - (i % 8)));
}
// Extract the received CRC-12.
uint16_t crc = (static_cast<uint16_t>(bits[9] & 0x1F) << 7) |
(static_cast<uint16_t>(bits[10]) >> 1);
// Clear bits that correspond to the CRC in the last bytes.
bits[9] &= 0xE0;
bits[10] = 0x00;
// Compute CRC and indicate if we have a match.
return crc == CRC12(bits);
}
std::string
extractmessage174(std::array<int8_t, KK> const & decoded)
{
std::string message;
// Ensure received CRC matches computed CRC.
if (checkCRC12(decoded))
{
message.reserve(12);
// Decode the message from the 72 data bits
std::array<uint8_t, 12> words;
for (std::size_t i = 0; i < 12; ++i)
{
words[i] = (decoded[i * 6 + 0] << 5) |
(decoded[i * 6 + 1] << 4) |
(decoded[i * 6 + 2] << 3) |
(decoded[i * 6 + 3] << 2) |
(decoded[i * 6 + 4] << 1) |
(decoded[i * 6 + 5] << 0);
}
// Map 6-bit words to the alphabet
for (auto const word : words) message += alphabet[word];
}
return message;
}
// Parity matrix for JS8 message generation.
//
// This should be 952 bytes in size; to store an 87x87 matrix of bits,
// you need 7569 bits, which requires 119 64-bit values, or 952 bytes.
//
// Background here is that this is a low-density parity check code (LDPC),
// generated using the PEG algorithm. In short, true values in a row i of
// the matrix define which of the 87 message bits must be summed, modulo
// 2, to produce the ith parity check bit. Decent references on this are:
//
// 1. https://wsjt.sourceforge.io/FT4_FT8_QEX.pdf
// 2. https://inference.org.uk/mackay/PEG_ECC.html
// 3. https://github.com/Lcrypto/classic-PEG-
//
// The data used was harvested from the original 'ldpc_174_87_params.f90',
// but you'll note that the rows have been reordered here, because this
// isn't Fortran; C++ is row-major, not column-major.
constexpr auto parity = []()
{
constexpr std::size_t Rows = 87;
constexpr std::size_t Cols = 87;
using ElementType = std::uint64_t;
constexpr std::size_t ElementSize = std::numeric_limits<ElementType>::digits;
constexpr auto matrix = []()
{
constexpr std::array<std::string_view, Rows> Data =
{
"23bba830e23b6b6f50982e", "1f8e55da218c5df3309052", "ca7b3217cd92bd59a5ae20",
"56f78313537d0f4382964e", "6be396b5e2e819e373340c", "293548a138858328af4210",
"cb6c6afcdc28bb3f7c6e86", "3f2a86f5c5bd225c961150", "849dd2d63673481860f62c",
"56cdaec6e7ae14b43feeee", "04ef5cfa3766ba778f45a4", "c525ae4bd4f627320a3974",
"41fd9520b2e4abeb2f989c", "7fb36c24085a34d8c1dbc4", "40fc3e44bb7d2bb2756e44",
"d38ab0a1d2e52a8ec3bc76", "3d0f929ef3949bd84d4734", "45d3814f504064f80549ae",
"f14dbf263825d0bd04b05e", "db714f8f64e8ac7af1a76e", "8d0274de71e7c1a8055eb0",
"51f81573dd4049b082de14", "d8f937f31822e57c562370", "b6537f417e61d1a7085336",
"ecbd7c73b9cd34c3720c8a", "3d188ea477f6fa41317a4e", "1ac4672b549cd6dba79bcc",
"a377253773ea678367c3f6", "0dbd816fba1543f721dc72", "ca4186dd44c3121565cf5c",
"29c29dba9c545e267762fe", "1616d78018d0b4745ca0f2", "fe37802941d66dde02b99c",
"a9fa8e50bcb032c85e3304", "83f640f1a48a8ebc0443ea", "3776af54ccfbae916afde6",
"a8fc906976c35669e79ce0", "f08a91fb2e1f78290619a8", "cc9da55fe046d0cb3a770c",
"d36d662a69ae24b74dcbd8", "40907b01280f03c0323946", "d037db825175d851f3af00",
"1bf1490607c54032660ede", "0af7723161ec223080be86", "eca9afa0f6b01d92305edc",
"7a8dec79a51e8ac5388022", "9059dfa2bb20ef7ef73ad4", "6abb212d9739dfc02580f2",
"f6ad4824b87c80ebfce466", "d747bfc5fd65ef70fbd9bc", "612f63acc025b6ab476f7c",
"05209a0abb530b9e7e34b0", "45b7ab6242b77474d9f11a", "6c280d2a0523d9c4bc5946",
"f1627701a2d692fd9449e6", "8d9071b7e7a6a2eed6965e", "bf4f56e073271f6ab4bf80",
"c0fc3ec4fb7d2bb2756644", "57da6d13cb96a7689b2790", "a9fa2eefa6f8796a355772",
"164cc861bdd803c547f2ac", "cc6de59755420925f90ed2", "a0c0033a52ab6299802fd2",
"b274db8abd3c6f396ea356", "97d4169cb33e7435718d90", "81cfc6f18c35b1e1f17114",
"481a2a0df8a23583f82d6c", "081c29a10d468ccdbcecb6", "2c4142bf42b01e71076acc",
"a6573f3dc8b16c9d19f746", "c87af9a5d5206abca532a8", "012dee2198eba82b19a1da",
"b1ca4ea2e3d173bad4379c", "b33ec97be83ce413f9acc8", "5b0f7742bca86b8012609a",
"37d8e0af9258b9e8c5f9b2", "35ad3fb0faeb5f1b0c30dc", "6114e08483043fd3f38a8a",
"cd921fdf59e882683763f6", "95e45ecd0135aca9d6e6ae", "2e547dd7a05f6597aac516",
"14cd0f642fc0c5fe3a65ca", "3a0a1dfd7eee29c2e827e0", "c8b5dffc335095dcdcaf2a",
"3dd01a59d86310743ec752", "8abdb889efbe39a510a118", "3f231f212055371cf3e2a2"
};
constexpr std::size_t Total = (Rows * Cols + ElementSize - 1);
constexpr std::size_t Count = Total / ElementSize;
constexpr std::array<std::uint8_t, 4> Masks = {0x8, 0x4, 0x2, 0x1};
std::array<ElementType, Count> data{};
for (std::size_t row = 0; row < Rows; ++row)
{
std::size_t col = 0;
for (auto const c : Data[row])
{
std::uint8_t const value = (c >= '0' && c <= '9') ? c - '0' :
(c >= 'a' && c <= 'f') ? c - 'a' + 10 :
(c >= 'A' && c <= 'F') ? c - 'A' + 10 : throw "Invalid hex";
for (auto const mask : Masks)
{
if (col >= Cols) break;
if (value & mask)
{
auto const index = row * Cols + col;
data[index / ElementSize] |= (ElementType(1) << (index % ElementSize));
}
++col;
}
}
}
return data;
}();
return [matrix](std::size_t const row,
std::size_t const col)
{
auto const index = row * Cols + col;
return (matrix[index / ElementSize] >>
(index % ElementSize)) & 1;
};
}();
}
/******************************************************************************/
// DecodeMode Template Class
/******************************************************************************/
// Mode-parameterized decode class.
namespace
{
template <typename Mode>
class DecodeMode
{
// Data members
std::array<float, Mode::NFFT1> nuttal;
std::array<std::array<std::array<std::complex<float>, Mode::NDOWNSPS>, 7>, 3> csyncs;
alignas(64) std::array<std::complex<float>, Mode::NDOWNSPS> csymb;
alignas(64) std::array<std::complex<float>, Mode::NMAX> filter;
alignas(64) std::array<std::complex<float>, Mode::NMAX> cfilt;
alignas(64) std::array<std::complex<float>, Mode::NDFFT1 / 2 + 1> ds_cx;
alignas(64) std::array<std::complex<float>, Mode::NFFT1 / 2 + 1> sd;
alignas(64) std::array<std::complex<float>, NP> cd0;
std::array<float, Mode::NMAX> dd;
std::array<std::array<float, Mode::NHSYM>, Mode::NSPS> s;
std::array<float, Mode::NSPS> savg;
FFTWPlanManager plans;
SyncIndex sync;
using Plan = FFTWPlanManager::Type;
static constexpr auto Costas = JS8::Costas::array(Mode::NCOSTAS);
// Fore and aft tapers to reduce spectral leakage during the
// downsampling process. We can compute these at compile time.
static constexpr auto Taper = []
{
std::array<std::array<float, Mode::NDD + 1>, 2> taper{};
for (size_t i = 0; i <= Mode::NDD; ++i)
{
float const value = 0.5f * (1.0f + cos(i * std::numbers::pi_v<float> / Mode::NDD));
taper[1][ i] = value; // TailTaper (original taper)
taper[0][Mode::NDD - i] = value; // HeadTaper (reversed taper)
}
return taper;
}();
// Baseline computation support.
using Points = Eigen::Matrix<double, BASELINE_NODES.size(), 2>;
using Vandermonde = Eigen::Matrix<double, BASELINE_NODES.size(),
BASELINE_NODES.size()>;
using Coefficients = Eigen::Vector<double, BASELINE_NODES.size()>;
Points p;
Vandermonde V;
Coefficients c;
// Polynomial evaluation using Estrin's method, loop is unrolled at
// compile time. A compiler should emit SIMD instructions from what
// it sees here when the optimizer is involved, but even without it,
// we'll likely see fused multiply-add instructions.
inline auto
evaluate(float const x) const
{
return [this]<Eigen::Index... I>(float const x,
std::integer_sequence<Eigen::Index, I...>)
{
auto baseline = 0.0;
auto exponent = 1.0;
((baseline += (c[I * 2] + c[I * 2 + 1] * x) * exponent, exponent *= x * x), ...);
return static_cast<float>(baseline);
}(x, std::make_integer_sequence<Eigen::Index,
Coefficients::SizeAtCompileTime / 2>{});
}
std::optional<Decode>
js8dec(bool const syncStats,
bool const lsubtract,
float & f1,
float & xdt,
int & nharderrors,
float & xsnr,
JS8::Event::Emitter emitEvent)
{
constexpr float FR = 12000.0f / Mode::NFFT1; // Frequency resolution
constexpr float FS2 = 12000.0f / Mode::NDOWN;
constexpr float DT2 = 1.0f / FS2;
auto const index = static_cast<int>(std::round(f1 / FR)); // Closest index
float const scaled_value = 0.1f * (savg[index] - Mode::BASESUB); // Adjust and scale
float const xbase = std::pow(10.0f, scaled_value); // Convert to linear scale
float delfbest = 0.0f;
int ibest = 0;
// Downsample the signal and prepare for processing.
js8_downsample(f1);
// Initial guess for the start of the signal.
int i0 = static_cast<int>(std::round((xdt + Mode::ASTART) * FS2));
float smax = 0.0f;
// Search for the best synchronization offset.
for (int idt = i0 - Mode::NQSYMBOL;
idt <= i0 + Mode::NQSYMBOL;
++idt)
{
float const sync = syncjs8d(idt, 0.0f);
if (sync > smax) {
smax = sync;
ibest = idt;
}
}
// Improved estimate for DT.
float const xdt2 = ibest * DT2;
// Fine frequency synchronization
i0 = static_cast<int>(std::round(xdt2 * FS2));
smax = 0.0f;
for (int ifr = -NFSRCH;
ifr <= NFSRCH;
++ifr)
{
float const delf = ifr * 0.5f;
float const sync = syncjs8d(i0, delf);
if (sync > smax) {
smax = sync;
delfbest = delf;
}
}
// Frequency tweaking.
float const dphi = -delfbest * ((2.0f * std::numbers::pi_v<float>) / FS2); // Phase increment
std::complex<float> const wstep = std::polar(1.0f, dphi); // Step for phase rotation
std::complex<float> w = {1.0f, 0.0f}; // Cumlative phase
for (int i = 0; i < NP2; ++i)
{
w *= wstep; // Update cumulative phase
cd0[i] *= w; // Apply phase shift
}
// Adjust the frequency and time offset.
xdt = xdt2;
f1 += delfbest;
float const sync = syncjs8d(i0, 0.0f);
std::array<std::array<float, NN>, NROWS> s2;
for (int k = 0; k < NN; ++k)
{
// Calculate the starting index for the current symbol.
int const i1 = ibest + k * Mode::NDOWNSPS;
csymb.fill(ZERO);
if (i1 >= 0 && i1 + Mode::NDOWNSPS <= NP2)
{
std::copy(cd0.begin() + i1,
cd0.begin() + i1 + Mode::NDOWNSPS,
csymb.begin());
}
fftwf_execute(plans[Plan::CS]);
// Normalize and take the magnitude of the first 8 points.
for (int i = 0; i < NROWS; ++i)
{
s2[i][k] = std::abs(csymb[i]) / 1000.0f;
}
}
// Sync quality check using Costas tone patterns.
int nsync = 0;
for (std::size_t costas = 0; costas < Costas.size(); ++costas)
{
auto const offset = costas * 36;
for (std::size_t column = 0; column < 7; ++column)
{
// Find the row containing the maximum value in the
// current column.
auto const max_row = std::distance(
s2.begin(),
std::max_element(s2.begin(),
s2.end(),
[index = offset + column]
(auto const & rowA,
auto const & rowB)
{
return rowA[index] < rowB[index];
}));
// Check if the max row matches the Costas pattern.
if (Costas[costas][column] == max_row) ++nsync;
}
}
// If the sync quality isn't at least 7, this one's a loser.
if (nsync <= 6) return std::nullopt;
if (syncStats) emitEvent(JS8::Event::SyncState{JS8::Event::SyncState::Type::CANDIDATE,
Mode::NSUBMODE,
f1,
xdt,
{.candidate = nsync}});
std::array<std::array<float, ND>, NROWS> s1;
// Fill s1 from s2, excluding the Costas arrays.
for (int row = 0; row < NROWS; ++row)
{
std::copy(s2[row].begin() + 7, s2[row].begin() + 36, s1[row].begin());
std::copy(s2[row].begin() + 43, s2[row].begin() + 72, s1[row].begin() + 29);
}
// Temporary variables for metrics
std::array<float, 3 * ND> llr0 = {};
std::array<float, 3 * ND> llr1 = {};
// Compute metrics for each row in `s1`
for (int j = 0; j < ND; ++j)
{
int const i1 = 3 * j; // First column (matches Fortran's i1)
int const i2 = 3 * j + 1; // Second column (matches Fortran's i2)
int const i4 = 3 * j + 2; // Third column (matches Fortran's i4)
std::array<float, NROWS> ps;
for (int i = 0; i < NROWS; ++i) ps[i] = s1[i][j];
// Assign to `bmeta` in column order, with correct values
llr0[i1] = std::max({ps[4], ps[5], ps[6], ps[7]}) - std::max({ps[0], ps[1], ps[2], ps[3]}); // r4
llr0[i2] = std::max({ps[2], ps[3], ps[6], ps[7]}) - std::max({ps[0], ps[1], ps[4], ps[5]}); // r2
llr0[i4] = std::max({ps[1], ps[3], ps[5], ps[7]}) - std::max({ps[0], ps[2], ps[4], ps[6]}); // r1
for (auto & x : ps) x = std::log(x + 1e-32f);
// Assign to `bmetb` in column order, with correct values
llr1[i1] = std::max({ps[4], ps[5], ps[6], ps[7]}) - std::max({ps[0], ps[1], ps[2], ps[3]}); // r4
llr1[i2] = std::max({ps[2], ps[3], ps[6], ps[7]}) - std::max({ps[0], ps[1], ps[4], ps[5]}); // r2
llr1[i4] = std::max({ps[1], ps[3], ps[5], ps[7]}) - std::max({ps[0], ps[2], ps[4], ps[6]}); // r1
}
auto const normalizeLLR = [](auto & llr)
{
float sum = 0.0f;
float sum_of_squares = 0.0f;
for (auto const value : llr)
{
sum += value;
sum_of_squares += value * value;
}
float const llrav = sum / llr.size();
float const llr2av = sum_of_squares / llr.size();
float const variance = llr2av - llrav * llrav;
float const llrsig = std::sqrt(variance > 0.0f ? variance : llr2av);
for (float & val : llr) val = (val / llrsig) * 2.83f;
};
// Normalize and process metrics
normalizeLLR(llr0);
normalizeLLR(llr1);
std::array<int8_t, K> decoded;
std::array<int8_t, N> cw;
// Loop over decoding passes
for (int ipass = 1; ipass <= 4; ++ipass)
{
// LLR 0 used on passes 1, 3, and 4; LLR 1 used on pass 2.
auto const & llr = ipass == 2 ? llr1 : llr0;
// Zero the first 24 bytes of LLR 0 on the third pass;
// the first 48 bytes of LLR 0 on the fourth pass;
if (ipass == 3) std::fill(llr0.begin(), llr0.begin() + 24, 0.0f);
else if (ipass == 4) std::fill(llr0.begin() + 24, llr0.begin() + 48, 0.0f);
// Decode using belief propagation.
nharderrors = bpdecode174(llr, decoded, cw);
xsnr = -99.0f;
// Check for all-zero codeword
if (std::all_of(cw.begin(), cw.end(), [](int x) { return x == 0; }))
{
continue;
}
if (nharderrors >= 0 && nharderrors < 60 &&
!(sync < 2.0f && nharderrors > 35) &&
!(ipass > 2 && nharderrors > 39) &&
!(ipass == 4 && nharderrors > 30))
{
if (checkCRC12(decoded))
{
if (syncStats) emitEvent(JS8::Event::SyncState{JS8::Event::SyncState::Type::DECODED,
Mode::NSUBMODE,
f1,
xdt2,
{.decoded = sync}});
auto message = extractmessage174(decoded);
int const i3bit = (decoded[72] << 2) |
(decoded[73] << 1) |
decoded[74];
std::array<int, NN> itone;
JS8::encode(i3bit, Costas, message.data(), itone.data());
// Subtract signal if needed.
if (lsubtract) subtractjs8(genjs8refsig(itone, f1), xdt2);
// Compute the signal power.
float xsig = 0.0f;
for (std::size_t i = 0; i < itone.size(); ++i)
{
xsig += std::pow(s2[itone[i]][i], 2);
}
// Compute SNR, clamping results lower than -28 to -28.
// Note that std::log10(1.259e-10) is about -9.9; we're
// avoiding undefined behavior in the log10 computation.
xsnr = std::max(
10.0f * std::log10(std::max(
xsig / xbase - 1.0f,
1.259e-10f)) - 32.0f,
-60.0f); // XXX was -28.0f in Fortran
return std::make_optional<Decode>(i3bit, message);
}
}
else
{
nharderrors = -1;
}
}
return std::nullopt;
}
// Compute noise baseline. We differ quite a bit from the Fortran
// implementation here.
//
// The Fortran version took `savg` as input; power scaled data from
// `syncjs8`, and produced `sbase`, the noise baseline. To accomplish
// that, it used up to 1000 lower envelope points for the polynomial
// determination, which caused some oddities when the matrix was
// ill-conditioned; we're just looking for a low-order polynomial
// here, and a massively tall matrix isn't in general going to be
// helpful there. Additionally, the methodology seemed to be very
// susceptible to Runge's phenomenon.
//
// This approach instead uses a number of Chebyshev nodes proportional
// to the polynomial degree, and we evaluate 2 kHz of `savg`, centered
// around 1.5 kHz, to determine the polynomial, figuring that that's
// going to be an optimal place to measure the 10% noise floor. We then
// map the [ia, ib] range to the domain of the polynomial to compute
// the baseline. Since `savg` would otherwise no longer be referenced
// beyond this function, we dispense with `sbase` and instead overwrite
// `savg` with the baseline.
void
baselinejs8(int const ia,
int const ib)
{
// Data referenced in savg is defined by the closed range [bmin, bmax].
// From this we can derive the size of the closed range and the number
// of points in each of the arms on either side of a node. All of these
// values can be computed at compile time.
using boost::math::ccmath::round;
constexpr auto bmin = static_cast<std::size_t>(round(BASELINE_MIN / Mode::DF));
constexpr auto bmax = static_cast<std::size_t>(round(BASELINE_MAX / Mode::DF));
constexpr auto size = bmax - bmin + 1;
constexpr auto arm = size / (2 * BASELINE_NODES.size());
// Loop invariants; beginning of the data range, sentinel one past the
// end of the range.
auto const data = savg.begin() + bmin;
auto const end = data + size;
// Convert savg range of interest from power scale to dB scale.
std::transform(data,
end,
data,
[](float const value)
{
return 10.0f * std::log10(value);
});
// Collect lower envelope points; use Chebyshev node interpolants
// to reduce Runge's phenomenon oscillations.
for (std::size_t i = 0; i < BASELINE_NODES.size(); ++i)
{
auto const node = size * BASELINE_NODES[i];
auto const base = data + static_cast<int>(std::round(node));
auto span = std::vector<float>(std::clamp(base - arm, data, end),
std::clamp(base + arm, data, end));
auto const n = span.size() * BASELINE_SAMPLE / 100;
std::nth_element(span.begin(), span.begin() + n, span.end());
p.row(i) << node, span[n];
}
// Extract x and y values from points and prepare the Vandermonde
// matrix, initializing the first column with 1 (x^0); remaining
// columns are filled with the Schur product.
Eigen::VectorXd x = p.col(0);
Eigen::VectorXd y = p.col(1);
V.col(0).setOnes();
for (Eigen::Index i = 1; i < V.cols(); ++i)
{
V.col(i) = V.col(i - 1).cwiseProduct(x);
}
// Solve the least squares problem for polynomial coefficients.
c = V.colPivHouseholderQr().solve(y);
// To map an index i in the range [ia, ib] to the polynomial's
// input domain [0, size - 1]:
//
// i - ia
// x = ------- * (size - 1)
// ib - ia
auto const mapIndex = [ia, ib, last = size - 1](int const i)
{
return (i - ia) * last / float(ib - ia);
};
// Replace savg with a computed baseline in the range [ia, ib].
// This might be interpolation, which should be quite accurate,
// or extrapolation, likely somewhat less so the further we get
// from the polynomial fitting domain, but hopefully still good
// enough for our purposes here.
savg.fill(0.0f);
for (int i = ia; i <= ib; ++i)
{
savg[i] = evaluate(mapIndex(i)) + 0.65f;
}
}
// Extracted from the downsampling process; this step is part of the
// frequency-domain filtering process for downsampling the JS8 signal.
// After the FFT, the resulting frequency-domain data (ds_cx) can be
// manipulated (e.g., band-pass filtered or shifted). Subsequent inverse
// FFT operations convert the filtered data back to the time domain at
// a lower sample rate, achieving the desired downsampling.
void
computeBasebandFFT()
{
// ds_dx is an array of complex<float>; we're going to do an in-place
// FFT, so we'll interpret the first half of the array as if they were
// floats, which they are.
float * fftw_real = reinterpret_cast<float *>(ds_cx.data());
// Copy in data and zero-pad any remainder; not all modes will have
// a remainder.
std::copy(dd.begin(), dd.end(), fftw_real);
std::fill(fftw_real + dd.size(), fftw_real + Mode::NDFFT1, 0.0f);
fftwf_execute(plans[Plan::BB]);
}
// This function extracts a narrow frequency band around the target frequency f0,
// applies tapering to reduce spectral artifacts, aligns the signal to the center
// frequency, performs an inverse FFT to convert the data back into the time domain,
// and normalizes the result for further processing in the JS8 decoding pipeline.
void
js8_downsample(float const f0)
{
// Frequency band extraction; identifies a narrow frequency band around the
// target frequency (f0) based on a predefined range (8.5 baud above and 1.5
// baud below). The indices of this range in the frequency-domain representation
// (ds_cx) are calculated (ib and it), and the relevant frequency-domain samples
// are extracted into cd0.
constexpr float DF = 12000.0f / Mode::NDFFT1;
constexpr float BAUD = 12000.0f / Mode::NSPS;
float const ft = f0 + 8.5f * BAUD;
float const fb = f0 - 1.5f * BAUD;
int const i0 = static_cast<int>(std::round(f0 / DF));
int const it = std::min( static_cast<int>(std::round(ft / DF)), Mode::NDFFT1 / 2);
int const ib = std::max(0, static_cast<int>(std::round(fb / DF)));
std::size_t const NDD_SIZE = Mode::NDD + 1;
std::size_t const RANGE_SIZE = it - ib + 1;
std::fill_n(cd0.begin(), Mode::NDFFT2, ZERO);
std::copy(ds_cx.begin() + ib,
ds_cx.begin() + ib + RANGE_SIZE,
cd0.begin());
// Tapering is applied to smooth the edges of the frequency band, reducing
// spectral leakage during the inverse FFT. Reversed taper at the beginning,
// normal taper at the end.
auto const head = cd0.begin();
auto const tail = cd0.begin() + RANGE_SIZE;
std::transform(head, head + NDD_SIZE, Taper[0].begin(), head, std::multiplies<>());
std::transform(tail - NDD_SIZE, tail, Taper[1].begin(), tail - NDD_SIZE, std::multiplies<>());
// The extracted frequency band is aligned to the center of the frequency domain
// representation (i0 - ib) via a cyclic shift using std::rotate. This centers
// the desired signal.
std::rotate(cd0.begin(), cd0.begin() + (i0 - ib), cd0.begin() + Mode::NDFFT2);
// An inverse FFT is performed on the frequency-domain data (cd0) to transform it
// back into the time domain, effectively yielding a downsampled, time-domain signal
// focused on the extracted narrow frequency band.
fftwf_execute(plans[Plan::DS]);
// The resulting time-domain samples are normalized by a factor derived from the
// input and output FFT sizes (Mode::NDFFT1 and Mode::NDFFT2), ensuring consistency
// in the signal’s amplitude.
float const factor = 1.0f / std::sqrt(static_cast<float>(Mode::NDFFT1) * Mode::NDFFT2);
std::transform(cd0.begin(),
cd0.end(),
cd0.begin(),
[factor](auto & value) { return value * factor; });
}
// Evaluate the synchronization power of signal segments, ranks potential candidates, and
// extracts the most promising ones for further decoding.
//
// Detailed Steps:
//
// 1. Compute Symbol Spectra:
//
// - The signal is processed in overlapping segments, with each segment multiplied by
// a Nuttall window to reduce spectral leakage.
// - An FFT is performed on each windowed segment to obtain the frequency-domain
// representation.
// - The power spectrum of each segment is computed, and the average spectrum is
// accumulated across segments.
//
// 2. Filter Edge Adjustments:
//
// - Adjusts the frequency bounds (nfa and nfb) to ensure the analysis remains
// within valid and meaningful regions of the signal.
//
// 3. Baseline Computation:
//
// - The average spectrum is converted to a dB scale.
// - Baseline is computed to distinguish significant signal components from
// background noise.
//
// 4. Synchronization Metric Calculation:
//
// - For each frequency bin in the specified range, evaluates synchronization
// power using a Costas waveform.
// - Sync metric is computed over the index range, considering all combinations
// of Costas patterns.
// - The maximum sync value and its corresponding offset are recorded for each
// frequency bin.
//
// 5. Normalization:
//
// - The sync values are normalized to the 40th percentile value using a ranked
// index. This ensures a consistent scaling across different signals and noise
// levels.
//
// 6. Candidate Extraction:
//
// - Candidates with a strong sync metric (above a defined threshold) are extracted.
// - Near-duplicate candidates of lesser synchronization power, based on frequency
// proximity, are eliminated.
//
// 7. Output:
//
// - Returns a vector of the most promising signal candidates, sorted by their
// synchronization power. It's expected that these will be re-sorted by the
// caller into a desirable order, but synchronization power order facilitates
// debugging this function.
//
// Note: The Fortran version of this routine would normalize `s` at the end of this
// function, but I'm unsure why; nothing beyond this function references `s`,
// so it was effectively a somewhat expensive dead store. It's been eliminated
// in this version.
std::vector<Sync>
syncjs8(int nfa,
int nfb)
{
// Compute symbol spectra
savg.fill(0.0f);
for (int j = 0; j < Mode::NHSYM; ++j)
{
int const ia = j * Mode::NSTEP;
int const ib = ia + Mode::NFFT1;
if (ib > Mode::NMAX) break;
std::transform(dd.begin() + ia,
dd.begin() + ib,
nuttal.begin(),
reinterpret_cast<float *>(sd.data()),
std::multiplies<float>{});
fftwf_execute(plans[Plan::SD]);
// Compute power spectrum
for (int i = 0; i < Mode::NSPS; ++i)
{
auto const power = std::norm(sd[i]);
s[i][j] = power;
savg[i] += power;
}
}
// Filter edge sanity measures
int const nwin = nfb - nfa;
if (nfa < 100)
{
nfa = 100;
if (nwin < 100) nfb = nfa + nwin;
}
if (nfb > 4910)
{
nfb = 4910;
if (nwin < 100) nfa = nfb - nwin;
}
auto const ia = std::max(0, static_cast<int>(std::round(nfa / Mode::DF)));
auto const ib = static_cast<int>(std::round(nfb / Mode::DF));
// Convert average spectrum from power to db scale and compute
// baseline from it; baseline replaces average spectrum.
baselinejs8(ia, ib);
// Compute and populate the sync index.
sync.clear();
for (int i = ia; i <= ib; ++i)
{
float max_value = -std::numeric_limits<float>::infinity();
int max_index = -Mode::JZ;
for (int j = -Mode::JZ; j <= Mode::JZ; ++j)
{
std::array<std::array<float, 3>, 2> t{};
for (int p = 0; p < 3; ++p)
{
for (int n = 0; n < 7; ++n)
{
int const offset = j + Mode::JSTRT + NSSY * n + p * 36 * NSSY;
if (offset >= 0 && offset < Mode::NHSYM)
{
// Accumulate Costas pattern contributions.
t[0][p] += s[i + NFOS * Costas[p][n]][offset];
// Accumulate sum over all frequencies for this block.
for (int freq = 0; freq < 7; ++freq)
{
t[1][p] += s[i + NFOS * freq][offset];
}
}
}
}
// Compute sync metric over the index range. We are at the moment
// maintaining the Fortran summation methodology for compatibility
// testing; there are more efficient ways to do this, but IEEE 754
// addition is a touchy thing, so we'll need to ensure that any
// changes don't negatively affect result precision.
auto const compute_sync = [&t](int start, int end)
{
float tx = 0.0f;
float t0 = 0.0f;
for (int i = start; i <= end; ++i)
{
tx += t[0][i];
t0 += t[1][i];
}
return tx / ((t0 - tx) / 6.0f);
};
if (auto const sync_value = std::max({
compute_sync(0, 2),
compute_sync(0, 1),
compute_sync(1, 2)
}); sync_value > max_value)
{
max_value = sync_value;
max_index = j;
}
}
sync.emplace(Mode::DF * i,
Mode::TSTEP * (max_index + 0.5f),
max_value);
}
// If we found nothing, we're done here.
if (sync.empty()) return {};
// Access the sync indices.
auto & freqIndex = sync.get<Tag::Freq>();
auto & rankIndex = sync.get<Tag::Rank>();
auto & syncIndex = sync.get<Tag::Sync>();
// Normalize to the 40th percentile using the frequency index,
// which is stable under sync value mutation. One thing to note
// here is that the Fortran version didn't seem to reliably
// calculate the 40th percentile rank; sometimes high, other
// times low, infrequently actually the 40th percentile value.
// This method should be perfectly accurate in all cases.
auto const normalize =
[
sync = rankIndex.nth(rankIndex.size() * 4 / 10)->sync
]
(Sync & entry)
{
entry.sync /= sync;
};
for (auto it = freqIndex.begin();
it != freqIndex.end();
++it)
{
freqIndex.modify(it, normalize);
}
// Extract candidates.
std::vector<Sync> candidates;
for (auto it = syncIndex.begin();
it != syncIndex.end() && candidates.size() < NMAXCAND;
it = syncIndex.begin())
{
// Stop iteration if below threshold or invalid; as the
// index is sorted by sync, any subsequent entries will
// also be below the threshold or invalid.
if (it->sync < ASYNCMIN || std::isnan(it->sync)) break;
// Good value, relatively strong; save the candidate.
candidates.push_back(*it);
// Remove the candidate and any near-duplicates based
// on frequency. This invalidates `it`, so we reset it
// to the index begin in the loop increment condition.
freqIndex.erase(
freqIndex.lower_bound(it->freq - Mode::AZ),
freqIndex.upper_bound(it->freq + Mode::AZ));
}
return candidates;
}
// Returns the total synchronization power, which is a measure of how well
// the signal aligns with the Costas sequence after accounting for the
// frequency adjustment. Used to identify the best alignment for further
// decoding.
float
syncjs8d(int const i0,
float const delf)
{
constexpr float BASE_DPHI = TAU * (1.0f / (12000.0f / Mode::NDOWN));
// If delta frequency is non-zero, compute the frequency
// adjustment array, otherwise, use what'll be an identity
// transfrom when multiplied.
std::array<std::complex<float>, Mode::NDOWNSPS> freqAdjust;
if (delf != 0.0f)
{
float const dphi = BASE_DPHI * delf;
float phi = 0.0f;
// std::fmod() is almost like Fortran's mod(), but not quite;
// Since delf can be negative, we must ensure that phi stays
// within [0, TAU), which Fortran's mod() handles by itself.
for (int i = 0; i < Mode::NDOWNSPS; ++i)
{
freqAdjust[i] = std::polar(1.0f, phi);
if (phi = std::fmod(phi + dphi, TAU);
phi < 0.0f)
{
phi += TAU;
}
}
}
else
{
freqAdjust.fill(std::complex<float>{1.0f, 0.0f});
}
// Compute sync power by looping over the Costas indices for
// each of the 3 Costas blocks, accumulating as we go.
float sync = 0.0f;
for (int i = 0; i < 3; ++i)
{
for (int j = 0; j < 7; ++j)
{
if (auto const offset = 36 * i * Mode::NDOWNSPS
+ i0 + j * Mode::NDOWNSPS; offset >= 0 &&
offset + Mode::NDOWNSPS <= Mode::NP2)
{
sync += std::norm(
std::transform_reduce(
freqAdjust.begin(), // Range start
freqAdjust.end(), // Range end
cd0.begin() + offset, // Data start
std::complex<float>{}, // Initial reduction value
std::plus<>{}, // Reduction by accumulation
[&](auto const & fa, // Conjugate and multiply
auto const & cd)
{
return cd * std::conj(fa * csyncs[i][j][&fa - &freqAdjust[0]]);
}
));
}
}
}
return sync;
}
// Generate a reference signal, based on the provided tone sequence and
// base frequency. The output is a vector of complex values representing
// the signal in the time domain.
std::vector<std::complex<float>>
genjs8refsig(std::array<int, NN> const & itone,
float const f0)
{
// Precompute the base frequency contribution; full circle in
// radians, multipled by the base frequency, multiplied by the
// sampling interval, i.e., the time step between samples, which
// results in the base frequency phase increment. Start the
// phase accumulator off at zero.
float const BFPI = TAU * f0 * (1.0f / 12000.0f);
auto phi = 0.0f;
std::vector<std::complex<float>> cref;
cref.reserve(NN * Mode::NSPS);
for (int i = 0; i < NN; ++i)
{
// Compute phase increment for the tone; frequency offset is
// determined by the tone value.
float const dphi = BFPI + TAU * static_cast<float>(itone[i]) / Mode::NSPS;
// Iterate over the samples per symbol to generate the time
// domain signal.
for (std::size_t is = 0; is < Mode::NSPS; ++is)
{
cref.push_back(std::polar(1.0f, phi));
phi = std::fmod(phi + dphi, TAU);
}
}
return cref;
}
// Subtract a JS8 signal
//
// Measured signal : dd(t) = a(t)cos(2*pi*f0*t+theta(t))
// Reference signal : cref(t) = exp( j*(2*pi*f0*t+phi(t)) )
// Complex amp : cfilt(t) = LPF[ dd(t)*CONJG(cref(t)) ]
// Subtract : dd(t) = dd(t) - 2*REAL{cref*cfilt}
//
// Important to note that dt can be negative here.
void
subtractjs8(std::vector<std::complex<float>> const & cref,
float const dt)
{
auto const nstart = static_cast<int>(dt * 12000.0f);
std::size_t const cref_start = (nstart < 0) ? static_cast<std::size_t>(-nstart) : 0;
std::size_t const dd_start = (nstart > 0) ? static_cast<std::size_t>( nstart) : 0;
auto const size = std::min(cref.size() - cref_start, dd.size() - dd_start);
// Populate complex filter with the conjugate of the reference signal.
for (std::size_t i = 0; i < size; ++i)
{
cfilt[i] = dd[dd_start + i] * std::conj(cref[cref_start + i]);
}
// Zero-fill the remainder, if any.
std::fill(cfilt.begin() + size, cfilt.end(), ZERO);
// FFT to the frequency domain.
fftwf_execute(plans[Plan::CF]);
// Apply the filter in the frequency domain.
std::transform(cfilt.begin(),
cfilt.end(),
filter.begin(),
cfilt.begin(),
std::multiplies<>());
// Inverse FFT to return to the time domain.
fftwf_execute(plans[Plan::CB]);
// Subtract the reconstructed signal.
for (std::size_t i = 0; i < size; ++i)
{
dd[dd_start + i] -= 2.0f * std::real(cfilt[i] * cref[cref_start + i]);
}
}
public:
// Constructor
DecodeMode()
{
// Intialize the Nuttal window. In theory, we can do this as a
// constexpr function at compile time, but doing so yield results
// slightly different than the Fortran version did, so for sanity
// while testing, we'll opt for consistency. IEEE 754 is always a
// bit brittle.
constexpr float a0 = 0.3635819f;
constexpr float a1 = -0.4891775f;
constexpr float a2 = 0.1365995f;
constexpr float a3 = -0.0106411f;
// Computed Pi constant to match the Fortran version; we could
// probably use std::numbers::pi_v<float> here, but for the
// moment, matching Fortran exactly.
float const pi = 4.0f * std::atan(1.0f);
float sum = 0.0f;
for (std::size_t i = 0; i < nuttal.size(); ++i)
{
// Naive summation here will exhibit substantial precision loss
// relative to the Fortran version; we use Kahan summation to
// compensate, which should yield results identical to Fortran.
KahanSum value = a0;
value += a1 * std::cos(2 * pi * i / nuttal.size());
value += a2 * std::cos(4 * pi * i / nuttal.size());
value += a3 * std::cos(6 * pi * i / nuttal.size());
nuttal[i] = value;
sum += value;
}
// Normalize the Nuttal window.
for (auto & value : nuttal) value = value / sum * nuttal.size() / 300.0f;
// Initialize Costas waveforms.
for (int i = 0; i < 7; ++i)
{
float const dphia = TAU * Costas[0][i] / Mode::NDOWNSPS;
float const dphib = TAU * Costas[1][i] / Mode::NDOWNSPS;
float const dphic = TAU * Costas[2][i] / Mode::NDOWNSPS;
float phia = 0.0f;
float phib = 0.0f;
float phic = 0.0f;
for (int j = 0; j < Mode::NDOWNSPS; ++j)
{
csyncs[0][i][j] = std::polar(1.0f, phia);
csyncs[1][i][j] = std::polar(1.0f, phib);
csyncs[2][i][j] = std::polar(1.0f, phic);
phia = std::fmod(phia + dphia, TAU);
phib = std::fmod(phib + dphib, TAU);
phic = std::fmod(phic + dphic, TAU);
}
}
// Compute a Hann-like window directly into the real part of the
// first NFILT + 1 elements in the filter, accumulating the sum
// as we go.
sum = 0.0f;
for (int j = -NFILT / 2; j <= NFILT / 2; ++j)
{
int const index = j + NFILT / 2;
float const value = std::pow(std::cos(pi * j / NFILT), 2);
filter[index].real(value);
sum += value;
}
// Now that we've got the sum, create actual complex numbers using
// the normalized real values that we just populated and zero the
// rest of the filter.
std::fill(std::transform(filter.begin(),
filter.begin() + NFILT + 1,
filter.begin(),
[sum](auto const value)
{
return std::complex<float>(value.real() / sum, 0.0f);
}),
filter.end(),
ZERO);
// Shift to position the window.
std::rotate(filter.begin(),
filter.begin() + NFILT / 2,
filter.begin() + NFILT + 1);
// Transform the filter into the frequency domain.
fftwf_plan fftw_plan;
{
std::lock_guard<std::mutex> lock(fftw_mutex);
fftw_plan = fftwf_plan_dft_1d(Mode::NMAX,
reinterpret_cast<fftwf_complex *>(filter.data()),
reinterpret_cast<fftwf_complex *>(filter.data()),
FFTW_FORWARD,
FFTW_ESTIMATE_PATIENT);
if (!fftw_plan)
{
throw std::runtime_error("Failed to create FFT plan");
}
}
fftwf_execute(fftw_plan);
{
std::lock_guard<std::mutex> lock(fftw_mutex);
fftwf_destroy_plan(fftw_plan);
}
// Normalize the frequency domain representation.
std::transform(filter.begin(),
filter.end(),
filter.begin(),
[factor = 1.0f / Mode::NMAX](auto value)
{
return value * factor;
});
// The rest of our FFT plans are always the same size and operate on the
// same data, so we can reuse them as long as we're alive.
std::lock_guard<std::mutex> lock(fftw_mutex);
plans[Plan::DS] = fftwf_plan_dft_1d(Mode::NDFFT2,
reinterpret_cast<fftwf_complex *>(cd0.data()),
reinterpret_cast<fftwf_complex *>(cd0.data()),
FFTW_BACKWARD,
FFTW_ESTIMATE_PATIENT);
plans[Plan::BB] = fftwf_plan_dft_r2c_1d(Mode::NDFFT1,
reinterpret_cast<float *>(ds_cx.data()),
reinterpret_cast<fftwf_complex *>(ds_cx.data()),
FFTW_ESTIMATE_PATIENT);
plans[Plan::CF] = fftwf_plan_dft_1d(Mode::NMAX,
reinterpret_cast<fftwf_complex *>(cfilt.data()),
reinterpret_cast<fftwf_complex *>(cfilt.data()),
FFTW_FORWARD,
FFTW_ESTIMATE_PATIENT);
plans[Plan::CB] = fftwf_plan_dft_1d(Mode::NMAX,
reinterpret_cast<fftwf_complex *>(cfilt.data()),
reinterpret_cast<fftwf_complex *>(cfilt.data()),
FFTW_BACKWARD,
FFTW_ESTIMATE_PATIENT);
plans[Plan::SD] = fftwf_plan_dft_r2c_1d(Mode::NFFT1,
reinterpret_cast<float *>(sd.data()),
reinterpret_cast<fftwf_complex *>(sd.data()),
FFTW_ESTIMATE_PATIENT);
plans[Plan::CS] = fftwf_plan_dft_1d(Mode::NDOWNSPS,
reinterpret_cast<fftwf_complex *>(csymb.data()),
reinterpret_cast<fftwf_complex *>(csymb.data()),
FFTW_FORWARD,
FFTW_ESTIMATE_PATIENT);
for (auto plan : plans)
{
if (!plan) throw std::runtime_error("Failed to create FFT plan");
}
}
// Decode entry point.
std::size_t
operator()(struct dec_data const & data,
int const kpos,
int const ksz,
JS8::Event::Emitter emitEvent)
{
// Copy the relevant frames for decoding
auto const pos = std::max(0, kpos);
auto const sz = std::max(0, ksz);
assert(sz <= Mode::NMAX);
if (data.params.syncStats) emitEvent(JS8::Event::SyncStart{pos, sz});
auto const ddCopy = [](auto const begin,
auto const end,
auto const to)
{
std::transform(begin, end, to, [](auto const value)
{
return static_cast<float>(value);
});
};
dd.fill(0.0f);
if ((JS8_RX_SAMPLE_SIZE - pos) < sz)
{
// Wrap case; split into two parts.
int const firstsize = JS8_RX_SAMPLE_SIZE - pos;
int const secondsize = sz - firstsize;
ddCopy(std::begin(data.d2) + pos, std::begin(data.d2) + pos + firstsize, dd.begin());
ddCopy(std::begin(data.d2), std::begin(data.d2) + secondsize, dd.begin() + firstsize);
}
else
{
// Non-wrapping case; copy directly.
ddCopy(std::begin(data.d2) + pos, std::begin(data.d2) + pos + sz, dd.begin());
}
Decode::Map decodes;
for (int ipass = 1; ipass <= 3; ++ipass)
{
// Determine if there's anything worth considering in the signal.
// If not, then we can just bail completely; more passes will not
// yield more results. If we do have some candidates, sort them
// by frequency, but put any that are close to nfqso up front.
auto candidates = syncjs8(data.params.nfa,
data.params.nfb);
if (candidates.empty()) break;
std::sort(candidates.begin(),
candidates.end(),
[nfqso = data.params.nfqso](auto const & a,
auto const & b)
{
auto const a_dist = std::abs(a.freq - nfqso);
auto const b_dist = std::abs(b.freq - nfqso);
if (a_dist < 10.0f && b_dist >= 10.0f) return true;
if (b_dist < 10.0f && a_dist >= 10.0f) return false;
return std::tie(a_dist, a.freq) <
std::tie(b_dist, b.freq);
});
// Recompute the baseband signal; subtraction during the last
// pass might have changed the landscape.
computeBasebandFFT();
bool const subtract = ipass < 3;
bool improved = false;
for (auto [f1, xdt, sync] : candidates)
{
float xsnr = 0.0f;
int nharderrors = -1;
if (auto decode = js8dec(data.params.syncStats,
subtract,
f1,
xdt,
nharderrors,
xsnr,
emitEvent))
{
// We don't need to be emitting duplicate events for something
// that's effectively the same SNR as a previous event.
auto const snr = static_cast<int>(std::round(xsnr));
// If this decode is new, or it's a duplicate with a better SNR
// than what we had before, then our situation has improved and
// we must announce that we've had some success.
if (auto [it, inserted] = decodes.try_emplace(std::move(*decode), snr);
inserted || it->second < snr)
{
improved = true;
// Update the SNR if this is an improved decode.
if (!inserted) it->second = snr;
// Emit decoded events on new or improved decodes.
emitEvent(JS8::Event::Decoded{data.params.nutc,
snr,
xdt - Mode::ASTART,
f1,
it->first.data,
it->first.type,
1.0f - nharderrors / 60.0f,
Mode::NSUBMODE});
}
}
}
// If nothing from this pass improved our situation, there's no
// point in trying any remaining passes.
if (!improved) break;
}
// Let the caller know how many unique decodes we discovered, if any.
return decodes.size();
}
};
// Explicit template class instantiations; avoids compiler complaints
// about unused variables.
template class DecodeMode<ModeA>;
template class DecodeMode<ModeB>;
template class DecodeMode<ModeC>;
template class DecodeMode<ModeE>;
template class DecodeMode<ModeI>;
}
/******************************************************************************/
// Worker
/******************************************************************************/
namespace JS8
{
class Worker : public QObject
{
Q_OBJECT
// Initialization of the decoders, in that they're heavy with
// FFT plan creations, is non-trivial, so a handle-body class
// to avoid initializing them on the main thread.
class Impl
{
// To avoid data races, decode data is referenced here but is
// actually located in the Worker that instantiates us, as it
// must be possible to copy data for us before we're ready to
// process it.
struct dec_data & m_data;
// Mode-specific decode strategy; we'll instantiate one of
// these for each of the 5 modes; this class is an aggregate
// of the 5 modes.
struct DecodeEntry
{
std::variant<
DecodeMode<ModeA>,
DecodeMode<ModeB>,
DecodeMode<ModeC>,
DecodeMode<ModeE>,
DecodeMode<ModeI>
> decode;
int mode;
int & kpos;
int & ksz;
template <typename DecodeModeType>
DecodeEntry(std::in_place_type_t<DecodeModeType>,
int mode,
int & kpos,
int & ksz)
: decode(std::in_place_type<DecodeModeType>)
, mode (mode)
, kpos (kpos)
, ksz (ksz)
{}
};
// Since a strategy can be neither moved nor copied, we must
// instantiate them in-place. Note that with the advent of the
// multi-decoder, mode identifiers became a bitset instead of
// integral values. The order defined here is the order that
// the decode loop will run in; we're matching the Fortran
// version here in terms of faster modes first.
template <typename ModeType>
DecodeEntry makeDecodeEntry(int shift,
int & kpos,
int & ksz)
{
return DecodeEntry(std::in_place_type<DecodeMode<ModeType>>,
1 << shift,
kpos,
ksz);
}
std::array<DecodeEntry, 5> m_decodes =
{{
makeDecodeEntry<ModeI>(4, m_data.params.kposI, m_data.params.kszI),
makeDecodeEntry<ModeE>(3, m_data.params.kposE, m_data.params.kszE),
makeDecodeEntry<ModeC>(2, m_data.params.kposC, m_data.params.kszC),
makeDecodeEntry<ModeB>(1, m_data.params.kposB, m_data.params.kszB),
makeDecodeEntry<ModeA>(0, m_data.params.kposA, m_data.params.kszA)
}};
public:
// Constructor
explicit Impl(struct dec_data & data)
: m_data(data)
{}
// Execute a decoding pass, using the supplied event emitter to
// emit events as they occur.
void operator()(Event::Emitter emitEvent)
{
// The multi-decoder can provide data for multiple modes at
// the same time; specific decodes to be performed for this
// pass are in the `nsubmodes` bitset.
auto const set = m_data.params.nsubmodes;
std::size_t sum = 0;
// Let any interested parties know that we've started a run
// for the set of modes requested.
emitEvent(Event::DecodeStarted{set});
// Iterate through all the modes we're aware of, performing
// a mode-specific decode pass if the mode is scheduled for
// decoding during this pass.
for (auto & entry : m_decodes)
{
if ((set & entry.mode) == entry.mode)
{
std::visit([&](auto && decode) {
sum += decode(m_data,
entry.kpos,
entry.ksz,
emitEvent);
}, entry.decode);
}
}
// Let any interested parties know the total number of decodes
// performed during this run.
emitEvent(Event::DecodeFinished{sum});
}
};
// Data members
QSemaphore * m_semaphore;
std::atomic<bool> m_quit = false;
struct dec_data m_data;
public:
// Constructor
explicit Worker(QSemaphore * semaphore,
QObject * parent = nullptr)
: QObject (parent)
, m_semaphore(semaphore)
{}
// Used to inform the worker that it's time to go; the next
// time it wakes up due to the semaphore being released, it
// will exit the runloop.
void stop()
{
m_quit = true;
}
// Called by the owning Decoder to refresh the copy of the
// decode data that the Worker implementation references.
void copy()
{
m_data = dec_data;
};
signals:
// Signal used to indicate that something of interest has
// occurred during a decoding pass.
void decodeEvent(Event::Variant const &);
public slots:
// Runloop for the thread that the worker is scheduled on; this
// is started by the Decoder when it's informed that the thread
// has started. Performs decoding runs each time the semaphore
// is released, until it's informed that it should quit.
void run()
{
// Our thread has started, and we're now running on it, so
// we're good to now allocate our implementation; we didn't
// want that to happen on the main thread, as the FFT plans
// can take a while. We only need the implementation while
// we're running.
std::unique_ptr<Impl> impl = std::make_unique<Impl>(m_data);
// Wait until there's something that requires our attention,
// which is going to either be needing to quit or needing to
// perform a decoding pass.
while (true)
{
m_semaphore->acquire();
if (m_quit) break;
(*impl)([this](Event::Variant const & event)
{
emit decodeEvent(event);
});
}
}
};
}
/******************************************************************************/
// Public Interface - Decoding
/******************************************************************************/
#include "JS8.moc"
namespace JS8
{
Decoder::Decoder(QObject * parent)
: QObject(parent)
, m_semaphore(0)
, m_worker(new Worker(&m_semaphore))
{
m_worker->moveToThread(&m_thread);
connect(&m_thread, &QThread::started, m_worker, &Worker::run);
connect(&m_thread, &QThread::finished, m_worker, &QObject::deleteLater);
connect(m_worker, &Worker::decodeEvent, this, &Decoder::decodeEvent);
}
void
Decoder::start(QThread::Priority priority)
{
m_thread.start(priority);
}
void
Decoder::quit()
{
m_worker->stop();
m_semaphore.release();
m_thread.quit();
m_thread.wait();
}
void
Decoder::decode()
{
m_worker->copy();
m_semaphore.release();
}
}
/******************************************************************************/
// Public Interface - Encoding
/******************************************************************************/
namespace JS8
{
// Port of the Fortran `genjs8` subroutine; from the 12 bytes of `message`,
// construct an 87-bit JS8 message and encode it into tones. Costas array
// to use supplied by the caller, as is the type of message, indicated by
// the lower 3 bits of `type`.
void
encode(int const type,
Costas::Array const & costas,
const char * const message,
int * const tones)
{
// Our initial goal here is an 87-bit message, for which a std::bitset
// would be the obvious choice, but we've got to compute a checksum of
// the first 75 bits; thus, an array instead.
//
// Message structure:
//
// +----------+----------+----------+
// | | | 72 bits | 12 6-bit words
// | | +==========+
// | | 87 bits | 3 bits | Frame type
// | 11 bytes | +==========+
// | | | 12 bits | 12-bit BE checksum
// | |----------+==========+
// | | 1 bit | 1 bit | Leftover bit in array
// +----------+----------+==========+
std::array<std::uint8_t, 11> bytes = {};
// Convert the 12 characters we've been handed to 6-bit words and pack
// them into the byte array, 4 characters, 24 bits at a time, into the
// 9 bytes [0,8], 72 bits total. Throws if handed an invalid character.
for (int i = 0, j = 0; i < 12; i += 4, j += 3)
{
std::uint32_t words = (alphabetWord(message[i ]) << 18) |
(alphabetWord(message[i + 1]) << 12) |
(alphabetWord(message[i + 2]) << 6) |
alphabetWord(message[i + 3]);
bytes[j ] = words >> 16;
bytes[j + 1] = words >> 8;
bytes[j + 2] = words;
}
// The bottom 3 bits of type are the frame type; these go into the
// next 3 bits in the byte array, i.e., the first 3 bits of byte 9,
// after which we'll be at 75 bits in total.
bytes[9] = (type & 0b111) << 5;
// We now need to compute the augmented CRC-12 of the complete
// byte array, including the trailing zero bits that we've not
// set yet.
auto const crc = CRC12(bytes);
// That CRC needs to occupy the next 12 bits of the array, i.e.,
// the final 5 bits of byte 9, and the first 7 bits of byte 10.
bytes[9] |= (crc >> 7) & 0x1F;
bytes[10] = (crc & 0x7F) << 1;
// That's it for our 87-bit message; we're now going to turn it
// into two blocks of 29 3-bit words, which will in turn become
// tones, the first block being parity for the second, bracketed
// by the Costas arrays.
//
// Output structure:
//
// +----------+----------+
// | | 7 bytes | Costas array A
// | +==========+
// | | 29 bytes | Parity data
// | +==========+
// | 79 bytes | 7 bytes | Costas array B
// | +==========+
// | | 29 bytes | Output data
// | +==========+
// | | 7 bytes | Costas array C
// +----------+==========+
auto costasData = tones;
auto parityData = tones + 7;
auto outputData = tones + 43;
// Output the 3 Costas arrays at offsets 0, 36, and 72.
for (auto const & array : costas)
{
std::copy(array.begin(), array.end(), costasData);
costasData += 36;
}
// Our 87 bits are going to be morphed into two sets of 29 3-bit
// words, the first one parity for the second; we're going to do
// this in parallel.
std::size_t outputBits = 0;
std::size_t outputByte = 0;
std::uint8_t outputMask = 0x80;
std::uint8_t outputWord = 0;
std::uint8_t parityWord = 0;
for (std::size_t i = 0; i < 87; ++i)
{
// Compute parity for the current bit; inputs for parity computation
// are the corresponding parity matrix row and each bit in the message;
// the parity matrix row, referenced by `i`, contains 87 boolean values.
// Each `true` value defines a message bit that must be summed, modulo
// 2, to produce the parity check bit for the bit we're working on now.
//
// In short, if the parity matrix bit `(i, j)` and the message bit `j`
// are both set, then we add 1 to the parity bits accumulator. If, after
// processing all message bits the accumulated result is odd, then the
// parity bit should be set for the current bit.
std::size_t parityBits = 0;
std::size_t parityByte = 0;
std::uint8_t parityMask = 0x80;
for (std::size_t j = 0; j < 87; ++j)
{
parityBits += parity(i, j) && (bytes[parityByte] & parityMask);
parityMask = (parityMask == 1) ? (++parityByte, 0x80) : (parityMask >> 1);
}
// Accumulate the parity and output bits; this is the point at which
// we perform the modulo 2 operation on the summed parity bits.
parityWord = (parityWord << 1) | (parityBits & 1);
outputWord = (outputWord << 1) | ((bytes[outputByte] & outputMask) != 0);
outputMask = (outputMask == 1) ? (++outputByte, 0x80) : (outputMask >> 1);
// If we're at a 3-bit boundary, output the words and reset.
if (++outputBits == 3)
{
*parityData++ = parityWord;
*outputData++ = outputWord;
parityWord = 0;
outputWord = 0;
outputBits = 0;
}
}
}
}
/******************************************************************************/
Reference in a new issue View git blame Copy permalink