Fix Acquisition dump readers

Update changelog Add Joel to the CITATION.cff file
2026-10-09 17:21:51 +00:00 · 2026-09-01 09:22:38 +02:00
parent b980b55ed2
commit 0af508713a
6 changed files with 242 additions and 159 deletions
@@ -20,9 +20,30 @@
 #include <matio.h>
 #include <algorithm>
 #include <cmath>
+#include <cstdint>
 #include <iostream>
 #include <utility>

+namespace
+{
+// Reads a 1x1 variable from the dump, writing fallback instead when the
+// variable is absent (e.g. num_dwells in pcps_acquisition_fine_doppler_cc
+// dumps).
+template <typename T>
+void read_scalar_or(mat_t* matfile, const char* name, T& value, T fallback)
+{
+    matvar_t* var = Mat_VarRead(matfile, name);
+    if (var == nullptr)
+        {
+            value = fallback;
+            return;
+        }
+    value = *static_cast<T*>(var->data);
+    Mat_VarFree(var);
+}
+}  // namespace
+
+
 bool Acquisition_Dump_Reader::read_binary_acq()
 {
    mat_t* matfile = Mat_Open(d_dump_filename.c_str(), MAT_ACC_RDONLY);
@@ -45,6 +66,72 @@ bool Acquisition_Dump_Reader::read_binary_acq()
            Mat_Close(matfile);
            return false;
        }
+    if (var_->data_type != MAT_T_SINGLE)
+        {
+            std::cout << "Invalid Acquisition dump file: data type error\n";
+            Mat_VarFree(var_);
+            Mat_Close(matfile);
+            return false;
+        }
+
+    // Read the Doppler metadata before validating the grid dimensions: dumps
+    // written by an assisted (narrowed) acquisition carry two grid columns (the
+    // assisted candidate and its noise reference) encoded as doppler_max = 0,
+    // doppler_step = <configured doppler_max>, so the expected column count
+    // follows from the file's own metadata rather than from the full-grid
+    // parameters this reader was constructed with.
+    matvar_t* var2_ = Mat_VarRead(matfile, "doppler_max");
+    if (var2_ == nullptr)
+        {
+            std::cout << "Unreachable doppler_max variable in Acquisition dump file.\n";
+            Mat_VarFree(var_);
+            Mat_Close(matfile);
+            return false;
+        }
+    d_doppler_max = *static_cast<unsigned int*>(var2_->data);
+    Mat_VarFree(var2_);
+
+    var2_ = Mat_VarRead(matfile, "doppler_step");
+    if (var2_ == nullptr)
+        {
+            std::cout << "Unreachable doppler_step variable in Acquisition dump file.\n";
+            Mat_VarFree(var_);
+            Mat_Close(matfile);
+            return false;
+        }
+    d_doppler_step = *static_cast<unsigned int*>(var2_->data);
+    Mat_VarFree(var2_);
+
+    // Optional variables, absent from dumps written by older versions and by
+    // blocks that never narrow (e.g. pcps_acquisition_fine_doppler_cc): their
+    // absence means a full grid centered at 0 Hz.
+    read_scalar_or(matfile, "doppler_center", doppler_center, 0);
+    int32_t narrowed_flag = 0;
+    read_scalar_or(matfile, "doppler_narrowed", narrowed_flag, static_cast<int32_t>(0));
+    doppler_narrowed = (narrowed_flag != 0);
+
+    if (d_doppler_step == 0)
+        {
+            d_doppler_step = 1;
+        }
+    if (doppler_narrowed)
+        {
+            d_num_doppler_bins = 2U;
+        }
+    else
+        {
+            // pcps_acquisition sizes its grid with ceil(2*doppler_max/doppler_step),
+            // but pcps_acquisition_fine_doppler_cc sizes it with floor(): accept
+            // either count when they differ.
+            const double bins_exact = static_cast<double>(2 * d_doppler_max) / static_cast<double>(d_doppler_step);
+            d_num_doppler_bins = static_cast<unsigned int>(std::ceil(bins_exact));
+            const auto bins_floor = static_cast<unsigned int>(std::floor(bins_exact));
+            if ((var_->dims[1] == bins_floor) && (bins_floor != 0U))
+                {
+                    d_num_doppler_bins = bins_floor;
+                }
+        }
+
    if ((var_->dims[0] != d_samples_per_code) || (var_->dims[1] != d_num_doppler_bins))
        {
            std::cout << "Invalid Acquisition dump file: dimension matrix error\n";
@@ -60,44 +147,25 @@ bool Acquisition_Dump_Reader::read_binary_acq()
            Mat_Close(matfile);
            return false;
        }
