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")