From fd274fb8919f9f6c00444aa880c40beb04fafe60 Mon Sep 17 00:00:00 2001 From: Carles Fernandez Date: Wed, 23 Sep 2026 12:46:37 +0200 Subject: [PATCH] Make acquisition dump readers skip noise-reference columns With set_doppler_num_bins(), acq_grid holds the Doppler candidate columns followed by 0, 1 or 2 noise-reference columns placed at doppler_center +/- the configured doppler_max, outside the linear Doppler axis. Nothing in the dump said how many columns were candidates, so the readers either rejected the file (Acquisition_Dump_Reader) or plotted the reference columns on a wrong Doppler axis (plot_acq_grid.py, plot_acq_grid.m). Write a new doppler_num_candidates variable to the dump and use it in the readers to keep only the candidate columns. Dumps without it are read as before, except that the reference column of older narrowed dumps is also skipped. The Python and MATLAB scripts plot a single-candidate (narrowed) grid as a code-delay line instead of an empty surface. Add a regression test that reads back a full-grid dump with two reference columns and a narrowed N = 5 dump. --- .../gnuradio_blocks/pcps_acquisition.cc | 5 ++ ...pcps_acquisition_doppler_narrowing_test.cc | 61 ++++++++++++++ .../libs/acquisition_dump_reader.cc | 42 ++++++---- .../libs/acquisition_dump_reader.h | 2 +- utils/matlab/plot_acq_grid.m | 38 +++++++-- utils/python/plot_acq_grid.py | 81 ++++++++++++++----- 6 files changed, 189 insertions(+), 40 deletions(-) diff --git a/src/algorithms/acquisition/gnuradio_blocks/pcps_acquisition.cc b/src/algorithms/acquisition/gnuradio_blocks/pcps_acquisition.cc index c17fff78d..8dfdc6179 100644 --- a/src/algorithms/acquisition/gnuradio_blocks/pcps_acquisition.cc +++ b/src/algorithms/acquisition/gnuradio_blocks/pcps_acquisition.cc @@ -582,6 +582,11 @@ void pcps_acquisition::dump_results(const AcquisitionResult& result) write_matlab_var<1>("doppler_step", dump_doppler_step, matfp, dims_1d); write_matlab_var<1>("doppler_center", d_doppler_center, matfp, dims_1d); write_matlab_var<1>("doppler_narrowed", static_cast(dump_narrowed ? 1 : 0), matfp, dims_1d); + // Number of leading acq_grid columns that are real Doppler candidates. + // Any remaining trailing columns (0, 1 or 2) are noise-reference rows + // at doppler_center +/- the configured doppler_max, outside the linear + // Doppler axis above, and must be skipped by readers. + write_matlab_var<1>("doppler_num_candidates", static_cast(dump_candidate_count), matfp, dims_1d); write_matlab_var<1>("positive_acq", static_cast(result.positive_acq ? 1 : 0), matfp, dims_1d); write_matlab_var<1>("acq_doppler_hz", static_cast(d_gnss_synchro->Acq_doppler_hz), matfp, dims_1d); write_matlab_var<1>("acq_delay_samples", static_cast(d_gnss_synchro->Acq_delay_samples), matfp, dims_1d); diff --git a/tests/unit-tests/signal-processing-blocks/acquisition/pcps_acquisition_doppler_narrowing_test.cc b/tests/unit-tests/signal-processing-blocks/acquisition/pcps_acquisition_doppler_narrowing_test.cc index 516143aa2..80536b0e2 100644 --- a/tests/unit-tests/signal-processing-blocks/acquisition/pcps_acquisition_doppler_narrowing_test.cc +++ b/tests/unit-tests/signal-processing-blocks/acquisition/pcps_acquisition_doppler_narrowing_test.cc @@ -18,6 +18,7 @@ #include "GPS_L1_CA.h" +#include "acquisition_dump_reader.h" #include "concurrent_queue.h" #include "gnss_block_interface.h" #include "gnss_sdr_filesystem.h" @@ -29,6 +30,7 @@ #include #include #include +#include #include #include #include @@ -705,3 +707,62 @@ TEST_F(PcpsAcquisitionDopplerNarrowingTest /*unused*/, NumBinsJustBelowFullGridS ASSERT_EQ(1, result.rx_message) << "Acquisition failure with N = d_num_doppler_bins - 1."; EXPECT_LE(result.doppler_error_hz, 1000); } + + +TEST_F(PcpsAcquisitionDopplerNarrowingTest /*unused*/, DumpReaderSkipsReferenceColumns /*unused*/) +{ + // Dumps carry doppler_num_candidates so readers can tell the candidate + // columns from the trailing noise-reference columns. Covers a full grid + // with two reference columns (CFAR, doppler_max = 3200, step = 125: the + // wraparound distance is too small, 52 candidates + 2 references) and a + // narrowed N = 5 search (5 candidates + 2 references). + struct Case + { + unsigned int doppler_max; + unsigned int doppler_step; + int center; + unsigned int num_bins; + unsigned int expected_candidates; + }; + const std::array cases{{{3200, 125, 0, 0, 52}, {5000, 200, static_cast(kTrueDopplerHz) - 400, 5, 5}}}; + for (const auto &c : cases) + { + config = std::make_shared(); + init(c.doppler_max, c.doppler_step, /*use_cfar=*/true); + const std::string data_str = "./tmp-acq-reader-dump"; + if (fs::exists(data_str)) + { + fs::remove_all(data_str); + } + fs::create_directory(data_str); + config->supersede_property("Acquisition_1C.dump", "true"); + config->supersede_property("Acquisition_1C.dump_filename", data_str + "/acquisition"); + config->supersede_property("Acquisition_1C.dump_channel", "1"); + + const RunResult result = run_acquisition(top_block, config.get(), gnss_synchro, c.center, c.num_bins); + ASSERT_EQ(1, result.rx_message) << "Acquisition failure while producing the dump."; + + Acquisition_Dump_Reader reader(data_str + "/acquisition_G_1C", 1, 1); + ASSERT_TRUE(reader.read_binary_acq()) << "Dump not readable by Acquisition_Dump_Reader"; + ASSERT_EQ(c.expected_candidates, reader.doppler.size()); + ASSERT_EQ(c.expected_candidates, reader.mag.size()); + const int half_span = (static_cast(c.expected_candidates) - 1) / 2; + EXPECT_EQ(c.center - half_span * static_cast(c.doppler_step), reader.doppler.front()); + EXPECT_EQ(static_cast(c.doppler_step), reader.doppler[1] - reader.doppler[0]); + + // The grid maximum must sit at the Doppler the acquisition reported. + size_t best_bin = 0; + float best_value = -1.0F; + for (size_t k = 0; k < reader.mag.size(); k++) + { + const float value = *std::max_element(reader.mag[k].begin(), reader.mag[k].end()); + if (value > best_value) + { + best_value = value; + best_bin = k; + } + } + EXPECT_EQ(static_cast(reader.acq_doppler_hz), reader.doppler[best_bin]); + fs::remove_all(data_str); + } +} diff --git a/tests/unit-tests/signal-processing-blocks/libs/acquisition_dump_reader.cc b/tests/unit-tests/signal-processing-blocks/libs/acquisition_dump_reader.cc index 462d4a14b..b674d1425 100644 --- a/tests/unit-tests/signal-processing-blocks/libs/acquisition_dump_reader.cc +++ b/tests/unit-tests/signal-processing-blocks/libs/acquisition_dump_reader.cc @@ -74,12 +74,9 @@ bool Acquisition_Dump_Reader::read_binary_acq() 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 = , so the expected column count - // follows from the file's own metadata rather than from the full-grid - // parameters this reader was constructed with. + // Read the Doppler metadata before validating the grid dimensions: 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) { @@ -109,14 +106,28 @@ bool Acquisition_Dump_Reader::read_binary_acq() int32_t narrowed_flag = 0; read_scalar_or(matfile, "doppler_narrowed", narrowed_flag, static_cast(0)); doppler_narrowed = (narrowed_flag != 0); + // Number of leading acq_grid columns that are Doppler candidates. Up to two + // trailing columns are noise-reference bins outside the Doppler axis, which + // are skipped. Absent from older dumps: there, a narrowed dump holds one + // candidate plus one noise-reference column, and a full-grid dump has no + // reference columns. + int32_t num_candidates = 0; + read_scalar_or(matfile, "doppler_num_candidates", num_candidates, static_cast(0)); if (d_doppler_step == 0) { d_doppler_step = 1; } - if (doppler_narrowed) + bool columns_ok = false; + if (num_candidates > 0) { - d_num_doppler_bins = 2U; + d_num_doppler_bins = static_cast(num_candidates); + columns_ok = (var_->dims[1] >= d_num_doppler_bins) && (var_->dims[1] <= d_num_doppler_bins + 2U); + } + else if (doppler_narrowed) + { + d_num_doppler_bins = 1U; + columns_ok = (var_->dims[1] == 2U); } else { @@ -130,18 +141,19 @@ bool Acquisition_Dump_Reader::read_binary_acq() { d_num_doppler_bins = bins_floor; } + columns_ok = (var_->dims[1] == d_num_doppler_bins); } - if ((var_->dims[0] != d_samples_per_code) || (var_->dims[1] != d_num_doppler_bins)) + if ((var_->dims[0] != d_samples_per_code) || !columns_ok) { std::cout << "Invalid Acquisition dump file: dimension matrix error\n"; if (var_->dims[0] != d_samples_per_code) { std::cout << "Expected " << d_samples_per_code << " samples per code. Obtained " << var_->dims[0] << '\n'; } - if (var_->dims[1] != d_num_doppler_bins) + if (!columns_ok) { - std::cout << "Expected " << d_num_doppler_bins << " Doppler bins. Obtained " << var_->dims[1] << '\n'; + std::cout << "Expected " << d_num_doppler_bins << " Doppler candidate bins (plus up to 2 reference columns). Obtained " << var_->dims[1] << " columns\n"; } Mat_VarFree(var_); Mat_Close(matfile); @@ -149,10 +161,10 @@ bool Acquisition_Dump_Reader::read_binary_acq() } // 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}. + // metadata, for the candidate columns only (acq_grid is column-major, so + // they are its first d_num_doppler_bins * d_samples_per_code values). + // doppler(i) = -doppler_max + doppler_center + doppler_step * i is the same + // decoding pcps_acquisition::compute_statistics() uses. doppler.clear(); for (unsigned int doppler_index = 0; doppler_index < d_num_doppler_bins; doppler_index++) { diff --git a/tests/unit-tests/signal-processing-blocks/libs/acquisition_dump_reader.h b/tests/unit-tests/signal-processing-blocks/libs/acquisition_dump_reader.h index 28440d395..e1ceb9c0c 100644 --- a/tests/unit-tests/signal-processing-blocks/libs/acquisition_dump_reader.h +++ b/tests/unit-tests/signal-processing-blocks/libs/acquisition_dump_reader.h @@ -59,7 +59,7 @@ public: 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} + bool doppler_narrowed{}; //!< true if the dump was written by a narrowed search (fewer candidate bins than the full configured grid) unsigned int PRN{}; unsigned int num_dwells{}; uint64_t sample_counter{}; diff --git a/utils/matlab/plot_acq_grid.m b/utils/matlab/plot_acq_grid.m index 9052ab042..9039057a8 100644 --- a/utils/matlab/plot_acq_grid.m +++ b/utils/matlab/plot_acq_grid.m @@ -95,16 +95,38 @@ else % Acquisition dumps written before doppler_center was added were centered at 0 Hz. doppler_center = 0; end -[n_fft, n_dop_bins] = size(acq_grid); -[d_max, f_max] = find(acq_grid == max(max(acq_grid))); +[n_fft, n_columns] = size(acq_grid); +if isfield(dump, 'doppler_num_candidates') + % The leading doppler_num_candidates columns are the Doppler candidates. + % Up to two trailing columns are noise-reference bins placed outside the + % Doppler axis, which are never acquisition candidates. + n_dop_bins = min(double(dump.doppler_num_candidates), n_columns); +elseif isfield(dump, 'doppler_narrowed') && dump.doppler_narrowed + % Older narrowed dumps: one candidate plus one noise-reference column. + n_dop_bins = 1; +else + n_dop_bins = n_columns; +end +acq_grid = acq_grid(:, 1 : n_dop_bins); +[d_max, f_max] = find(acq_grid == max(max(acq_grid)), 1); freq = double(doppler_center) + (0 : n_dop_bins - 1) * double(doppler_step) - double(doppler_max); delay = (0 : n_fft - 1) / n_fft * n_chips; +if n_dop_bins > 1 + freq_limits = [min(freq) max(freq)]; +else + freq_limits = freq(1) + [-0.5 0.5] * max(abs(double(doppler_step)), 1); +end %% Plot data %--- Acquisition grid (3D) figure(1) -if(lite_view == false) +if(n_dop_bins == 1) + % A single Doppler bin (narrowed search) cannot be drawn as a surface. + plot3(freq(1) * ones(size(delay)), delay, acq_grid(:, 1)) + grid on + ylim([min(delay) max(delay)]) +elseif(lite_view == false) surf(freq, delay, acq_grid, 'FaceColor', 'interp', 'LineStyle', 'none') ylim([min(delay) max(delay)]) else @@ -114,15 +136,19 @@ else ylim([min(delay_interp) max(delay_interp)]) end xlabel('Doppler shift (Hz)') -xlim([min(freq) max(freq)]) +xlim(freq_limits) ylabel('Code delay (chips)') zlabel('Test Statistics') %--- Acquisition grid (2D) figure(2) subplot(2,1,1) -plot(freq, acq_grid(d_max, :)) -xlim([min(freq) max(freq)]) +if n_dop_bins > 1 + plot(freq, acq_grid(d_max, :)) +else + plot(freq, acq_grid(d_max, :), 'o') +end +xlim(freq_limits) xlabel('Doppler shift (Hz)') ylabel('Test statistics') title(['Fixed code delay to ' num2str((d_max - 1) / n_fft * n_chips) ' chips']) diff --git a/utils/python/plot_acq_grid.py b/utils/python/plot_acq_grid.py index b680e1e45..8271c5579 100755 --- a/utils/python/plot_acq_grid.py +++ b/utils/python/plot_acq_grid.py @@ -343,6 +343,7 @@ def configured_doppler_axis( file_path, doppler_center=0.0, doppler_narrowed=False, + exact=False, ): import numpy as np @@ -352,14 +353,12 @@ def configured_doppler_axis( if doppler_step == 0: raise ValueError(f"{file_path}: doppler_step must be non-zero") - if doppler_narrowed: - if n_dop_bins != 2: - print( - f"Warning: {file_path.name}: narrowed acquisition metadata " - f"expects 2 Doppler bins, but acq_grid has {n_dop_bins} " - "columns. Plotting the available columns from the stored " - "center and step." - ) + if exact or doppler_narrowed: + # Candidate column k decodes to + # doppler_center - doppler_max + k * doppler_step, the same rule + # pcps_acquisition uses. This is exact for dumps that store + # doppler_num_candidates, and for legacy narrowed dumps (a single + # candidate at doppler_center, stored with doppler_max = 0). return ( doppler_center - doppler_max @@ -437,19 +436,42 @@ def read_acquisition_dump(file_path, n_chips, positive_only): doppler_narrowed = bool( dump_scalar_value(data, "doppler_narrowed", 0) ) - acq_grid = data["acq_grid"][:] - expected_doppler_bins = ( - {2} - if doppler_narrowed - else expected_doppler_bin_counts(doppler_max, doppler_step) + # Dumps that store doppler_num_candidates may carry up to two + # trailing noise-reference columns after the candidate columns. + # Older narrowed dumps have exactly one candidate plus one reference + # column; older full-grid dumps have no reference columns. + doppler_num_candidates = dump_scalar_value( + data, "doppler_num_candidates" ) + acq_grid = data["acq_grid"][:] + if doppler_num_candidates is not None: + doppler_num_candidates = int(doppler_num_candidates) + expected_doppler_bins = set( + range(doppler_num_candidates, doppler_num_candidates + 3) + ) + elif doppler_narrowed: + expected_doppler_bins = {2} + else: + expected_doppler_bins = expected_doppler_bin_counts( + doppler_max, doppler_step + ) if ( acq_grid.shape[1] not in expected_doppler_bins and acq_grid.shape[0] in expected_doppler_bins ): acq_grid = acq_grid.T - n_fft, n_dop_bins = acq_grid.shape + n_fft, n_columns = acq_grid.shape + if doppler_num_candidates is not None: + n_dop_bins = max(1, min(doppler_num_candidates, n_columns)) + elif doppler_narrowed: + n_dop_bins = 1 + else: + n_dop_bins = n_columns + if n_dop_bins < n_columns: + # Noise-reference columns are not on the Doppler axis and are + # never acquisition candidates: leave them out of the plots. + acq_grid = acq_grid[:, :n_dop_bins] freq = configured_doppler_axis( doppler_max, doppler_step, @@ -457,6 +479,7 @@ def read_acquisition_dump(file_path, n_chips, positive_only): file_path, doppler_center, doppler_narrowed, + exact=doppler_num_candidates is not None, ) d_max, f_max = np.unravel_index(np.argmax(acq_grid), acq_grid.shape) delay = np.arange(n_fft) / n_fft * n_chips @@ -570,9 +593,27 @@ def plot_dump(file_path, fig_path, image_name_root, n_chips, args, metadata=None metadata = parse_dump_name(file_path) plot_title = acquisition_plot_title(metadata) if metadata else str(file_path) + # A narrowed search can have a single candidate Doppler bin: no surface + # can be drawn over one Doppler value, so draw it as a code-delay line. + single_doppler_bin = len(freq) < 2 + if single_doppler_bin: + half_step = max(abs(float(doppler_step)), 1.0) / 2 + freq_limits = [freq[0] - half_step, freq[0] + half_step] + else: + freq_limits = [min(freq), max(freq)] + fig = plt.figure() plt.gcf().canvas.manager.set_window_title(str(file_path)) - if not args.lite_view: + if single_doppler_bin: + ax = fig.add_subplot(111, projection="3d") + ax.plot( + np.full(len(delay), freq[0]), + delay, + acq_grid[:, 0], + linewidth=0.8, + ) + ax.set_ylim([min(delay), max(delay)]) + elif not args.lite_view: ax = fig.add_subplot(111, projection="3d") x_axis, y_axis = np.meshgrid(freq, delay, indexing="xy") # rcount/ccount default to 50, which would silently decimate large @@ -624,7 +665,7 @@ def plot_dump(file_path, fig_path, image_name_root, n_chips, args, metadata=None label="Peak", ) ax.set_xlabel("Doppler shift (Hz)") - ax.set_xlim([min(freq), max(freq)]) + ax.set_xlim(freq_limits) ax.set_ylabel("Code delay (chips)") ax.set_zlabel("Test Statistics") ax.set_title(plot_title) @@ -637,9 +678,13 @@ def plot_dump(file_path, fig_path, image_name_root, n_chips, args, metadata=None fig2, axes = plt.subplots(2, 1, figsize=(8, 6)) plt.gcf().canvas.manager.set_window_title(str(file_path)) fig2.suptitle(plot_title) - axes[0].plot(freq, acq_grid[d_max, :]) + axes[0].plot( + freq, + acq_grid[d_max, :], + marker="o" if single_doppler_bin else None, + ) axes[0].axvline(peak_doppler, color="red", linestyle="--", linewidth=0.8) - axes[0].set_xlim([min(freq), max(freq)]) + axes[0].set_xlim(freq_limits) axes[0].set_xlabel("Doppler shift (Hz)") axes[0].set_ylabel("Test statistics") axes[0].set_title(f"Doppler cut at code delay = {peak_delay:g} chips")