-    if (var_->data_type != MAT_T_SINGLE)
+
+    // Rebuild the Doppler axis and the magnitude storage from the file's
+    // metadata. doppler(i) = -doppler_max + doppler_center + doppler_step * i
+    // is the same decoding pcps_acquisition::compute_statistics() uses, and it
+    // maps a narrowed dump's two columns to {doppler_center, doppler_center +
+    // configured doppler_max}.
+    doppler.clear();
+    for (unsigned int doppler_index = 0; doppler_index < d_num_doppler_bins; doppler_index++)
        {
-            std::cout << "Invalid Acquisition dump file: data type error\n";
-            Mat_VarFree(var_);
-            Mat_Close(matfile);
-            return false;
+            doppler.push_back(-static_cast<int>(d_doppler_max) + doppler_center + static_cast<int>(d_doppler_step) * static_cast<int>(doppler_index));
        }
-    matvar_t* var2_ = Mat_VarRead(matfile, "doppler_max");
-    d_doppler_max = *static_cast<unsigned int*>(var2_->data);
-    Mat_VarFree(var2_);
+    mag.assign(d_num_doppler_bins, std::vector<float>(d_samples_per_code));

-    var2_ = Mat_VarRead(matfile, "doppler_step");
-    d_doppler_step = *static_cast<unsigned int*>(var2_->data);
-    Mat_VarFree(var2_);
-
-    var2_ = Mat_VarRead(matfile, "input_power");
-    input_power = *static_cast<float*>(var2_->data);
-    Mat_VarFree(var2_);
-
-    var2_ = Mat_VarRead(matfile, "acq_doppler_hz");
-    acq_doppler_hz = *static_cast<float*>(var2_->data);
-    Mat_VarFree(var2_);
-
-    var2_ = Mat_VarRead(matfile, "acq_delay_samples");
-    acq_delay_samples = *static_cast<float*>(var2_->data);
-    Mat_VarFree(var2_);
-
-    var2_ = Mat_VarRead(matfile, "test_statistic");
-    test_statistic = *static_cast<float*>(var2_->data);
-    Mat_VarFree(var2_);
-
-    var2_ = Mat_VarRead(matfile, "threshold");
-    threshold = *static_cast<float*>(var2_->data);
-    Mat_VarFree(var2_);
-
-    var2_ = Mat_VarRead(matfile, "sample_counter");
-    sample_counter = *static_cast<uint64_t*>(var2_->data);
-    Mat_VarFree(var2_);
+    read_scalar_or(matfile, "input_power", input_power, 0.0F);
+    read_scalar_or(matfile, "acq_doppler_hz", acq_doppler_hz, 0.0F);
+    read_scalar_or(matfile, "acq_delay_samples", acq_delay_samples, 0.0F);
+    read_scalar_or(matfile, "test_statistic", test_statistic, 0.0F);
+    read_scalar_or(matfile, "threshold", threshold, 0.0F);
+    read_scalar_or(matfile, "sample_counter", sample_counter, static_cast<uint64_t>(0));

    var2_ = Mat_VarRead(matfile, "positive_acq");
    if (var2_ == nullptr)
@@ -114,19 +182,20 @@ bool Acquisition_Dump_Reader::read_binary_acq()
    positive_acq = *static_cast<int*>(var2_->data);
    Mat_VarFree(var2_);

-    var2_ = Mat_VarRead(matfile, "num_dwells");
-    num_dwells = *static_cast<int*>(var2_->data);
-    Mat_VarFree(var2_);
-
-    var2_ = Mat_VarRead(matfile, "PRN");
-    PRN = *static_cast<int*>(var2_->data);
-    Mat_VarFree(var2_);
+    read_scalar_or(matfile, "num_dwells", num_dwells, 0U);
+    read_scalar_or(matfile, "PRN", PRN, 0U);

    std::vector<std::vector<float> >::iterator it1;
    std::vector<float>::iterator it2;
    auto* aux = static_cast<float*>(var_->data);
    int k = 0;
    float normalization_factor = std::pow(d_samples_per_code, 4) * input_power;
