Merge branch 'vladisslav2011-b1c_enable_sumproduct' into next

This commit is contained in:
Carles Fernandez
2026-09-30 15:07:13 +02:00
16 changed files with 722 additions and 508 deletions
@@ -124,32 +124,6 @@ bool build_rotated_bpsk_frame_from_iq(
}
bool apply_llr_scale_and_decode(
std::vector<float>& frame,
Beidou_Cnav1_Navigation_Message& nav,
int32_t expected_prn,
float cn0_db_hz,
int32_t* fail_stage)
{
const int32_t n = BEIDOU_CNAV1_FRAME_SYMBOLS;
double rms_acc = 0.0;
for (float sample : frame)
{
const auto v = static_cast<double>(sample);
rms_acc += v * v;
}
const double rms = std::sqrt(rms_acc / static_cast<double>(n)) + 1e-6;
const double cn0_linear = std::pow(10.0, static_cast<double>(cn0_db_hz) / 10.0);
const double snr_sym = std::max(cn0_linear * (BEIDOU_B1C_CODE_PERIOD_MS * 1e-3), 1e-4);
const double llr_scale = std::min(std::max(std::sqrt(snr_sym) / rms, 0.05), 50.0);
for (float& sample : frame)
{
sample = static_cast<float>(static_cast<double>(sample) * llr_scale);
}
return nav.decode_frame_symbols(frame.data(), BEIDOU_CNAV1_FRAME_SYMBOLS, expected_prn, fail_stage);
}
bool probe_frame_from_iq_window(
const std::vector<float>& history_i,
const std::vector<float>& history_q,
@@ -174,7 +148,6 @@ bool decode_frame_from_iq_window(
int32_t expected_prn,
int32_t start_offset,
bool invert,
float cn0_db_hz,
int32_t* fail_stage = nullptr)
{
std::vector<float> frame;
@@ -182,7 +155,8 @@ bool decode_frame_from_iq_window(
{
return false;
}
return apply_llr_scale_and_decode(frame, nav, expected_prn, cn0_db_hz, fail_stage);
// LLR scaling is done in Beidou_Cnav1_Navigation_Message::decode_frame().
return nav.decode_frame_symbols(frame.data(), BEIDOU_CNAV1_FRAME_SYMBOLS, expected_prn, fail_stage);
}
@@ -475,7 +449,6 @@ int beidou_b1c_telemetry_decoder_gs::general_work(
// Exact secondary matches need no lag search; loose matches get a small span.
const int32_t local_offset_span =
(sec_candidates.first >= 1799.0F) ? 0 : B1C_SCAN_OFFSET_SPAN_LOOSE;
const auto cn0_db_hz = static_cast<float>(current_symbol.CN0_dB_hz);
int32_t full_decode_attempts = 0;
if (!candidate_offsets.empty())
@@ -506,7 +479,7 @@ int beidou_b1c_telemetry_decoder_gs::general_work(
++full_decode_attempts;
if (decode_frame_from_iq_window(
history_i, history_q, d_nav, expected_prn_scan, start_offset, inv,
cn0_db_hz, &fail_stage))
&fail_stage))
{
frame_decoded = true;
invert = inv;
@@ -26,6 +26,7 @@
#include "tow_to_trk.h"
#include <pmt/pmt.h>
#include <pmt/pmt_sugar.h>
#include <algorithm>
#include <array>
#include <cmath>
#include <iomanip>
@@ -46,7 +47,9 @@ constexpr float SECOND_CODE[BEIDOU_B2A_SECONDARY_CODE_LENGTH] = {1.0F, 1.0F, 1.0
// First CRC failure at the expected frame boundary drops sync immediately so
// the next 1 ms samples resume a sliding preamble search. Waiting for several
// consecutive CRC failures (the B1C long-file failure mode) would keep
// publishing a coasted TOW on a slipped alignment.
// publishing a coasted TOW on a slipped alignment. With the soft two-frame
// preamble detector the search re-acquires the next preamble, so holding the
// alignment would not recover more frames.
} // namespace
@@ -71,15 +74,6 @@ beidou_b2a_telemetry_decoder_gs::beidou_b2a_telemetry_decoder_gs(const Gnss_Sate
d_tow_to_trk(conf.tow_to_trk)
{
configure_basic_outputs();
// kron(preamble antipodal, 5-chip data secondary code). 1 → −1, 0 → +1.
for (int32_t i = 0; i < BEIDOU_CNAV2_PREAMBLE_SYMBOLS; i++)
{
const float pbit = (BEIDOU_B2A_PREAMBLE_SYMBOLS_STR[i] == '1') ? -1.0F : 1.0F;
for (int32_t k = 0; k < BEIDOU_B2A_SECONDARY_CODE_LENGTH; k++)
{
d_preamble_ms[static_cast<size_t>(i * BEIDOU_B2A_SECONDARY_CODE_LENGTH + k)] = pbit * SECOND_CODE[k];
}
}
if (d_enable_navdata_monitor)
{
this->message_port_register_out(pmt::mp("Nav_msg_from_TLM"));
@@ -163,21 +157,24 @@ void beidou_b2a_telemetry_decoder_gs::publish_navigation(double cn0_db_hz)
}
int32_t beidou_b2a_telemetry_decoder_gs::preamble_correlation() const
bool beidou_b2a_telemetry_decoder_gs::preamble_detected() const
{
// Soft, scale-independent statistic on the 1 ms samples of the candidate
// frame start, combined with the preamble one frame earlier when available.
if (d_symbol_history.size() < static_cast<size_t>(BEIDOU_CNAV2_FRAME_MS))
{
return 0;
return false;
}
auto it = d_symbol_history.end() - BEIDOU_CNAV2_FRAME_MS;
int32_t corr = 0;
for (int32_t i = 0; i < BEIDOU_CNAV2_PREAMBLE_MS; i++)
std::array<float, BEIDOU_CNAV2_PREAMBLE_MS> current{};
std::array<float, BEIDOU_CNAV2_PREAMBLE_MS> previous{};
std::copy_n(d_symbol_history.end() - BEIDOU_CNAV2_FRAME_MS, BEIDOU_CNAV2_PREAMBLE_MS, current.begin());
const bool have_previous = d_symbol_history.size() >= static_cast<size_t>(2 * BEIDOU_CNAV2_FRAME_MS);
if (have_previous)
{
const float s = (*it >= 0.0F) ? 1.0F : -1.0F;
corr += static_cast<int32_t>(s * d_preamble_ms[static_cast<size_t>(i)]);
++it;
std::copy_n(d_symbol_history.end() - 2 * BEIDOU_CNAV2_FRAME_MS, BEIDOU_CNAV2_PREAMBLE_MS, previous.begin());
}
return corr;
return Beidou_Cnav2_Navigation_Message::preamble_detection_statistic(
current.data(), have_previous ? previous.data() : nullptr) >= BEIDOU_CNAV2_PREAMBLE_DETECTION_THRESHOLD;
}
@@ -187,14 +184,6 @@ bool beidou_b2a_telemetry_decoder_gs::try_decode_frame()
{
return false;
}
if (!d_flag_frame_sync)
{
const int32_t corr = preamble_correlation();
if (std::abs(corr) <= BEIDOU_CNAV2_PREAMBLE_CORR_THRESHOLD)
{
return false;
}
}
std::array<float, BEIDOU_CNAV2_FRAME_SYMBOLS> nav{};
auto it = d_symbol_history.end() - BEIDOU_CNAV2_FRAME_MS;
@@ -240,20 +229,11 @@ int beidou_b2a_telemetry_decoder_gs::general_work(
const bool history_ready =
d_symbol_history.size() >= static_cast<size_t>(BEIDOU_CNAV2_FRAME_MS);
bool try_now = false;
if (history_ready)
{
if (!d_flag_frame_sync)
{
try_now = true;
}
else if ((d_sample_counter > d_preamble_index) &&
((d_sample_counter - d_preamble_index) %
static_cast<uint64_t>(BEIDOU_CNAV2_FRAME_MS)) == 0U)
{
try_now = true;
}
}
const bool at_frame_boundary =
d_flag_frame_sync && (d_sample_counter > d_preamble_index) &&
((d_sample_counter - d_preamble_index) % static_cast<uint64_t>(BEIDOU_CNAV2_FRAME_MS)) == 0U;
const bool try_now = history_ready &&
(at_frame_boundary || (!d_flag_frame_sync && preamble_detected()));
if (try_now)
{
@@ -50,13 +50,12 @@ private:
beidou_b2a_telemetry_decoder_gs(const Gnss_Satellite& satellite, const Tlm_Conf& conf);
void publish_navigation(double cn0_db_hz);
bool try_decode_frame();
int32_t preamble_correlation() const;
bool preamble_detected() const;
Beidou_Cnav2_Navigation_Message d_nav;
Nav_Message_Packet d_nav_msg_packet;
Gnss_Satellite d_satellite;
std::deque<float> d_symbol_history;
std::array<float, BEIDOU_CNAV2_PREAMBLE_MS> d_preamble_ms{};
std::string d_dump_filename;
std::ofstream d_dump_file;
std::unique_ptr<Tlm_CRC_Stats> d_Tlm_CRC_Stats;
+1 -1
View File
@@ -38,7 +38,7 @@ constexpr int32_t BEIDOU_CNAV2_CRC_BITS = 24;
constexpr int32_t BEIDOU_CNAV2_DATA_BITS = BEIDOU_CNAV2_INFO_BITS - BEIDOU_CNAV2_CRC_BITS;
constexpr int32_t BEIDOU_CNAV2_MES_TYPE_BITS = 6;
constexpr int32_t BEIDOU_CNAV2_SOW_LSB_S = 3;
constexpr int32_t BEIDOU_CNAV2_PREAMBLE_CORR_THRESHOLD = 115;
constexpr double BEIDOU_CNAV2_PREAMBLE_DETECTION_THRESHOLD = 5.2; //!< on the normalized soft preamble statistic, about |N(0,1)| on noise
constexpr int32_t BEIDOU_CNAV2_MSG_EPH1 = 10;
constexpr int32_t BEIDOU_CNAV2_MSG_EPH2 = 11;
@@ -75,23 +75,23 @@ const BeidouCnav1LdpcGraph& beidou_cnav1_ldpc_graph_88_44()
bool beidou_cnav1_ldpc_decode_200_100(const float* symbol_llr, int32_t num_bits, uint8_t* info_bits600)
{
return beidou_ldpc_decode(beidou_cnav1_ldpc_graph_200_100(), symbol_llr, num_bits, info_bits600, nullptr);
return beidou_ldpc_decode(beidou_cnav1_ldpc_graph_200_100(), symbol_llr, num_bits, info_bits600, nullptr, true);
}
bool beidou_cnav1_ldpc_decode_88_44(const float* symbol_llr, int32_t num_bits, uint8_t* info_bits264)
{
return beidou_ldpc_decode(beidou_cnav1_ldpc_graph_88_44(), symbol_llr, num_bits, info_bits264, nullptr);
return beidou_ldpc_decode(beidou_cnav1_ldpc_graph_88_44(), symbol_llr, num_bits, info_bits264, nullptr, true);
}
bool beidou_cnav1_ldpc_decode_200_100_codeword(const float* symbol_llr, int32_t num_bits, uint8_t* codeword_bits1200)
{
return beidou_ldpc_decode(beidou_cnav1_ldpc_graph_200_100(), symbol_llr, num_bits, nullptr, codeword_bits1200);
return beidou_ldpc_decode(beidou_cnav1_ldpc_graph_200_100(), symbol_llr, num_bits, nullptr, codeword_bits1200, true);
}
bool beidou_cnav1_ldpc_decode_88_44_codeword(const float* symbol_llr, int32_t num_bits, uint8_t* codeword_bits528)
{
return beidou_ldpc_decode(beidou_cnav1_ldpc_graph_88_44(), symbol_llr, num_bits, nullptr, codeword_bits528);
return beidou_ldpc_decode(beidou_cnav1_ldpc_graph_88_44(), symbol_llr, num_bits, nullptr, codeword_bits528, true);
}
@@ -47,11 +47,20 @@ void append_bits_as_chars(std::string& out, const uint8_t* bits, int32_t n)
}
// Minimum normalized soft-correlation margin between the best and second-best
// BCH codeword, (C1 - C2) / sum|y|. It keeps the rate of wrongly decoded SOH
// words below that of hard-decision bounded-distance decoding at t = 11, while
// rejecting fewer words (simulated over AWGN, Es/N0 from -3 to -1 dB).
constexpr float BCH_MIN_SOFT_MARGIN = 0.1F;
// ICD Table 7-2: SOH ranges from 0 to 3582 s in steps of 18 s.
constexpr uint32_t BEIDOU_CNAV1_SOH_MAX = 3582U / BEIDOU_CNAV1_SOH_LSB_S;
struct BchCodebook
{
int32_t n = 0;
int32_t k = 0;
int32_t corr_threshold = 0;
std::vector<int8_t> bipolar_codewords;
std::vector<uint16_t> messages;
};
@@ -94,13 +103,11 @@ std::vector<int8_t> encode_bch_hypothesis_bipolar(
BchCodebook build_bch_codebook(
int32_t n,
int32_t k,
const std::vector<int32_t>& feedback_pos_1based,
int32_t corr_threshold)
const std::vector<int32_t>& feedback_pos_1based)
{
BchCodebook cb{};
cb.n = n;
cb.k = k;
cb.corr_threshold = corr_threshold;
const uint32_t num_messages = 1U << static_cast<uint32_t>(k);
cb.bipolar_codewords.resize(static_cast<size_t>(num_messages) * static_cast<size_t>(n));
cb.messages.resize(num_messages);
@@ -115,27 +122,44 @@ BchCodebook build_bch_codebook(
}
bool decode_bch_codeword(const int32_t* in, int32_t* out_msg, const BchCodebook& cb)
// Soft-decision maximum-likelihood decoding: pick the codeword with the largest
// correlation with the soft symbols (positive values favor bit 1, i.e. bipolar -1).
// With min_margin > 0, reject words whose best and second-best correlations are
// closer than min_margin * sum|y|. The metric is independent of the input scale.
bool decode_bch_soft(const float* in, int32_t* out_msg, const BchCodebook& cb, float min_margin)
{
int32_t best_corr = std::numeric_limits<int32_t>::min();
double magnitude = 0.0;
for (int32_t j = 0; j < cb.n; j++)
{
magnitude += std::abs(static_cast<double>(in[j]));
}
if (!(magnitude > 0.0) || !std::isfinite(magnitude))
{
return false;
}
double best_corr = -std::numeric_limits<double>::max();
double second_corr = -std::numeric_limits<double>::max();
uint16_t best_message = 0U;
for (size_t i = 0; i < cb.messages.size(); i++)
{
int32_t corr = 0;
double corr = 0.0;
const auto* cw = cb.bipolar_codewords.data() + static_cast<ptrdiff_t>(i) * cb.n;
for (int32_t j = 0; j < cb.n; j++)
{
const int32_t bit = in[j] > 0 ? 1 : 0;
const int32_t bipolar = (bit == 0) ? 1 : -1;
corr += static_cast<int32_t>(cw[j]) * bipolar;
corr -= static_cast<double>(cw[j]) * static_cast<double>(in[j]);
}
if (corr > best_corr)
{
second_corr = best_corr;
best_corr = corr;
best_message = cb.messages[i];
}
else if (corr > second_corr)
{
second_corr = corr;
}
}
if (best_corr < cb.corr_threshold)
if (best_corr - second_corr < static_cast<double>(min_margin) * magnitude)
{
return false;
}
@@ -147,19 +171,31 @@ bool decode_bch_codeword(const int32_t* in, int32_t* out_msg, const BchCodebook&
}
static bool decode_bch_21_6(const int32_t* in, int32_t* out_msg)
// PRN word. A wrong decode cannot pass an expected-PRN check, so the margin
// test is only needed when the PRN is not known in advance.
bool decode_bch_21_6(const float* in, int32_t* out_msg, bool prn_known)
{
static const std::vector<int32_t> feedback_pos = {2, 4, 5, 6};
static const BchCodebook cb = build_bch_codebook(21, 6, feedback_pos, 15);
return decode_bch_codeword(in, out_msg, cb);
static const BchCodebook cb = build_bch_codebook(21, 6, feedback_pos);
return decode_bch_soft(in, out_msg, cb, prn_known ? 0.0F : BCH_MIN_SOFT_MARGIN);
}
static bool decode_bch_51_8(const int32_t* in, int32_t* out_msg)
// SOH word. It feeds the TOW directly, so always require the margin.
bool decode_bch_51_8(const float* in, int32_t* out_msg)
{
static const std::vector<int32_t> feedback_pos = {1, 4, 5, 6, 7, 8};
static const BchCodebook cb = build_bch_codebook(51, 8, feedback_pos, 29);
return decode_bch_codeword(in, out_msg, cb);
static const BchCodebook cb = build_bch_codebook(51, 8, feedback_pos);
if (!decode_bch_soft(in, out_msg, cb, BCH_MIN_SOFT_MARGIN))
{
return false;
}
uint32_t soh = 0U;
for (int32_t bit = 0; bit < 8; bit++)
{
soh = (soh << 1U) | static_cast<uint32_t>(out_msg[bit] & 1);
}
return soh <= BEIDOU_CNAV1_SOH_MAX;
}
@@ -613,26 +649,41 @@ bool Beidou_Cnav1_Navigation_Message::decode_frame(
return false;
}
// Soft symbols arrive with an arbitrary tracking correlator gain (and, from
// some telemetry paths, a heuristic pre-scaling). The LDPC decoder expects
// bit LLRs on a known scale, so normalize the mean magnitude to 4, as done
// for B-CNAV2. This matches 2|y|/sigma^2 near the decoding threshold.
constexpr double mean_llr_magnitude = 4.0;
double magnitude_sum = 0.0;
for (int32_t i = 0; i < BEIDOU_CNAV1_FRAME_SYMBOLS; i++)
{
if (!std::isfinite(symbols[i]))
{
set_fail(0);
return false;
}
magnitude_sum += std::abs(static_cast<double>(symbols[i]));
}
if (magnitude_sum <= 0.0)
{
set_fail(0);
return false;
}
const double llr_scale = mean_llr_magnitude * BEIDOU_CNAV1_FRAME_SYMBOLS / magnitude_sum;
std::array<float, BEIDOU_CNAV1_FRAME_SYMBOLS> bit_llr{};
for (int32_t i = 0; i < BEIDOU_CNAV1_FRAME_SYMBOLS; i++)
{
bit_llr[i] = symbols[i];
}
std::array<int32_t, BEIDOU_CNAV1_SUBFRAME1_SYMBOLS> hard_bits{};
for (int32_t i = 0; i < BEIDOU_CNAV1_SUBFRAME1_SYMBOLS; i++)
{
hard_bits[i] = (symbols[i] >= 0.0F) ? 1 : 0;
bit_llr[i] = static_cast<float>(static_cast<double>(symbols[i]) * llr_scale);
}
int32_t prn_bits[6];
int32_t soh_bits[8];
if (!decode_bch_21_6(hard_bits.data(), prn_bits))
if (!decode_bch_21_6(bit_llr.data(), prn_bits, expected_prn > 0))
{
set_fail(1);
return false;
}
if (!decode_bch_51_8(hard_bits.data() + 21, soh_bits))
if (!decode_bch_51_8(bit_llr.data() + 21, soh_bits))
{
set_fail(2);
return false;
@@ -792,19 +843,13 @@ bool Beidou_Cnav1_Navigation_Message::probe_subframe1_prn(const float* symbols,
return false;
}
std::array<int32_t, BEIDOU_CNAV1_SUBFRAME1_SYMBOLS> hard_bits{};
for (int32_t i = 0; i < BEIDOU_CNAV1_SUBFRAME1_SYMBOLS; i++)
{
hard_bits[i] = (symbols[i] >= 0.0F) ? 1 : 0;
}
int32_t prn_bits[6];
int32_t soh_bits[8];
if (!decode_bch_21_6(hard_bits.data(), prn_bits))
if (!decode_bch_21_6(symbols, prn_bits, expected_prn > 0))
{
return false;
}
if (!decode_bch_51_8(hard_bits.data() + 21, soh_bits))
if (!decode_bch_51_8(symbols + 21, soh_bits))
{
return false;
}
@@ -15,10 +15,12 @@
*/
#include "beidou_cnav2_navigation_message.h"
#include "Beidou_B2a.h"
#include "Beidou_CNAV1.h"
#include "Beidou_CNAV2.h"
#include "MATH_CONSTANTS.h"
#include "beidou_cnav2_ldpc.h"
#include <algorithm>
#include <cmath>
#include <cstdint>
@@ -98,25 +100,26 @@ bool Beidou_Cnav2_Navigation_Message::decode_frame_symbols(const float* symbols,
return false;
}
int corr = 0;
// Resolve carrier polarity from the soft preamble correlation. Positive NAV
// symbols are bit 0. No exact preamble match is required: frame alignment
// is validated by the LDPC parity checks and the CRC, and a hard 24/24
// match rejects most frames near the LDPC decoding threshold.
double corr = 0.0;
for (int32_t i = 0; i < BEIDOU_CNAV2_PREAMBLE_SYMBOLS; i++)
{
if (!std::isfinite(symbols[i]))
{
return false;
}
// Positive NAV symbol is bit 0; negative is bit 1.
const bool bit = symbols[i] < 0.0F;
corr += (bit == (PREAMBLE[i] == '1')) ? 1 : -1;
corr += (PREAMBLE[i] == '1') ? -static_cast<double>(symbols[i]) : static_cast<double>(symbols[i]);
}
if (std::abs(corr) != BEIDOU_CNAV2_PREAMBLE_SYMBOLS)
if (corr == 0.0)
{
return false;
}
// Resolve carrier polarity from the preamble and preserve soft magnitudes.
// The shared LDPC decoder uses the opposite sign: positive LLR favors bit 1.
const float polarity = corr < 0 ? 1.0F : -1.0F;
const float polarity = corr < 0.0 ? 1.0F : -1.0F;
std::array<float, BEIDOU_CNAV2_LDPC_SYMBOLS> bit_llr{};
for (int32_t i = 0; i < BEIDOU_CNAV2_LDPC_SYMBOLS; i++)
{
@@ -172,6 +175,49 @@ bool Beidou_Cnav2_Navigation_Message::decode_frame_symbols(const float* symbols,
}
double Beidou_Cnav2_Navigation_Message::preamble_detection_statistic(const float* current_1ms, const float* previous_1ms)
{
auto correlate = [](const float* samples, double& corr, double& energy) {
corr = 0.0;
energy = 0.0;
for (int32_t i = 0; i < BEIDOU_CNAV2_PREAMBLE_SYMBOLS; i++)
{
const double preamble_chip = (PREAMBLE[i] == '1') ? -1.0 : 1.0;
for (int32_t k = 0; k < BEIDOU_B2A_SECONDARY_CODE_LENGTH; k++)
{
const double chip = (BEIDOU_B2A_SECONDARY_CODE_STR[k] == '1') ? -preamble_chip : preamble_chip;
const auto y = static_cast<double>(samples[i * BEIDOU_B2A_SECONDARY_CODE_LENGTH + k]);
corr += chip * y;
energy += y * y;
}
}
};
if (current_1ms == nullptr)
{
return 0.0;
}
double corr = 0.0;
double energy = 0.0;
correlate(current_1ms, corr, energy);
if (!(energy > 0.0) || !std::isfinite(energy))
{
return 0.0;
}
double statistic = std::abs(corr) / std::sqrt(energy);
if (previous_1ms != nullptr)
{
double prev_corr = 0.0;
double prev_energy = 0.0;
correlate(previous_1ms, prev_corr, prev_energy);
if (std::isfinite(prev_energy))
{
statistic = std::max(statistic, std::abs(corr + prev_corr) / std::sqrt(energy + prev_energy));
}
}
return statistic;
}
void Beidou_Cnav2_Navigation_Message::reset()
{
d_published = {};
@@ -41,6 +41,19 @@ public:
uint32_t last_frame_prn() const { return d_last_frame_prn; }
const std::string& get_last_nav_bits() const { return d_last_nav_bits; }
/*!
* \brief Soft preamble detection statistic on 1 ms prompt samples.
*
* current_1ms points to BEIDOU_CNAV2_PREAMBLE_MS samples at a candidate frame
* start, previous_1ms (optional) to the samples one frame (3 s) earlier.
* With C the correlation with the preamble spread by the data secondary code
* and E the sample energy, it returns max(|C1| / sqrt(E1), |C1 + C2| / sqrt(E1 + E2)).
* The result does not depend on the input scale and is about |N(0,1)| on
* noise. The single-frame term keeps detection working across a carrier
* polarity flip between the two frames.
*/
static double preamble_detection_statistic(const float* current_1ms, const float* previous_1ms = nullptr);
private:
void parse_info_bits(const uint8_t* bits, uint32_t channel_prn);
void parse_clock_common(const uint8_t* bits, int32_t toc_off, int32_t iodc_off);
+103 -255
View File
@@ -27,8 +27,11 @@ namespace
constexpr int32_t BITS_PER_SYMBOL = 6;
constexpr int32_t Q = 64;
constexpr int32_t NM = BEIDOU_LDPC_NM;
constexpr int32_t FP_CANDIDATES = 8 + 2 * NM;
constexpr float LLR_FALLBACK_GAP = 25.0F;
// LLR assigned to field elements missing from a truncated C2V message is its
// largest kept LLR plus this offset (ICD Annex, EMS variable node rule).
// Tuned by simulation for bit LLRs on the 2|y|/sigma^2 scale (B-CNAV1
// LDPC(200,100) over AWGN, 1.5 to 3 dB Eb/N0).
constexpr float EMS_TRUNCATION_OFFSET = 2.0F;
struct TruncMessage
{
@@ -37,15 +40,6 @@ struct TruncMessage
};
struct FixedPathCandidate
{
FixedPathCandidate() = default;
FixedPathCandidate(uint8_t sym_, float llr_) : sym(sym_), llr(llr_) {}
uint8_t sym = 0;
float llr = std::numeric_limits<float>::max();
};
uint8_t bits_to_symbol_msb_first(const float* bit_llr, int32_t bit_offset)
{
uint8_t symbol = 0;
@@ -57,50 +51,33 @@ uint8_t bits_to_symbol_msb_first(const float* bit_llr, int32_t bit_offset)
}
void normalize_message(TruncMessage& msg)
{
const float min_llr = msg.llr[0];
for (int32_t i = 0; i < NM; i++)
{
msg.llr[i] -= min_llr;
}
}
float lookup_message_llr_var_domain(const TruncMessage& c2v_msg, uint8_t h_inv, uint8_t x_var)
{
for (int32_t i = 0; i < NM; i++)
{
const uint8_t mapped = GaloisField64::mul(h_inv, c2v_msg.sym[i]);
if (mapped == x_var)
{
return c2v_msg.llr[i];
}
}
return c2v_msg.llr[NM - 1] + LLR_FALLBACK_GAP;
}
// Keep the NM most reliable (smallest-LLR) distinct field elements of a
// full-alphabet metric vector, in ascending order, with the first LLR at zero.
void truncate_metric_vector(const std::array<float, Q>& metric, TruncMessage& out)
{
std::array<bool, Q> used{};
for (int32_t kept = 0; kept < NM; kept++)
int32_t kept = 0;
for (int32_t x = 0; x < Q; x++)
{
float best_llr = std::numeric_limits<float>::max();
int32_t best_sym = 0;
for (int32_t x = 0; x < Q; x++)
const float llr = metric[static_cast<size_t>(x)];
if (kept == NM && llr >= out.llr[NM - 1])
{
if (!used[static_cast<size_t>(x)] && metric[static_cast<size_t>(x)] < best_llr)
{
best_llr = metric[static_cast<size_t>(x)];
best_sym = x;
}
continue;
}
used[static_cast<size_t>(best_sym)] = true;
out.sym[kept] = static_cast<uint8_t>(best_sym);
out.llr[kept] = best_llr;
int32_t pos = (kept < NM) ? kept++ : NM - 1;
while (pos > 0 && out.llr[pos - 1] > llr)
{
out.llr[pos] = out.llr[pos - 1];
out.sym[pos] = out.sym[pos - 1];
pos--;
}
out.llr[pos] = llr;
out.sym[pos] = static_cast<uint8_t>(x);
}
const float min_llr = out.llr[0];
for (auto& llr : out.llr)
{
llr -= min_llr;
}
normalize_message(out);
}
@@ -200,17 +177,20 @@ void initialize_channel_llr(
}
// Extended Min-Sum (ICD Annex, section 2.(1)).
// V2C/C2V messages live in the edge domain y = h_ij * x_j, so a check node
// only needs GF(64) additions of the incoming symbols.
void initialize_v2c_messages(
const BeidouLdpcGraph& graph,
const std::vector<float>& channel_llr,
std::vector<TruncMessage>& v2c)
{
const int32_t edge_count = graph.num_checks * graph.row_weight;
std::array<float, Q> mapped_metric{};
for (int32_t edge = 0; edge < edge_count; edge++)
{
const uint16_t var = graph.check_to_var[static_cast<size_t>(edge)];
const uint8_t h_ij = graph.check_to_h[static_cast<size_t>(edge)];
std::array<float, Q> mapped_metric{};
for (int32_t x = 0; x < Q; x++)
{
const uint8_t y = GaloisField64::mul(h_ij, static_cast<uint8_t>(x));
@@ -222,243 +202,112 @@ void initialize_v2c_messages(
}
void variable_node_update(
const BeidouLdpcGraph& graph,
const std::vector<float>& channel_llr,
std::vector<TruncMessage>& v2c,
const std::vector<TruncMessage>& c2v,
std::vector<uint8_t>& hard_codeword)
// Elementary check-node operation V = U (+) W: every pair (d, delta) proposes
// the element Us[d] + Ws[delta] with LLR U[d] + W[delta]; keep the NM
// smallest distinct elements. With NM = 8 the full NM x NM matrix is cheap,
// so this is exact rather than a bubble-sort approximation.
void ems_combine(const TruncMessage& u, const TruncMessage& w, TruncMessage& out)
{
std::array<float, Q> posterior{};
std::array<float, Q> outgoing_metric{};
for (int32_t var = 0; var < graph.num_variables; var++)
std::array<float, Q> best{};
best.fill(std::numeric_limits<float>::max());
for (int32_t d = 0; d < NM; d++)
{
const uint32_t begin = graph.var_offsets[static_cast<size_t>(var)];
const uint32_t end = graph.var_offsets[static_cast<size_t>(var) + 1];
for (int32_t x = 0; x < Q; x++)
for (int32_t delta = 0; delta < NM; delta++)
{
float metric = channel_llr[static_cast<size_t>(var) * Q + static_cast<size_t>(x)];
for (uint32_t p = begin; p < end; p++)
const auto sym = static_cast<size_t>(u.sym[d] ^ w.sym[delta]);
const float llr = u.llr[d] + w.llr[delta];
if (llr < best[sym])
{
const uint32_t edge = graph.var_to_edge[p];
const uint8_t h_fj_inv = graph.var_to_h_inv[p];
metric += lookup_message_llr_var_domain(c2v[edge], h_fj_inv, static_cast<uint8_t>(x));
best[sym] = llr;
}
posterior[static_cast<size_t>(x)] = metric;
}
int32_t hard_symbol = 0;
float best_hard_llr = posterior[0];
for (int32_t x = 1; x < Q; x++)
{
if (posterior[static_cast<size_t>(x)] < best_hard_llr)
{
best_hard_llr = posterior[static_cast<size_t>(x)];
hard_symbol = x;
}
}
hard_codeword[static_cast<size_t>(var)] = static_cast<uint8_t>(hard_symbol);
for (uint32_t p = begin; p < end; p++)
{
const uint32_t edge_out = graph.var_to_edge[p];
const uint8_t h_ij = graph.var_to_h[p];
for (int32_t y = 0; y < Q; y++)
{
outgoing_metric[static_cast<size_t>(y)] = std::numeric_limits<float>::max();
}
for (int32_t x = 0; x < Q; x++)
{
float metric = channel_llr[static_cast<size_t>(var) * Q + static_cast<size_t>(x)];
for (uint32_t q = begin; q < end; q++)
{
const uint32_t edge_in = graph.var_to_edge[q];
if (edge_in == edge_out)
{
continue;
}
const uint8_t h_fj_inv = graph.var_to_h_inv[q];
metric += lookup_message_llr_var_domain(c2v[edge_in], h_fj_inv, static_cast<uint8_t>(x));
}
const uint8_t y_out = GaloisField64::mul(h_ij, static_cast<uint8_t>(x));
outgoing_metric[static_cast<size_t>(y_out)] = metric;
}
truncate_metric_vector(outgoing_metric, v2c[edge_out]);
}
}
// Fixing delta = 0 already yields NM distinct elements, so NM finite entries exist.
truncate_metric_vector(best, out);
}
void build_fixed_path_candidates(
const std::array<TruncMessage, 4>& in,
std::array<FixedPathCandidate, FP_CANDIDATES>& cand)
{
for (auto& c : cand)
{
c.sym = 0U;
c.llr = std::numeric_limits<float>::max();
}
const uint8_t s0 = in[0].sym[0];
const uint8_t s1 = in[1].sym[0];
const uint8_t s2 = in[2].sym[0];
const uint8_t s3 = in[3].sym[0];
const float r0 = in[0].llr[0];
const float r1 = in[1].llr[0];
const float r2 = in[2].llr[0];
const float r3 = in[3].llr[0];
const uint8_t base_sym = GaloisField64::add(GaloisField64::add(s0, s1), GaloisField64::add(s2, s3));
const float base_llr = r0 + r1 + r2 + r3;
int32_t idx = 0;
cand[idx++] = {base_sym, base_llr};
for (int32_t l = 0; l < 4; l++)
{
const uint8_t sym = GaloisField64::add(base_sym, GaloisField64::add(in[l].sym[0], in[l].sym[1]));
const float llr = base_llr - in[l].llr[0] + in[l].llr[1];
cand[idx++] = {sym, llr};
}
static const std::array<std::pair<int32_t, int32_t>, 6> pairs = {
std::pair<int32_t, int32_t>{0, 1},
{0, 2},
{0, 3},
{1, 2},
{1, 3},
{2, 3},
};
for (const auto& pair : pairs)
{
const int32_t a = pair.first;
const int32_t b = pair.second;
const uint8_t sym = GaloisField64::add(
base_sym,
GaloisField64::add(
GaloisField64::add(in[a].sym[0], in[a].sym[1]),
GaloisField64::add(in[b].sym[0], in[b].sym[1])));
const float llr = base_llr - in[a].llr[0] - in[b].llr[0] + in[a].llr[1] + in[b].llr[1];
cand[idx++] = {sym, llr};
}
for (int32_t rank = 2; rank < NM && idx < FP_CANDIDATES; rank++)
{
for (int32_t l = 0; l < 4 && idx < FP_CANDIDATES; l++)
{
const uint8_t sym = GaloisField64::add(base_sym, GaloisField64::add(in[l].sym[0], in[l].sym[rank]));
const float llr = base_llr - in[l].llr[0] + in[l].llr[rank];
cand[idx++] = {sym, llr};
}
}
}
void fixed_path_check_update(
// C2V_{i->j} = (+)_{j' in N(i), j' != j} V2C_{j'->i}, by forward-backward recursion.
void ems_check_update(
const BeidouLdpcGraph& graph,
const std::vector<TruncMessage>& v2c,
std::vector<TruncMessage>& c2v)
{
std::array<int32_t, 4> edge_idx{};
std::array<int32_t, 4> order{};
std::array<TruncMessage, 4> in_sorted{};
std::array<FixedPathCandidate, FP_CANDIDATES> candidates{};
std::array<uint8_t, FP_CANDIDATES> T{};
std::array<uint8_t, FP_CANDIDATES> Tbar{};
const int32_t dc = graph.row_weight;
std::vector<TruncMessage> forward(static_cast<size_t>(dc));
std::vector<TruncMessage> backward(static_cast<size_t>(dc));
for (int32_t check = 0; check < graph.num_checks; check++)
{
const uint32_t begin = graph.check_offsets[static_cast<size_t>(check)];
const uint32_t end = graph.check_offsets[static_cast<size_t>(check) + 1];
if (static_cast<int32_t>(end - begin) != 4)
forward[0] = v2c[begin];
for (int32_t l = 1; l < dc - 1; l++)
{
continue;
ems_combine(forward[static_cast<size_t>(l) - 1], v2c[begin + static_cast<uint32_t>(l)], forward[static_cast<size_t>(l)]);
}
for (int32_t l = 0; l < 4; l++)
backward[static_cast<size_t>(dc) - 1] = v2c[begin + static_cast<uint32_t>(dc) - 1];
for (int32_t l = dc - 2; l > 0; l--)
{
edge_idx[l] = static_cast<int32_t>(begin) + l;
order[l] = l;
ems_combine(v2c[begin + static_cast<uint32_t>(l)], backward[static_cast<size_t>(l) + 1], backward[static_cast<size_t>(l)]);
}
std::sort(order.begin(), order.end(), [&](int32_t a, int32_t b) {
return v2c[static_cast<size_t>(edge_idx[a])].llr[1] < v2c[static_cast<size_t>(edge_idx[b])].llr[1];
});
for (int32_t l = 0; l < 4; l++)
c2v[begin] = backward[1];
c2v[begin + static_cast<uint32_t>(dc) - 1] = forward[static_cast<size_t>(dc) - 2];
for (int32_t l = 1; l < dc - 1; l++)
{
in_sorted[l] = v2c[static_cast<size_t>(edge_idx[order[l]])];
ems_combine(forward[static_cast<size_t>(l) - 1], backward[static_cast<size_t>(l) + 1], c2v[begin + static_cast<uint32_t>(l)]);
}
}
}
build_fixed_path_candidates(in_sorted, candidates);
int32_t theta = 0;
int32_t beta = 1;
float min2 = in_sorted[0].llr[1];
float second2 = in_sorted[1].llr[1];
if (second2 < min2)
// V2C_{j->i} = h_ij * (sum_{f != i} C2V_{f->j} * h_fj^-1 + L_j)_NM.
// Elements missing from a truncated C2V get its largest LLR plus a fixed
// offset. The hard decision uses the full a-posteriori metric.
void ems_variable_update(
const BeidouLdpcGraph& graph,
const std::vector<float>& channel_llr,
const std::vector<TruncMessage>& c2v,
std::vector<TruncMessage>& v2c,
std::vector<uint8_t>& hard_codeword)
{
std::vector<std::array<float, Q>> incoming;
std::array<float, Q> posterior{};
std::array<float, Q> outgoing{};
for (int32_t var = 0; var < graph.num_variables; var++)
{
const uint32_t begin = graph.var_offsets[static_cast<size_t>(var)];
const uint32_t end = graph.var_offsets[static_cast<size_t>(var) + 1];
incoming.resize(end - begin);
for (int32_t x = 0; x < Q; x++)
{
std::swap(min2, second2);
std::swap(theta, beta);
posterior[static_cast<size_t>(x)] = channel_llr[static_cast<size_t>(var) * Q + static_cast<size_t>(x)];
}
for (int32_t l = 2; l < 4; l++)
for (uint32_t p = begin; p < end; p++)
{
const float v = in_sorted[l].llr[1];
if (v < min2)
const TruncMessage& msg = c2v[graph.var_to_edge[p]];
auto& in = incoming[p - begin];
in.fill(msg.llr[NM - 1] + EMS_TRUNCATION_OFFSET);
for (int32_t k = 0; k < NM; k++)
{
second2 = min2;
beta = theta;
min2 = v;
theta = l;
const uint8_t x = GaloisField64::mul(graph.var_to_h_inv[p], msg.sym[k]);
in[x] = std::min(in[x], msg.llr[k]);
}
else if (v < second2)
for (int32_t x = 0; x < Q; x++)
{
second2 = v;
beta = l;
posterior[static_cast<size_t>(x)] += in[static_cast<size_t>(x)];
}
}
const int32_t ref_idx = NM / 2;
const float theta_ref = in_sorted[theta].llr[ref_idx];
const float beta_ref = in_sorted[beta].llr[ref_idx];
for (int32_t k = 0; k < FP_CANDIDATES; k++)
hard_codeword[static_cast<size_t>(var)] =
static_cast<uint8_t>(std::min_element(posterior.begin(), posterior.end()) - posterior.begin());
for (uint32_t p = begin; p < end; p++)
{
T[k] = (candidates[k].llr <= theta_ref) ? 1U : 0U;
Tbar[k] = (candidates[k].llr <= beta_ref) ? 1U : 0U;
}
for (int32_t l = 0; l < 4; l++)
{
TruncMessage out{};
int32_t filled = 0;
for (int32_t k = 0; k < FP_CANDIDATES && filled < NM; k++)
const auto& in = incoming[p - begin];
const uint8_t h_ij = graph.var_to_h[p];
for (int32_t x = 0; x < Q; x++)
{
bool pass_gate = false;
if (l == 0)
{
pass_gate = (T[k] != 0U);
}
else if (l == theta)
{
pass_gate = (Tbar[k] != 0U);
}
else
{
pass_gate = (T[k] != 0U) || (Tbar[k] != 0U);
}
if (!pass_gate)
{
continue;
}
out.sym[filled] = GaloisField64::add(in_sorted[l].sym[0], candidates[k].sym);
out.llr[filled] = candidates[k].llr - in_sorted[l].llr[0];
filled++;
outgoing[GaloisField64::mul(h_ij, static_cast<uint8_t>(x))] =
posterior[static_cast<size_t>(x)] - in[static_cast<size_t>(x)];
}
while (filled < NM)
{
out.sym[filled] = out.sym[filled - 1];
out.llr[filled] = out.llr[filled - 1] + 1.0F;
filled++;
}
normalize_message(out);
c2v[static_cast<size_t>(edge_idx[order[l]])] = out;
truncate_metric_vector(outgoing, v2c[graph.var_to_edge[p]]);
}
}
}
@@ -504,7 +353,7 @@ void normalize_probability(ProbabilityMessage& message)
// Full-alphabet sum-product retry for blocks that defeat the truncated
// fixed-path decoder. Reuse the same graph, channel metrics and GF mapping.
// EMS decoder. Reuse the same graph, channel metrics and GF mapping.
bool sum_product_decode(const BeidouLdpcGraph& graph,
const std::vector<float>& channel_llr, std::vector<uint8_t>& hard_codeword)
{
@@ -699,10 +548,10 @@ bool decode_block(
std::vector<TruncMessage> v2c(static_cast<size_t>(edge_count));
std::vector<TruncMessage> c2v(static_cast<size_t>(edge_count));
initialize_v2c_messages(graph, channel_llr, v2c);
c2v = v2c;
for (int32_t iteration = 0; iteration < BEIDOU_LDPC_MAX_ITER; iteration++)
{
variable_node_update(graph, channel_llr, v2c, c2v, hard_codeword);
ems_check_update(graph, v2c, c2v);
ems_variable_update(graph, channel_llr, c2v, v2c, hard_codeword);
if (syndrome_is_zero(graph, hard_codeword))
{
if (info_bits != nullptr)
@@ -715,7 +564,6 @@ bool decode_block(
}
return true;
}
fixed_path_check_update(graph, v2c, c2v);
}
if (enable_sum_product && sum_product_decode(graph, channel_llr, hard_codeword))
+6 -3
View File
@@ -104,7 +104,7 @@ struct BeidouLdpcGraph
// ICD tables are laid out in four bands of four entries per printed row.
// Only row weight four is supported by the fixed-path check update.
// Only row weight four is supported (ICD table layout).
bool beidou_ldpc_init_graph(
int32_t num_checks,
int32_t num_variables,
@@ -118,8 +118,11 @@ bool beidou_ldpc_init_graph(
// Positive bit LLRs favor 1. Bits within each GF(64) symbol are MSB first.
// Supply at least one output: n-m information symbols or all n codeword symbols.
// Outputs are written only after the parity checks pass; CRC validation is separate.
// The optional full-alphabet retry adds up to 15 iterations after fixed-path failure.
// Disabled by default to preserve the existing B-CNAV1 decoding behavior.
// Bit LLRs should be calibrated (2|y|/sigma^2 for BPSK, as in the ICD Annex): the
// EMS truncation offset is tuned for that scale and sum-product assumes it.
// Decoding uses the Extended Min-Sum algorithm of the ICD Annex. The optional
// full-alphabet sum-product retry adds up to BEIDOU_LDPC_MAX_ITER iterations
// after an EMS failure.
bool beidou_ldpc_decode(const BeidouLdpcGraph& graph, const float* bit_llr,
int32_t num_bits, uint8_t* info_bits, uint8_t* codeword_bits, bool enable_sum_product = false);
@@ -15,9 +15,14 @@
*/
#include "Beidou_CNAV1.h"
#include "beidou_cnav1_ldpc.h"
#include "beidou_cnav_test_helpers.h"
#include <gtest/gtest.h>
#include <algorithm>
#include <array>
#include <cstdint>
#include <cstring>
#include <random>
#include <vector>
namespace
{
@@ -40,6 +45,40 @@ std::array<float, BEIDOU_CNAV1_SUBFRAME3_SYMBOLS> make_sf3_bit_llrs(const uint8_
}
return llr;
}
// Encode random information blocks, send them over BPSK/AWGN and count
// blocks the Extended Min-Sum decoder (no sum-product retry) fails to recover.
template <int N>
int32_t count_ems_block_errors(const BeidouCnav1LdpcGraph& graph, double ebn0_db, int32_t blocks, uint32_t seed)
{
BeidouCnavTest::PortableGaussian gauss(seed);
std::mt19937 bit_rng(seed + 1U);
const double sigma = BeidouCnavTest::rate_half_bpsk_sigma(ebn0_db);
int32_t errors = 0;
for (int32_t block = 0; block < blocks; block++)
{
std::array<uint8_t, N / 2> info{};
for (auto& bit : info)
{
bit = static_cast<uint8_t>(bit_rng() & 1U);
}
const auto codeword = BeidouCnavTest::encode<N>(info);
std::vector<float> llr(N);
for (int32_t i = 0; i < N; i++)
{
// Positive LLR favors bit 1; 2y/sigma^2 is the ICD Annex scale.
const double y = (codeword[static_cast<size_t>(i)] != 0U ? 1.0 : -1.0) + sigma * gauss();
llr[static_cast<size_t>(i)] = static_cast<float>(2.0 * y / (sigma * sigma));
}
std::vector<uint8_t> decoded(N, 0U);
if (!beidou_ldpc_decode(graph, llr.data(), N, nullptr, decoded.data(), false) ||
!std::equal(decoded.begin(), decoded.end(), codeword.begin()))
{
errors++;
}
}
return errors;
}
} // namespace
TEST(BeidouCnav1LdpcTest, DecodeZeroCodeword200_100)
@@ -101,3 +140,14 @@ TEST(BeidouCnav1LdpcTest, CorrectSingleSymbolError88_44)
EXPECT_EQ(bit, 0U);
}
}
TEST(BeidouCnav1LdpcTest, EmsCorrectsNoisyBlocks200_100)
{
// About 115 hard-decision bit errors per block at this Eb/N0.
EXPECT_EQ(count_ems_block_errors<BEIDOU_CNAV1_SUBFRAME2_SYMBOLS>(beidou_cnav1_ldpc_graph_200_100(), 3.5, 20, 2026U), 0);
}
TEST(BeidouCnav1LdpcTest, EmsCorrectsNoisyBlocks88_44)
{
EXPECT_EQ(count_ems_block_errors<BEIDOU_CNAV1_SUBFRAME3_SYMBOLS>(beidou_cnav1_ldpc_graph_88_44(), 3.5, 20, 2027U), 0);
}
@@ -15,6 +15,7 @@
*/
#include "Beidou_CNAV1.h"
#include "beidou_cnav1_navigation_message.h"
#include "beidou_cnav_test_helpers.h"
#include "rtklib_conversions.h"
#include <gtest/gtest.h>
#include <array>
@@ -69,8 +70,8 @@ void append_bch_codeword(
const auto encoded = encode_bch_lfsr_bipolar(msg_bits, n, k, feedback_pos_1based);
for (const auto cw : encoded)
{
// Map bipolar codeword to float symbols expected by the hard-decision BCH decoder.
frame.push_back(cw > 0 ? -1.0F : 1.0F);
// Bipolar +1 is bit 0; positive soft symbols favor bit 1. Same amplitude as SF2/SF3.
frame.push_back(cw > 0 ? -4.0F : 4.0F);
}
}
@@ -288,7 +289,7 @@ TEST(BeidouCnav1NavigationMessageTest, MultiFrameSteadyStateKeepsConstantSohOffs
TEST(BeidouCnav1NavigationMessageTest, RejectsEphemerisWhenIodeIodcMismatch)
{
// SF2 with IODC/IODE mismatch (ICD §7.4.3).
std::vector<uint8_t> sf2_bits(static_cast<size_t>(BEIDOU_CNAV1_SUBFRAME2_SYMBOLS), 0U);
std::array<uint8_t, BEIDOU_CNAV1_SF2_DATA_BITS> sf2_bits{0U};
const int32_t iodc_offset = 13 + 8;
const uint32_t iodc = 0x101U; // low 8 = 0x01
for (int32_t b = 0; b < 10; b++)
@@ -326,8 +327,8 @@ TEST(BeidouCnav1NavigationMessageTest, RejectsEphemerisWhenIodeIodcMismatch)
sf2_bits[static_cast<size_t>(BEIDOU_CNAV1_SF2_DATA_BITS - BEIDOU_CNAV1_CRC_BITS + b)] =
static_cast<uint8_t>((crc >> static_cast<uint32_t>(23 - b)) & 1U);
}
const auto sf2_llr = encode_bits_to_llr(sf2_bits);
const auto sf2_encoded = BeidouCnavTest::encode<BEIDOU_CNAV1_SUBFRAME2_SYMBOLS>(sf2_bits);
const auto sf2_llr = encode_bits_to_llr(std::vector<uint8_t>(sf2_encoded.begin(), sf2_encoded.end()));
const auto sf3_llr = encode_bits_to_llr(std::vector<uint8_t>(BEIDOU_CNAV1_SUBFRAME3_SYMBOLS, 0U));
std::vector<float> interleaved;
deinterleave_like_icd(sf2_llr, sf3_llr, interleaved);
@@ -341,3 +342,69 @@ TEST(BeidouCnav1NavigationMessageTest, RejectsEphemerisWhenIodeIodcMismatch)
ASSERT_TRUE(nav.decode_frame_symbols(frame.data(), BEIDOU_CNAV1_FRAME_SYMBOLS, 19));
EXPECT_FALSE(nav.have_new_ephemeris());
}
TEST(BeidouCnav1NavigationMessageTest, DecodesNoisyFrameRegardlessOfInputScale)
{
// Zero-payload frame (symbol amplitude 4) with AWGN on every symbol at Eb/N0 = 3.5 dB.
std::vector<float> frame = build_valid_zero_payload_frame(19U, 1U);
BeidouCnavTest::PortableGaussian gauss(1135U);
const double sigma = 4.0 * BeidouCnavTest::rate_half_bpsk_sigma(3.5);
int32_t hard_errors = 0;
for (int32_t i = 0; i < BEIDOU_CNAV1_FRAME_SYMBOLS; i++)
{
const float clean = frame[static_cast<size_t>(i)];
frame[static_cast<size_t>(i)] = static_cast<float>(clean + sigma * gauss());
if ((clean >= 0.0F) != (frame[static_cast<size_t>(i)] >= 0.0F))
{
hard_errors++;
}
}
ASSERT_GT(hard_errors, 50);
// Tracking correlator gain is arbitrary: decoding must not depend on it.
for (const float scale : {1.0e-3F, 1.0F, 1.0e3F})
{
std::vector<float> scaled(frame);
for (auto& value : scaled)
{
value *= scale;
}
Beidou_Cnav1_Navigation_Message nav;
EXPECT_TRUE(nav.decode_frame_symbols(scaled.data(), BEIDOU_CNAV1_FRAME_SYMBOLS, 19)) << "input scale " << scale;
}
}
TEST(BeidouCnav1NavigationMessageTest, SoftBchDecodingCorrectsBeyondHardDecisionRadius)
{
// Weak sign errors: 4 in the BCH(21,6) PRN word and 12 in the BCH(51,8) SOH
// word, one more than each code's hard-decision radius (t = 3 and t = 11).
const uint32_t soh = 123U;
std::vector<float> frame = build_valid_zero_payload_frame(19U, soh);
for (const int32_t i : {0, 5, 10, 15})
{
frame[static_cast<size_t>(i)] *= -0.1F;
}
for (int32_t i = 0; i < 12; i++)
{
frame[static_cast<size_t>(21 + 4 * i)] *= -0.1F;
}
Beidou_Cnav1_Navigation_Message nav;
EXPECT_TRUE(nav.probe_subframe1_prn(frame.data(), BEIDOU_CNAV1_FRAME_SYMBOLS, 19));
ASSERT_TRUE(nav.decode_frame_symbols(frame.data(), BEIDOU_CNAV1_FRAME_SYMBOLS, 19));
// Zero SF2 payload: HOW = 0, so TOW is the SOH alone.
EXPECT_DOUBLE_EQ(nav.get_tow_s(), static_cast<double>(soh * BEIDOU_CNAV1_SOH_LSB_S));
}
TEST(BeidouCnav1NavigationMessageTest, RejectsSohOutsideIcdRange)
{
// ICD Table 7-2: SOH is at most 3582 s, i.e. 199 in units of 18 s.
const std::vector<float> valid = build_valid_zero_payload_frame(19U, 199U);
const std::vector<float> invalid = build_valid_zero_payload_frame(19U, 200U);
Beidou_Cnav1_Navigation_Message nav;
EXPECT_TRUE(nav.decode_frame_symbols(valid.data(), BEIDOU_CNAV1_FRAME_SYMBOLS, 19));
int32_t fail_stage = -1;
EXPECT_FALSE(nav.decode_frame_symbols(invalid.data(), BEIDOU_CNAV1_FRAME_SYMBOLS, 19, &fail_stage));
EXPECT_EQ(fail_stage, 2);
EXPECT_FALSE(nav.probe_subframe1_prn(invalid.data(), BEIDOU_CNAV1_FRAME_SYMBOLS, 19));
}
@@ -13,8 +13,9 @@
*
* -----------------------------------------------------------------------------
*/
#include "Beidou_CNAV2.h"
#include "beidou_cnav2_ldpc.h"
#include "beidou_cnav2_test_helpers.h"
#include "beidou_cnav_test_helpers.h"
#include <gtest/gtest.h>
#include <algorithm>
#include <array>
@@ -65,7 +66,7 @@ TEST(BeidouCnav2LdpcTest, ReferenceEncoderMatchesIcdExample)
const auto codeword = cnav2_icd_bits();
std::array<uint8_t, 288> info{};
std::copy_n(codeword.begin(), info.size(), info.begin());
EXPECT_EQ(BeidouCnav2Test::encode(info), codeword);
EXPECT_EQ(BeidouCnavTest::encode<BEIDOU_CNAV2_LDPC_SYMBOLS>(info), codeword);
}
TEST(BeidouCnav2LdpcTest, IcdExampleSatisfiesEveryParityCheck)
@@ -85,7 +86,7 @@ TEST(BeidouCnav2LdpcTest, IcdExampleSatisfiesEveryParityCheck)
{
symbol = static_cast<uint8_t>((symbol << 1U) | bits[graph.check_to_var[e] * 6 + b]);
}
syndrome ^= BeidouCnav2Test::multiply(graph.check_to_h[e], symbol);
syndrome ^= BeidouCnavTest::multiply(graph.check_to_h[e], symbol);
}
EXPECT_EQ(syndrome, 0U) << "check " << row;
}
@@ -16,9 +16,11 @@
#include "Beidou_CNAV1.h"
#include "Beidou_CNAV2.h"
#include "beidou_cnav2_navigation_message.h"
#include "beidou_cnav2_test_helpers.h"
#include "beidou_cnav_test_helpers.h"
#include <gtest/gtest.h>
#include <algorithm>
#include <array>
#include <cmath>
#include <cstdint>
#include <limits>
#include <vector>
@@ -69,7 +71,7 @@ std::vector<float> make_frame(const std::array<uint8_t, BEIDOU_CNAV2_INFO_BITS>&
{
symbols[static_cast<size_t>(i)] = (preamble[i] == '1') ? -1.0F : 1.0F;
}
const auto codeword = BeidouCnav2Test::encode(info);
const auto codeword = BeidouCnavTest::encode<BEIDOU_CNAV2_LDPC_SYMBOLS>(info);
for (int i = 0; i < BEIDOU_CNAV2_LDPC_SYMBOLS; i++)
{
symbols[static_cast<size_t>(BEIDOU_CNAV2_PREAMBLE_SYMBOLS + i)] =
@@ -438,9 +440,15 @@ TEST(BeidouCnav2NavigationMessageTest, RejectsInvalidFrames)
Beidou_Cnav2_Navigation_Message nav;
EXPECT_FALSE(nav.decode_frame_symbols(nullptr, frame.size(), 19));
EXPECT_FALSE(nav.decode_frame_symbols(frame.data(), frame.size() - 1, 19));
// A preamble symbol error no longer rejects the frame: the soft preamble
// only resolves polarity, and LDPC plus CRC validate the frame.
frame[0] *= -1.0F;
EXPECT_FALSE(nav.decode_frame_symbols(frame.data(), frame.size(), 19));
EXPECT_TRUE(nav.decode_frame_symbols(frame.data(), frame.size(), 19));
frame[0] *= -1.0F;
// A preamble that cannot resolve the carrier polarity is rejected.
auto no_preamble = frame;
std::fill_n(no_preamble.begin(), BEIDOU_CNAV2_PREAMBLE_SYMBOLS, 0.0F);
EXPECT_FALSE(nav.decode_frame_symbols(no_preamble.data(), no_preamble.size(), 19));
frame.back() = std::numeric_limits<float>::quiet_NaN();
EXPECT_FALSE(nav.decode_frame_symbols(frame.data(), frame.size(), 19));
EXPECT_FALSE(nav.last_crc_ok());
@@ -479,3 +487,135 @@ TEST(BeidouCnav2NavigationMessageTest, SoftDecodingIsIndependentOfTrackingGain)
EXPECT_EQ(nav.last_sow(), 12345 * 3);
}
}
namespace
{
// Unit-amplitude 1 ms samples of one frame start: preamble spread by the data secondary code.
std::vector<float> preamble_1ms(float polarity)
{
const char* preamble = "111000100100110111101000";
const float secondary[5] = {1.0F, 1.0F, 1.0F, -1.0F, 1.0F};
std::vector<float> samples;
for (int i = 0; i < BEIDOU_CNAV2_PREAMBLE_SYMBOLS; i++)
{
const float symbol = (preamble[i] == '1') ? -1.0F : 1.0F;
for (const float chip : secondary)
{
samples.push_back(polarity * symbol * chip);
}
}
return samples;
}
} // namespace
TEST(BeidouCnav2NavigationMessageTest, PreambleStatisticIsScaleAndPolarityIndependent)
{
const double single = std::sqrt(static_cast<double>(BEIDOU_CNAV2_PREAMBLE_MS));
for (const float scale : {1.0e-3F, 1.0F, 1.0e3F})
{
auto positive = preamble_1ms(scale);
auto negative = preamble_1ms(-scale);
EXPECT_NEAR(Beidou_Cnav2_Navigation_Message::preamble_detection_statistic(positive.data()), single, 1e-6);
EXPECT_NEAR(Beidou_Cnav2_Navigation_Message::preamble_detection_statistic(negative.data()), single, 1e-6);
// Two frames with the same polarity combine coherently (+3 dB).
EXPECT_NEAR(Beidou_Cnav2_Navigation_Message::preamble_detection_statistic(positive.data(), positive.data()),
std::sqrt(2.0) * single, 1e-6);
// A polarity flip between frames falls back to the single-frame term.
EXPECT_NEAR(Beidou_Cnav2_Navigation_Message::preamble_detection_statistic(positive.data(), negative.data()), single, 1e-6);
}
const std::vector<float> zeros(static_cast<size_t>(BEIDOU_CNAV2_PREAMBLE_MS), 0.0F);
EXPECT_EQ(Beidou_Cnav2_Navigation_Message::preamble_detection_statistic(zeros.data()), 0.0);
}
TEST(BeidouCnav2NavigationMessageTest, PreambleStatisticDetectsWeakPreamblesWithRareFalseAlarms)
{
// 1 ms samples at Eb/N0 = 2 dB for the rate-1/2 LDPC (5 ms symbol Es/N0 = -1 dB),
// where a hard-decision correlation almost never reaches its threshold.
const double sigma = BeidouCnavTest::rate_half_bpsk_sigma(2.0) * std::sqrt(5.0);
BeidouCnavTest::PortableGaussian gauss(2026U);
const auto clean = preamble_1ms(1.0F);
const int trials = 200;
int detections = 0;
int false_alarms = 0;
std::vector<float> current(clean.size());
std::vector<float> previous(clean.size());
for (int t = 0; t < trials; t++)
{
for (size_t i = 0; i < clean.size(); i++)
{
current[i] = static_cast<float>(clean[i] + sigma * gauss());
previous[i] = static_cast<float>(clean[i] + sigma * gauss());
}
if (Beidou_Cnav2_Navigation_Message::preamble_detection_statistic(current.data(), previous.data()) >= BEIDOU_CNAV2_PREAMBLE_DETECTION_THRESHOLD)
{
detections++;
}
for (size_t i = 0; i < clean.size(); i++)
{
current[i] = static_cast<float>(sigma * gauss());
previous[i] = static_cast<float>(sigma * gauss());
}
if (Beidou_Cnav2_Navigation_Message::preamble_detection_statistic(current.data(), previous.data()) >= BEIDOU_CNAV2_PREAMBLE_DETECTION_THRESHOLD)
{
false_alarms++;
}
}
EXPECT_GE(detections, 190);
EXPECT_EQ(false_alarms, 0);
}
TEST(BeidouCnav2NavigationMessageTest, DecodesFrameWithPreambleSymbolErrors)
{
auto frame = make_frame(make_info(19, BEIDOU_CNAV2_MSG_EPH1, 1000));
// Weak sign errors on 5 of the 24 preamble symbols.
for (const int i : {0, 5, 10, 15, 20})
{
frame[static_cast<size_t>(i)] *= -0.2F;
}
Beidou_Cnav2_Navigation_Message nav;
EXPECT_TRUE(nav.decode_frame_symbols(frame.data(), BEIDOU_CNAV2_FRAME_SYMBOLS, 19));
// Inverted carrier polarity is resolved from the soft preamble correlation.
for (auto& value : frame)
{
value = -value;
}
EXPECT_TRUE(nav.decode_frame_symbols(frame.data(), BEIDOU_CNAV2_FRAME_SYMBOLS, 19));
}
TEST(BeidouCnav2NavigationMessageTest, DecodesNoisyFrameNearLdpcThreshold)
{
// Symbol amplitude 4 everywhere, AWGN on every symbol at Eb/N0 = 3 dB.
auto frame = make_frame(make_info(19, BEIDOU_CNAV2_MSG_EPH1, 1000));
for (int i = 0; i < BEIDOU_CNAV2_PREAMBLE_SYMBOLS; i++)
{
frame[static_cast<size_t>(i)] *= 4.0F;
}
BeidouCnavTest::PortableGaussian gauss(1135U);
const double sigma = 4.0 * BeidouCnavTest::rate_half_bpsk_sigma(3.0);
int preamble_errors = 0;
for (int i = 0; i < BEIDOU_CNAV2_FRAME_SYMBOLS; i++)
{
const float clean = frame[static_cast<size_t>(i)];
frame[static_cast<size_t>(i)] = static_cast<float>(clean + sigma * gauss());
if (i < BEIDOU_CNAV2_PREAMBLE_SYMBOLS && ((clean >= 0.0F) != (frame[static_cast<size_t>(i)] >= 0.0F)))
{
preamble_errors++;
}
}
ASSERT_GT(preamble_errors, 0);
for (const float scale : {1.0e-3F, 1.0F, 1.0e3F})
{
std::vector<float> scaled(frame);
for (auto& value : scaled)
{
value *= scale;
}
Beidou_Cnav2_Navigation_Message nav;
EXPECT_TRUE(nav.decode_frame_symbols(scaled.data(), BEIDOU_CNAV2_FRAME_SYMBOLS, 19)) << "input scale " << scale;
}
}
@@ -1,120 +0,0 @@
/*!
* \file beidou_cnav2_test_helpers.h
* \brief Independent B-CNAV2 reference encoder for tests
* \author Carles Fernandez-Prades, 2026. cfernandez(at)cttc.es
*
* -----------------------------------------------------------------------------
*
* GNSS-SDR is a Global Navigation Satellite System software-defined receiver.
* This file is part of GNSS-SDR.
*
* Copyright (C) 2010-2026 (see AUTHORS file for a list of contributors)
* SPDX-License-Identifier: GPL-3.0-or-later
*
* -----------------------------------------------------------------------------
*/
#ifndef GNSS_SDR_BEIDOU_CNAV2_TEST_HELPERS_H
#define GNSS_SDR_BEIDOU_CNAV2_TEST_HELPERS_H
#include "beidou_cnav2_ldpc.h"
#include <array>
#include <stdexcept>
#include <utility>
namespace BeidouCnav2Test
{
// Polynomial arithmetic independent of the decoder's logarithm lookup tables.
inline uint8_t multiply(uint8_t a, uint8_t b)
{
uint8_t product = 0;
while (b != 0)
{
if ((b & 1U) != 0)
{
product ^= a;
}
b >>= 1U;
a <<= 1U;
if ((a & 0x40U) != 0)
{
a ^= 0x43U;
}
}
return product;
}
// Solve H2 * parity = H1 * information by Gaussian elimination over GF(64).
inline std::array<uint8_t, 576> encode(const std::array<uint8_t, 288>& bits)
{
std::array<uint8_t, 96> symbols{};
for (size_t i = 0; i < bits.size(); i++)
{
symbols[i / 6] = static_cast<uint8_t>((symbols[i / 6] << 1U) | bits[i]);
}
std::array<std::array<uint8_t, 49>, 48> augmented{};
for (size_t row = 0; row < 48; row++)
{
for (size_t edge = 0; edge < 4; edge++)
{
// Each printed row contains four groups of four entries.
const size_t offset = (row % 12) * 16 + (row / 12) * 4 + edge;
const auto col = BEIDOU_CNAV2_H48_96_INDEX[offset];
const auto h = BEIDOU_CNAV2_H48_96_ELEMENT[offset];
if (col < 48)
{
augmented[row][48] ^= multiply(h, symbols[col]);
}
else
{
augmented[row][col - 48] = h;
}
}
}
for (size_t col = 0; col < 48; col++)
{
size_t pivot = col;
while (pivot < 48 && augmented[pivot][col] == 0)
{
pivot++;
}
if (pivot == 48)
{
throw std::runtime_error("Singular B-CNAV2 parity matrix");
}
std::swap(augmented[col], augmented[pivot]);
uint8_t inverse = 1;
while (multiply(augmented[col][col], inverse) != 1)
{
inverse++;
}
for (size_t j = col; j <= 48; j++)
{
augmented[col][j] = multiply(augmented[col][j], inverse);
}
for (size_t row = 0; row < 48; row++)
{
if (row == col)
{
continue;
}
const auto factor = augmented[row][col];
for (size_t j = col; j <= 48; j++)
{
augmented[row][j] ^= multiply(factor, augmented[col][j]);
}
}
}
for (size_t row = 0; row < 48; row++)
{
symbols[48 + row] = augmented[row][48];
}
std::array<uint8_t, 576> codeword{};
for (size_t bit = 0; bit < codeword.size(); bit++)
{
codeword[bit] = (symbols[bit / 6] >> (5 - bit % 6)) & 1U;
}
return codeword;
}
} // namespace BeidouCnav2Test
#endif // GNSS_SDR_BEIDOU_CNAV2_TEST_HELPERS_H
@@ -0,0 +1,169 @@
/*!
* \file beidou_cnav_test_helpers.h
* \brief Independent B-CNAV1/B-CNAV2 reference encoder for tests
* \author Carles Fernandez-Prades, 2026. cfernandez(at)cttc.es
*
* -----------------------------------------------------------------------------
*
* GNSS-SDR is a Global Navigation Satellite System software-defined receiver.
* This file is part of GNSS-SDR.
*
* Copyright (C) 2010-2026 (see AUTHORS file for a list of contributors)
* SPDX-License-Identifier: GPL-3.0-or-later
*
* -----------------------------------------------------------------------------
*/
#ifndef GNSS_SDR_BEIDOU_CNAV_TEST_HELPERS_H
#define GNSS_SDR_BEIDOU_CNAV_TEST_HELPERS_H
#include "beidou_cnav1_ldpc.h"
#include "beidou_cnav2_ldpc.h"
#include <array>
#include <cmath>
#include <cstdint>
#include <random>
#include <stdexcept>
#include <utility>
namespace BeidouCnavTest
{
// Polynomial arithmetic independent of the decoder's logarithm lookup tables.
inline uint8_t multiply(uint8_t a, uint8_t b)
{
uint8_t product = 0;
while (b != 0)
{
if ((b & 1U) != 0)
{
product ^= a;
}
b >>= 1U;
a <<= 1U;
if ((a & 0x40U) != 0)
{
a ^= 0x43U;
}
}
return product;
}
// Solve H2 * parity = H1 * information by Gaussian elimination over GF(64).
template <int N>
inline std::array<uint8_t, N> encode(const std::array<uint8_t, N / 2>& bits)
{
if (N != 1200 && N != 528 && N != 576)
{
throw std::runtime_error("Invalid B-CNAV1/B-CNAV2 encoded subframe size");
}
std::array<uint8_t, N / 6> symbols{};
for (size_t i = 0; i < bits.size(); i++)
{
symbols[i / 6] = static_cast<uint8_t>((symbols[i / 6] << 1U) | bits[i]);
}
std::array<std::array<uint8_t, N / 12 + 1>, N / 12> augmented{};
for (size_t row = 0; row < N / 12; row++)
{
for (size_t edge = 0; edge < 4; edge++)
{
// Each printed row contains four groups of four entries.
const size_t offset = (row % (N / 48)) * 16 + (row / (N / 48)) * 4 + edge;
uint16_t col = 0;
uint8_t h = 0;
switch (N)
{
case 1200:
col = BEIDOU_CNAV1_H100_200_INDEX[offset];
h = BEIDOU_CNAV1_H100_200_ELEMENT[offset];
break;
case 528:
col = BEIDOU_CNAV1_H44_88_INDEX[offset];
h = BEIDOU_CNAV1_H44_88_ELEMENT[offset];
break;
case 576:
col = BEIDOU_CNAV2_H48_96_INDEX[offset];
h = BEIDOU_CNAV2_H48_96_ELEMENT[offset];
break;
}
if (col < N / 12)
{
augmented[row][N / 12] ^= multiply(h, symbols[col]);
}
else
{
augmented[row][col - N / 12] = h;
}
}
}
for (size_t col = 0; col < N / 12; col++)
{
size_t pivot = col;
while (pivot < N / 12 && augmented[pivot][col] == 0)
{
pivot++;
}
if (pivot == N / 12)
{
throw std::runtime_error("Singular B-CNAV1/B-CNAV2 parity matrix");
}
std::swap(augmented[col], augmented[pivot]);
uint8_t inverse = 1;
while (multiply(augmented[col][col], inverse) != 1)
{
inverse++;
}
for (size_t j = col; j <= N / 12; j++)
{
augmented[col][j] = multiply(augmented[col][j], inverse);
}
for (size_t row = 0; row < N / 12; row++)
{
if (row == col)
{
continue;
}
const auto factor = augmented[row][col];
for (size_t j = col; j <= N / 12; j++)
{
augmented[row][j] ^= multiply(factor, augmented[col][j]);
}
}
}
for (size_t row = 0; row < N / 12; row++)
{
symbols[N / 12 + row] = augmented[row][N / 12];
}
std::array<uint8_t, N> codeword{};
for (size_t bit = 0; bit < codeword.size(); bit++)
{
codeword[bit] = (symbols[bit / 6] >> (5 - bit % 6)) & 1U;
}
return codeword;
}
// Standard normal samples that are reproducible across standard libraries
// (std::normal_distribution is implementation-defined; std::mt19937 is not).
class PortableGaussian
{
public:
explicit PortableGaussian(uint32_t seed) : rng_(seed) {}
double operator()()
{
const double u1 = (static_cast<double>(rng_()) + 0.5) / 4294967296.0;
const double u2 = (static_cast<double>(rng_()) + 0.5) / 4294967296.0;
return std::sqrt(-2.0 * std::log(u1)) * std::cos(6.283185307179586 * u2);
}
private:
std::mt19937 rng_;
};
// Noise standard deviation for unit-amplitude BPSK carrying a rate-1/2 code.
inline double rate_half_bpsk_sigma(double ebn0_db)
{
const double esn0 = std::pow(10.0, (ebn0_db - 10.0 * std::log10(2.0)) / 10.0);
return std::sqrt(1.0 / (2.0 * esn0));
}
} // namespace BeidouCnavTest
#endif // GNSS_SDR_BEIDOU_CNAV_TEST_HELPERS_H