+    if (!(normalization_factor > 0.0F))
+        {
+            // Non-CFAR pcps dumps and fine-doppler dumps store input_power = 0:
+            // keep the raw grid values instead of dividing by zero.
+            normalization_factor = 1.0F;
+        }
    for (it1 = mag.begin(); it1 != mag.end(); it1++)
        {
            for (it2 = it1->begin(); it2 != it1->end(); it2++)
@@ -151,30 +220,47 @@ Acquisition_Dump_Reader::Acquisition_Dump_Reader(const std::string& basename,
    unsigned int doppler_step_ = 0;
    unsigned int samples_per_code_ = 0;

-    mat_t* matfile = Mat_Open(d_dump_filename.c_str(), MAT_ACC_RDONLY);
+    // The dump filename embeds the PRN, which is exactly what this constructor
+    // has to discover, so probe the candidate PRNs for an existing file
+    // (1..210 covers every supported system, QZSS included).
+    mat_t* matfile = nullptr;
+    for (unsigned int candidate_sat = 1; (candidate_sat <= 210) && (matfile == nullptr); candidate_sat++)
+        {
+            const std::string candidate_filename = basename + "_ch_" + std::to_string(channel) + "_" + std::to_string(execution) + "_sat_" + std::to_string(candidate_sat) + ".mat";
+            matfile = Mat_Open(candidate_filename.c_str(), MAT_ACC_RDONLY);
+            if (matfile != nullptr)
+                {
+                    sat_ = candidate_sat;
+                }
+        }
    if (matfile != nullptr)
        {
            matvar_t* var_ = Mat_VarRead(matfile, "doppler_max");
-            doppler_max_ = *static_cast<unsigned int*>(var_->data);
-            Mat_VarFree(var_);
+            if (var_ != nullptr)
+                {
+                    doppler_max_ = *static_cast<unsigned int*>(var_->data);
+                    Mat_VarFree(var_);
+                }

            var_ = Mat_VarRead(matfile, "doppler_step");
-            doppler_step_ = *static_cast<unsigned int*>(var_->data);
-            Mat_VarFree(var_);
+            if (var_ != nullptr)
+                {
+                    doppler_step_ = *static_cast<unsigned int*>(var_->data);
+                    Mat_VarFree(var_);
+                }

-            var_ = Mat_VarRead(matfile, "PRN");
-            sat_ = *static_cast<int*>(var_->data);
-            Mat_VarFree(var_);
-
-            var_ = Mat_VarRead(matfile, "grid");
-            samples_per_code_ = var_->dims[0];
-            Mat_VarFree(var_);
+            var_ = Mat_VarRead(matfile, "acq_grid");
+            if (var_ != nullptr)
+                {
+                    samples_per_code_ = var_->dims[0];
+                    Mat_VarFree(var_);
+                }

            Mat_Close(matfile);
        }
    else
        {
-            std::cout << "Unreachable Acquisition dump file " << d_dump_filename << '\n';
+            std::cout << "Unreachable Acquisition dump file " << basename << "_ch_" << channel << "_" << execution << "_sat_<PRN>.mat\n";
        }
    acq_doppler_hz = 0.0;
    acq_delay_samples = 0.0;
@@ -231,87 +317,3 @@ Acquisition_Dump_Reader::Acquisition_Dump_Reader(const std::string& basename,
            samples.push_back(k);
        }
 }
-
-
-// Copy assignment operator
-Acquisition_Dump_Reader& Acquisition_Dump_Reader::operator=(const Acquisition_Dump_Reader& other)
-{
-    if (this != &other)
-        {
-            doppler = other.doppler;
-            samples = other.samples;
-            mag = other.mag;
-            acq_doppler_hz = other.acq_doppler_hz;
-            acq_delay_samples = other.acq_delay_samples;
-            test_statistic = other.test_statistic;
-            input_power = other.input_power;
-            threshold = other.threshold;
-            positive_acq = other.positive_acq;
-            PRN = other.PRN;
-            num_dwells = other.num_dwells;
-            sample_counter = other.sample_counter;
-            d_basename = other.d_basename;
-            d_dump_filename = other.d_dump_filename;
-            d_sat = other.d_sat;
-            d_doppler_max = other.d_doppler_max;
-            d_doppler_step = other.d_doppler_step;
-            d_samples_per_code = other.d_samples_per_code;
-            d_num_doppler_bins = other.d_num_doppler_bins;
-        }
-    return *this;
-}
-
-
-// Move constructor
-Acquisition_Dump_Reader::Acquisition_Dump_Reader(Acquisition_Dump_Reader&& other) noexcept
-    : doppler(std::move(other.doppler)),
-      samples(std::move(other.samples)),
-      mag(std::move(other.mag)),
-      acq_doppler_hz(other.acq_doppler_hz),
-      acq_delay_samples(other.acq_delay_samples),
-      test_statistic(other.test_statistic),
-      input_power(other.input_power),
-      threshold(other.threshold),
-      positive_acq(other.positive_acq),
-      PRN(other.PRN),
-      num_dwells(other.num_dwells),
-      sample_counter(other.sample_counter),
-      d_basename(std::move(other.d_basename)),
-      d_dump_filename(std::move(other.d_dump_filename)),
-      d_sat(other.d_sat),
-      d_doppler_max(other.d_doppler_max),
-      d_doppler_step(other.d_doppler_step),
-      d_samples_per_code(other.d_samples_per_code),
-      d_num_doppler_bins(other.d_num_doppler_bins)
-{
-}
-
-
-// Move assignment operator
-Acquisition_Dump_Reader& Acquisition_Dump_Reader::operator=(Acquisition_Dump_Reader&& other) noexcept
-{
-    if (this != &other)  // Check for self-assignment
-        {
-            // Move member variables from the other object to this object
-            d_basename = std::move(other.d_basename);
-            d_dump_filename = std::move(other.d_dump_filename);
-            d_sat = other.d_sat;
-            d_doppler_max = other.d_doppler_max;
-            d_doppler_step = other.d_doppler_step;
-            d_samples_per_code = other.d_samples_per_code;
-            d_num_doppler_bins = other.d_num_doppler_bins;
-            doppler = std::move(other.doppler);
-            samples = std::move(other.samples);
-            mag = std::move(other.mag);
-            acq_doppler_hz = other.acq_doppler_hz;
-            acq_delay_samples = other.acq_delay_samples;
-            test_statistic = other.test_statistic;
-            input_power = other.input_power;
-            threshold = other.threshold;
-            positive_acq = other.positive_acq;
-            PRN = other.PRN;
-            num_dwells = other.num_dwells;
-            sample_counter = other.sample_counter;
-        }
-    return *this;
-}
@@ -37,10 +37,10 @@ public:
        int channel = 0,
        int execution = 1);

-    Acquisition_Dump_Reader(const Acquisition_Dump_Reader& other) = default;       //!< Copy constructor
-    Acquisition_Dump_Reader& operator=(const Acquisition_Dump_Reader& other);      //!< Copy assignment operator
-    Acquisition_Dump_Reader(Acquisition_Dump_Reader&& other) noexcept;             //!< Move constructor
-    Acquisition_Dump_Reader& operator=(Acquisition_Dump_Reader&& other) noexcept;  //!< Move assignment operator
+    Acquisition_Dump_Reader(const Acquisition_Dump_Reader& other) = default;                 //!< Copy constructor
+    Acquisition_Dump_Reader& operator=(const Acquisition_Dump_Reader& other) = default;      //!< Copy assignment operator
+    Acquisition_Dump_Reader(Acquisition_Dump_Reader&& other) noexcept = default;             //!< Move constructor
+    Acquisition_Dump_Reader& operator=(Acquisition_Dump_Reader&& other) noexcept = default;  //!< Move assignment operator

    bool read_binary_acq();

@@ -53,6 +53,8 @@ public:
    float input_power{};
    float threshold{};
    int positive_acq{};
+    int doppler_center{};     //!< Doppler grid center (Hz); 0 for dumps that predate this variable
+    bool doppler_narrowed{};  //!< true if the dump was written by an assisted (narrowed) acquisition: 2 columns, {center, center + doppler_max}
    unsigned int PRN{};
    unsigned int num_dwells{};
    uint64_t sample_counter{};