diff --git a/src/algorithms/PVT/adapters/rtklib_pvt.cc b/src/algorithms/PVT/adapters/rtklib_pvt.cc index dc8cf894a..13bcd6faf 100644 --- a/src/algorithms/PVT/adapters/rtklib_pvt.cc +++ b/src/algorithms/PVT/adapters/rtklib_pvt.cc @@ -506,12 +506,12 @@ bool RtklibPvt::save_assistance_to_XML() if (eph_map.empty() == false) { + std::ofstream ofs; try { std::ofstream ofs(eph_xml_filename_.c_str(), std::ofstream::trunc | std::ofstream::out); boost::archive::xml_oarchive xml(ofs); xml << boost::serialization::make_nvp("GNSS-SDR_ephemeris_map", eph_map); - ofs.close(); LOG(INFO) << "Saved GPS L1 Ephemeris map data"; } catch (const std::exception& e) @@ -519,13 +519,13 @@ bool RtklibPvt::save_assistance_to_XML() LOG(WARNING) << e.what(); return false; } - return true; // return variable (true == succeeded) } else { LOG(WARNING) << "Failed to save Ephemeris, map is empty"; return false; } + return true; // return variable (true == succeeded) } diff --git a/src/algorithms/PVT/gnuradio_blocks/rtklib_pvt_cc.cc b/src/algorithms/PVT/gnuradio_blocks/rtklib_pvt_cc.cc index 3a1228ec4..2f9bb177d 100644 --- a/src/algorithms/PVT/gnuradio_blocks/rtklib_pvt_cc.cc +++ b/src/algorithms/PVT/gnuradio_blocks/rtklib_pvt_cc.cc @@ -400,12 +400,12 @@ rtklib_pvt_cc::~rtklib_pvt_cc() if (d_ls_pvt->gps_cnav_ephemeris_map.empty() == false) { + std::ofstream ofs; try { - std::ofstream ofs(file_name.c_str(), std::ofstream::trunc | std::ofstream::out); + ofs.open(file_name.c_str(), std::ofstream::trunc | std::ofstream::out); boost::archive::xml_oarchive xml(ofs); xml << boost::serialization::make_nvp("GNSS-SDR_ephemeris_map", d_ls_pvt->gps_cnav_ephemeris_map); - ofs.close(); LOG(INFO) << "Saved GPS L2CM or L5 Ephemeris map data"; } catch (std::exception& e) @@ -423,12 +423,12 @@ rtklib_pvt_cc::~rtklib_pvt_cc() if (d_ls_pvt->gps_ephemeris_map.empty() == false) { + std::ofstream ofs; try { - std::ofstream ofs(file_name.c_str(), std::ofstream::trunc | std::ofstream::out); + ofs.open(file_name.c_str(), std::ofstream::trunc | std::ofstream::out); boost::archive::xml_oarchive xml(ofs); xml << boost::serialization::make_nvp("GNSS-SDR_ephemeris_map", d_ls_pvt->gps_ephemeris_map); - ofs.close(); LOG(INFO) << "Saved GPS L1 CA Ephemeris map data"; } catch (const std::exception& e) @@ -446,12 +446,12 @@ rtklib_pvt_cc::~rtklib_pvt_cc() if (d_ls_pvt->galileo_ephemeris_map.empty() == false) { + std::ofstream ofs; try { - std::ofstream ofs(file_name.c_str(), std::ofstream::trunc | std::ofstream::out); + ofs.open(file_name.c_str(), std::ofstream::trunc | std::ofstream::out); boost::archive::xml_oarchive xml(ofs); xml << boost::serialization::make_nvp("GNSS-SDR_ephemeris_map", d_ls_pvt->galileo_ephemeris_map); - ofs.close(); LOG(INFO) << "Saved Galileo E1 Ephemeris map data"; } catch (const std::exception& e) @@ -469,12 +469,12 @@ rtklib_pvt_cc::~rtklib_pvt_cc() if (d_ls_pvt->glonass_gnav_ephemeris_map.empty() == false) { + std::ofstream ofs; try { - std::ofstream ofs(file_name.c_str(), std::ofstream::trunc | std::ofstream::out); + ofs.open(file_name.c_str(), std::ofstream::trunc | std::ofstream::out); boost::archive::xml_oarchive xml(ofs); xml << boost::serialization::make_nvp("GNSS-SDR_ephemeris_map", d_ls_pvt->glonass_gnav_ephemeris_map); - ofs.close(); LOG(INFO) << "Saved GLONASS GNAV Ephemeris map data"; } catch (std::exception& e) @@ -515,12 +515,12 @@ bool rtklib_pvt_cc::save_gnss_synchro_map_xml(const std::string file_name) { if (gnss_observables_map.empty() == false) { + std::ofstream ofs; try { - std::ofstream ofs(file_name.c_str(), std::ofstream::trunc | std::ofstream::out); + ofs.open(file_name.c_str(), std::ofstream::trunc | std::ofstream::out); boost::archive::xml_oarchive xml(ofs); xml << boost::serialization::make_nvp("GNSS-SDR_gnss_synchro_map", gnss_observables_map); - ofs.close(); LOG(INFO) << "Saved gnss_sychro map data"; } catch (std::exception& e) @@ -537,17 +537,17 @@ bool rtklib_pvt_cc::save_gnss_synchro_map_xml(const std::string file_name) } } + bool rtklib_pvt_cc::load_gnss_synchro_map_xml(const std::string file_name) { - //load from xml (boost serialize) + // load from xml (boost serialize) + std::ifstream ifs; try { - std::ifstream ifs(file_name.c_str(), std::ifstream::binary | std::ifstream::in); + ifs.open(file_name.c_str(), std::ifstream::binary | std::ifstream::in); boost::archive::xml_iarchive xml(ifs); gnss_observables_map.clear(); xml >> boost::serialization::make_nvp("GNSS-SDR_gnss_synchro_map", gnss_observables_map); - ifs.close(); - return true; //std::cout << "Loaded gnss_synchro map data with " << gnss_synchro_map.size() << " pseudoranges" << std::endl; } catch (std::exception& e) @@ -555,6 +555,7 @@ bool rtklib_pvt_cc::load_gnss_synchro_map_xml(const std::string file_name) std::cout << e.what() << "File: " << file_name; return false; } + return true; } diff --git a/src/core/libs/gnss_sdr_supl_client.cc b/src/core/libs/gnss_sdr_supl_client.cc index fdd08b9a9..4d3668adc 100644 --- a/src/core/libs/gnss_sdr_supl_client.cc +++ b/src/core/libs/gnss_sdr_supl_client.cc @@ -384,7 +384,6 @@ bool gnss_sdr_supl_client::load_ephemeris_xml(const std::string file_name) LOG(WARNING) << e.what() << "File: " << file_name; return false; } - ifs.close(); return true; } @@ -399,7 +398,6 @@ bool gnss_sdr_supl_client::save_ephemeris_map_xml(const std::string file_name, s ofs.open(file_name.c_str(), std::ofstream::trunc | std::ofstream::out); boost::archive::xml_oarchive xml(ofs); xml << boost::serialization::make_nvp("GNSS-SDR_ephemeris_map", eph_map); - LOG(INFO) << "Saved Ephemeris map data"; } catch (std::exception& e) @@ -407,14 +405,13 @@ bool gnss_sdr_supl_client::save_ephemeris_map_xml(const std::string file_name, s LOG(WARNING) << e.what(); return false; } - ofs.close(); - return true; } else { LOG(WARNING) << "Failed to save Ephemeris, map is empty"; return false; } + return true; } @@ -426,7 +423,6 @@ bool gnss_sdr_supl_client::load_utc_xml(const std::string file_name) ifs.open(file_name.c_str(), std::ifstream::binary | std::ifstream::in); boost::archive::xml_iarchive xml(ifs); xml >> boost::serialization::make_nvp("GNSS-SDR_utc_map", this->gps_utc); - LOG(INFO) << "Loaded UTC model data"; } catch (std::exception& e) @@ -434,7 +430,6 @@ bool gnss_sdr_supl_client::load_utc_xml(const std::string file_name) LOG(WARNING) << e.what() << "File: " << file_name; return false; } - ifs.close(); return true; } @@ -456,14 +451,13 @@ bool gnss_sdr_supl_client::save_utc_map_xml(const std::string file_name, std::ma LOG(WARNING) << e.what(); return false; } - ofs.close(); - return true; } else { LOG(WARNING) << "Failed to save UTC model, map is empty"; return false; } + return true; } @@ -482,7 +476,6 @@ bool gnss_sdr_supl_client::load_iono_xml(const std::string file_name) LOG(WARNING) << e.what() << "File: " << file_name; return false; } - ifs.close(); return true; } @@ -504,14 +497,13 @@ bool gnss_sdr_supl_client::save_iono_map_xml(const std::string file_name, std::m LOG(WARNING) << e.what(); return false; } - ofs.close(); - return true; } else { LOG(WARNING) << "Failed to save IONO model, map is empty"; return false; } + return true; } @@ -530,7 +522,6 @@ bool gnss_sdr_supl_client::load_ref_time_xml(const std::string file_name) LOG(WARNING) << e.what() << "File: " << file_name; return false; } - ifs.close(); return true; } @@ -553,14 +544,13 @@ bool gnss_sdr_supl_client::save_ref_time_map_xml(const std::string file_name, st LOG(WARNING) << e.what(); return false; } - ofs.close(); - return true; } else { LOG(WARNING) << "Failed to save Ref Time, map is empty"; return false; } + return true; } @@ -579,7 +569,6 @@ bool gnss_sdr_supl_client::load_ref_location_xml(const std::string file_name) LOG(WARNING) << e.what() << "File: " << file_name; return false; } - ifs.close(); return true; } @@ -601,12 +590,11 @@ bool gnss_sdr_supl_client::save_ref_location_map_xml(const std::string file_name LOG(WARNING) << e.what(); return false; } - ofs.close(); - return true; } else { LOG(WARNING) << "Failed to save Ref Location, map is empty"; return false; } + return true; } diff --git a/src/tests/common-files/observable_tests_flags.h b/src/tests/common-files/observable_tests_flags.h index e20253f34..37d2bab99 100644 --- a/src/tests/common-files/observable_tests_flags.h +++ b/src/tests/common-files/observable_tests_flags.h @@ -37,5 +37,6 @@ DEFINE_double(skip_obs_transitory_s, 30.0, "Skip the initial observable outputs to avoid transitory results [s]"); DEFINE_bool(compute_single_diffs, false, "Compute also the single difference errors for Accumulated Carrier Phase and Carrier Doppler (requires LO synchronization between receivers)"); DEFINE_bool(compare_with_5X, false, "Compare the E5a Doppler and Carrier Phases with the E5 full bw in RINEX (expect discrepancy due to the center frequencies differences"); - +DEFINE_bool(duplicated_satellites_test, false, "Enable special observable test mode where the scenario contains duplicated satellite orbits"); +DEFINE_string(duplicated_satellites_prns, "1,2,3,4", "List of duplicated satellites PRN pairs (i.e. 1,2,3,4 indicates that the PRNs 1,2 share the same orbit. The same applies for PRNs 3,4)"); #endif diff --git a/src/tests/test_main.cc b/src/tests/test_main.cc index b62e1b428..79627ae9a 100644 --- a/src/tests/test_main.cc +++ b/src/tests/test_main.cc @@ -161,7 +161,6 @@ DECLARE_string(log_dir); #include "unit-tests/signal-processing-blocks/pvt/rtcm_printer_test.cc" #include "unit-tests/signal-processing-blocks/pvt/rinex_printer_test.cc" #include "unit-tests/signal-processing-blocks/pvt/nmea_printer_test.cc" -#include "unit-tests/signal-processing-blocks/pvt/rtklib_solver_test.cc" #if EXTRA_TESTS @@ -172,6 +171,7 @@ DECLARE_string(log_dir); #include "unit-tests/signal-processing-blocks/tracking/gps_l1_ca_dll_pll_tracking_test.cc" #include "unit-tests/signal-processing-blocks/tracking/gps_l1_ca_kf_tracking_test.cc" #include "unit-tests/signal-processing-blocks/tracking/tracking_pull-in_test.cc" +#include "unit-tests/signal-processing-blocks/pvt/rtklib_solver_test.cc" #if ENABLE_FPGA #include "unit-tests/signal-processing-blocks/tracking/tracking_pull-in_test_fpga.cc" #endif diff --git a/src/tests/unit-tests/signal-processing-blocks/observables/hybrid_observables_test.cc b/src/tests/unit-tests/signal-processing-blocks/observables/hybrid_observables_test.cc index 81bb788ee..be4f2fb4f 100644 --- a/src/tests/unit-tests/signal-processing-blocks/observables/hybrid_observables_test.cc +++ b/src/tests/unit-tests/signal-processing-blocks/observables/hybrid_observables_test.cc @@ -236,6 +236,11 @@ public: arma::mat& measured_ch1, std::string data_title); + void check_results_duplicated_satellite( + arma::mat& measured_sat1, + arma::mat& measured_sat2, + std::string data_title); + HybridObservablesTest() { factory = std::make_shared(); @@ -987,6 +992,207 @@ void HybridObservablesTest::check_results_carrier_doppler( ASSERT_LT(rmse_ch0, 30); } +void HybridObservablesTest::check_results_duplicated_satellite( + arma::mat& measured_sat1, + arma::mat& measured_sat2, + std::string data_title) +{ + //1. True value interpolation to match the measurement times + + double t0 = measured_sat1(0, 0); + int size1 = measured_sat1.col(0).n_rows; + double t1 = measured_sat1(size1 - 1, 0); + arma::vec t = arma::linspace(t0, t1, floor((t1 - t0) * 1e3)); + //conversion between arma::vec and std:vector + arma::vec t_from_start = arma::linspace(0, t1 - t0, floor((t1 - t0) * 1e3)); + std::vector time_vector(t_from_start.colptr(0), t_from_start.colptr(0) + t_from_start.n_rows); + //Doppler + arma::vec meas_sat1_doppler_interp; + arma::interp1(measured_sat1.col(0), measured_sat1.col(2), t, meas_sat1_doppler_interp); + arma::vec meas_sat2_doppler_interp; + arma::interp1(measured_sat2.col(0), measured_sat2.col(2), t, meas_sat2_doppler_interp); + + //Carrier Phase + arma::vec meas_sat1_carrier_phase_interp; + arma::vec meas_sat2_carrier_phase_interp; + arma::interp1(measured_sat1.col(0), measured_sat1.col(3), t, meas_sat1_carrier_phase_interp); + arma::interp1(measured_sat2.col(0), measured_sat2.col(3), t, meas_sat2_carrier_phase_interp); + + // generate double difference accumulated carrier phases + //compute error without the accumulated carrier phase offsets (which depends on the receiver starting time) + arma::vec delta_measured_carrier_phase_cycles = (meas_sat1_carrier_phase_interp - meas_sat1_carrier_phase_interp(0)) - (meas_sat2_carrier_phase_interp - meas_sat2_carrier_phase_interp(0)); + + //Pseudoranges + arma::vec meas_sat1_dist_interp; + arma::vec meas_sat2_dist_interp; + arma::interp1(measured_sat1.col(0), measured_sat1.col(4), t, meas_sat1_dist_interp); + arma::interp1(measured_sat2.col(0), measured_sat2.col(4), t, meas_sat2_dist_interp); + // generate delta pseudoranges + arma::vec delta_measured_dist_m = meas_sat1_dist_interp - meas_sat2_dist_interp; + + //Carrier Doppler error + //2. RMSE + arma::vec err_ch0_hz; + + //compute error + err_ch0_hz = meas_sat1_doppler_interp - meas_sat2_doppler_interp; + + arma::vec err2_ch0 = arma::square(err_ch0_hz); + double rmse_ch0 = sqrt(arma::mean(err2_ch0)); + + //3. Mean err and variance + double error_mean_ch0 = arma::mean(err_ch0_hz); + double error_var_ch0 = arma::var(err_ch0_hz); + + // 4. Peaks + double max_error_ch0 = arma::max(err_ch0_hz); + double min_error_ch0 = arma::min(err_ch0_hz); + + //5. report + std::streamsize ss = std::cout.precision(); + std::cout << std::setprecision(10) << data_title << "Carrier Doppler RMSE = " + << rmse_ch0 << ", mean = " << error_mean_ch0 + << ", stdev = " << sqrt(error_var_ch0) + << " (max,min) = " << max_error_ch0 + << "," << min_error_ch0 + << " [Hz]" << std::endl; + std::cout.precision(ss); + + //plots + if (FLAGS_show_plots) + { + Gnuplot g3("linespoints"); + g3.set_title(data_title + "Carrier Doppler error [Hz]"); + g3.set_grid(); + g3.set_xlabel("Time [s]"); + g3.set_ylabel("Carrier Doppler error [Hz]"); + //conversion between arma::vec and std:vector + std::vector error_vec(err_ch0_hz.colptr(0), err_ch0_hz.colptr(0) + err_ch0_hz.n_rows); + g3.cmd("set key box opaque"); + g3.plot_xy(time_vector, error_vec, + "Carrier Doppler error"); + g3.set_legend(); + g3.savetops(data_title + "Carrier_doppler_error"); + + g3.showonscreen(); // window output + } + + //check results against the test tolerance + ASSERT_LT(error_mean_ch0, 5); + ASSERT_GT(error_mean_ch0, -5); + //assuming PLL BW=35 + ASSERT_LT(error_var_ch0, 250); + ASSERT_LT(max_error_ch0, 100); + ASSERT_GT(min_error_ch0, -100); + ASSERT_LT(rmse_ch0, 30); + + //Carrier Phase error + //2. RMSE + arma::vec err_carrier_phase; + + err_carrier_phase = delta_measured_carrier_phase_cycles; + arma::vec err2_carrier_phase = arma::square(err_carrier_phase); + double rmse_carrier_phase = sqrt(arma::mean(err2_carrier_phase)); + + //3. Mean err and variance + double error_mean_carrier_phase = arma::mean(err_carrier_phase); + double error_var_carrier_phase = arma::var(err_carrier_phase); + + // 4. Peaks + double max_error_carrier_phase = arma::max(err_carrier_phase); + double min_error_carrier_phase = arma::min(err_carrier_phase); + + //5. report + ss = std::cout.precision(); + std::cout << std::setprecision(10) << data_title << "Carrier Phase RMSE = " + << rmse_carrier_phase << ", mean = " << error_mean_carrier_phase + << ", stdev = " << sqrt(error_var_carrier_phase) + << " (max,min) = " << max_error_carrier_phase + << "," << min_error_carrier_phase + << " [Cycles]" << std::endl; + std::cout.precision(ss); + + //plots + if (FLAGS_show_plots) + { + Gnuplot g3("linespoints"); + g3.set_title(data_title + "Carrier Phase error [Cycles]"); + g3.set_grid(); + g3.set_xlabel("Time [s]"); + g3.set_ylabel("Carrier Phase error [Cycles]"); + //conversion between arma::vec and std:vector + std::vector range_error_m(err_carrier_phase.colptr(0), err_carrier_phase.colptr(0) + err_carrier_phase.n_rows); + g3.cmd("set key box opaque"); + g3.plot_xy(time_vector, range_error_m, + "Carrier Phase error"); + g3.set_legend(); + g3.savetops(data_title + "duplicated_satellite_carrier_phase_error"); + + g3.showonscreen(); // window output + } + + //check results against the test tolerance + ASSERT_LT(rmse_carrier_phase, 0.25); + ASSERT_LT(error_mean_carrier_phase, 0.2); + ASSERT_GT(error_mean_carrier_phase, -0.2); + ASSERT_LT(error_var_carrier_phase, 0.5); + ASSERT_LT(max_error_carrier_phase, 0.5); + ASSERT_GT(min_error_carrier_phase, -0.5); + + //Pseudorange error + //2. RMSE + arma::vec err_pseudorange; + + err_pseudorange = delta_measured_dist_m; + arma::vec err2_pseudorange = arma::square(err_pseudorange); + double rmse_pseudorange = sqrt(arma::mean(err2_pseudorange)); + + //3. Mean err and variance + double error_mean_pseudorange = arma::mean(err_pseudorange); + double error_var_pseudorange = arma::var(err_pseudorange); + + // 4. Peaks + double max_error_pseudorange = arma::max(err_pseudorange); + double min_error_pseudorange = arma::min(err_pseudorange); + + //5. report + ss = std::cout.precision(); + std::cout << std::setprecision(10) << data_title << "Pseudorange RMSE = " + << rmse_pseudorange << ", mean = " << error_mean_pseudorange + << ", stdev = " << sqrt(error_var_pseudorange) + << " (max,min) = " << max_error_pseudorange + << "," << min_error_pseudorange + << " [meters]" << std::endl; + std::cout.precision(ss); + + //plots + if (FLAGS_show_plots) + { + Gnuplot g3("linespoints"); + g3.set_title(data_title + "Pseudorange error [m]"); + g3.set_grid(); + g3.set_xlabel("Time [s]"); + g3.set_ylabel("Pseudorange error [m]"); + //conversion between arma::vec and std:vector + std::vector range_error_m(err_pseudorange.colptr(0), err_pseudorange.colptr(0) + err_pseudorange.n_rows); + g3.cmd("set key box opaque"); + g3.plot_xy(time_vector, range_error_m, + "Pseudorrange error"); + g3.set_legend(); + g3.savetops(data_title + "duplicated_satellite_pseudorrange_error"); + + g3.showonscreen(); // window output + } + + //check results against the test tolerance + ASSERT_LT(rmse_pseudorange, 3.0); + ASSERT_LT(error_mean_pseudorange, 1.0); + ASSERT_GT(error_mean_pseudorange, -1.0); + ASSERT_LT(error_var_pseudorange, 10.0); + ASSERT_LT(max_error_pseudorange, 10.0); + ASSERT_GT(min_error_pseudorange, -10.0); +} + bool HybridObservablesTest::save_mat_xy(std::vector& x, std::vector& y, std::string filename) { try @@ -1535,7 +1741,6 @@ TEST_F(HybridObservablesTest, ValidationOfResults) } } - //Cut measurement tail zeros arma::uvec index; for (unsigned int n = 0; n < measured_obs_vec.size(); n++) @@ -1548,9 +1753,7 @@ TEST_F(HybridObservablesTest, ValidationOfResults) } //Cut measurement initial transitory of the measurements - double initial_transitory_s = FLAGS_skip_obs_transitory_s; - for (unsigned int n = 0; n < measured_obs_vec.size(); n++) { index = arma::find(measured_obs_vec.at(n).col(0) >= (measured_obs_vec.at(n)(0, 0) + initial_transitory_s), 1, "first"); @@ -1567,83 +1770,128 @@ TEST_F(HybridObservablesTest, ValidationOfResults) } - //Correct the clock error using true values (it is not possible for a receiver to correct - //the receiver clock offset error at the observables level because it is required the - //decoding of the ephemeris data and solve the PVT equations) - - //Find the reference satellite (the nearest) and compute the receiver time offset at observable level - double min_pr = std::numeric_limits::max(); - unsigned int min_pr_ch_id = 0; - for (unsigned int n = 0; n < measured_obs_vec.size(); n++) + if (FLAGS_duplicated_satellites_test) { - if (epoch_counters_vec.at(n) > 10) //discard non-valid channels + //special test mode for duplicated satellites + std::vector prn_pairs; + std::stringstream ss(FLAGS_duplicated_satellites_prns); + int i; + while (ss >> i) { - { - if (measured_obs_vec.at(n)(0, 4) < min_pr) - { - min_pr = measured_obs_vec.at(n)(0, 4); - min_pr_ch_id = n; - } - } + prn_pairs.push_back(i); + if (ss.peek() == ',') + ss.ignore(); + } + + if (prn_pairs.size() % 2 != 0) + { + std::cout << "Test settings error: duplicated_satellites_prns are even\n"; } else { - std::cout << "PRN " << gnss_synchro_vec.at(n).PRN << " has NO observations!\n"; + for (unsigned int n = 0; n < prn_pairs.size(); n = n + 2) + { + int sat1_ch_id = -1; + int sat2_ch_id = -1; + for (unsigned int ch = 0; ch < measured_obs_vec.size(); ch++) + { + if (epoch_counters_vec.at(ch) > 10) //discard non-valid channels + { + if (gnss_synchro_vec.at(ch).PRN == prn_pairs.at(n)) + { + sat1_ch_id = ch; + } + if (gnss_synchro_vec.at(ch).PRN == prn_pairs.at(n + 1)) + { + sat2_ch_id = ch; + } + } + } + if (sat1_ch_id != -1 and sat2_ch_id != -1) + { + //compute single differences for the duplicated satellite + + check_results_duplicated_satellite( + measured_obs_vec.at(sat1_ch_id), + measured_obs_vec.at(sat2_ch_id), + "Duplicated sat [CH " + std::to_string(sat1_ch_id) + "," + std::to_string(sat2_ch_id) + "] PRNs " + std::to_string(gnss_synchro_vec.at(sat1_ch_id).PRN) + "," + std::to_string(gnss_synchro_vec.at(sat2_ch_id).PRN) + " "); + } + else + { + std::cout << "Satellites PRNs " << prn_pairs.at(n) << "and " << prn_pairs.at(n) << " not found\n"; + } + } } } - - arma::vec receiver_time_offset_ref_channel_s; - receiver_time_offset_ref_channel_s = (true_obs_vec.at(min_pr_ch_id).col(1)(0) - measured_obs_vec.at(min_pr_ch_id).col(4)(0)) / GPS_C_m_s; - std::cout << "Ref. channel initial Receiver time offset " << receiver_time_offset_ref_channel_s(0) * 1e3 << " [ms]" << std::endl; - - for (unsigned int n = 0; n < measured_obs_vec.size(); n++) + else { - //debug save to .mat - std::vector tmp_vector_x(true_obs_vec.at(n).col(0).colptr(0), - true_obs_vec.at(n).col(0).colptr(0) + true_obs_vec.at(n).col(0).n_rows); - std::vector tmp_vector_y(true_obs_vec.at(n).col(1).colptr(0), - true_obs_vec.at(n).col(1).colptr(0) + true_obs_vec.at(n).col(1).n_rows); - save_mat_xy(tmp_vector_x, tmp_vector_y, std::string("true_pr_ch_" + std::to_string(n))); + //normal mode - std::vector tmp_vector_x2(measured_obs_vec.at(n).col(0).colptr(0), - measured_obs_vec.at(n).col(0).colptr(0) + measured_obs_vec.at(n).col(0).n_rows); - std::vector tmp_vector_y2(measured_obs_vec.at(n).col(4).colptr(0), - measured_obs_vec.at(n).col(4).colptr(0) + measured_obs_vec.at(n).col(4).n_rows); - save_mat_xy(tmp_vector_x2, tmp_vector_y2, std::string("measured_pr_ch_" + std::to_string(n))); + //Correct the clock error using true values (it is not possible for a receiver to correct + //the receiver clock offset error at the observables level because it is required the + //decoding of the ephemeris data and solve the PVT equations) - std::vector tmp_vector_x3(true_obs_vec.at(n).col(0).colptr(0), - true_obs_vec.at(n).col(0).colptr(0) + true_obs_vec.at(n).col(0).n_rows); - std::vector tmp_vector_y3(true_obs_vec.at(n).col(2).colptr(0), - true_obs_vec.at(n).col(2).colptr(0) + true_obs_vec.at(n).col(2).n_rows); - save_mat_xy(tmp_vector_x3, tmp_vector_y3, std::string("true_doppler_ch_" + std::to_string(n))); - - std::vector tmp_vector_x4(measured_obs_vec.at(n).col(0).colptr(0), - measured_obs_vec.at(n).col(0).colptr(0) + measured_obs_vec.at(n).col(0).n_rows); - std::vector tmp_vector_y4(measured_obs_vec.at(n).col(2).colptr(0), - measured_obs_vec.at(n).col(2).colptr(0) + measured_obs_vec.at(n).col(2).n_rows); - save_mat_xy(tmp_vector_x4, tmp_vector_y4, std::string("measured_doppler_ch_" + std::to_string(n))); - - - if (epoch_counters_vec.at(n) > 10) //discard non-valid channels + //Find the reference satellite (the nearest) and compute the receiver time offset at observable level + double min_pr = std::numeric_limits::max(); + unsigned int min_pr_ch_id = 0; + for (unsigned int n = 0; n < measured_obs_vec.size(); n++) { - arma::vec true_TOW_ref_ch_s = true_obs_vec.at(min_pr_ch_id).col(0) - receiver_time_offset_ref_channel_s(0); - arma::vec true_TOW_ch_s = true_obs_vec.at(n).col(0) - receiver_time_offset_ref_channel_s(0); - //Compare measured observables - if (min_pr_ch_id != n) + if (epoch_counters_vec.at(n) > 10) //discard non-valid channels { - check_results_code_pseudorange(true_obs_vec.at(n), - true_obs_vec.at(min_pr_ch_id), - true_TOW_ch_s, - true_TOW_ref_ch_s, - measured_obs_vec.at(n), - measured_obs_vec.at(min_pr_ch_id), - "[CH " + std::to_string(n) + "] PRN " + std::to_string(gnss_synchro_vec.at(n).PRN) + " "); + { + if (measured_obs_vec.at(n)(0, 4) < min_pr) + { + min_pr = measured_obs_vec.at(n)(0, 4); + min_pr_ch_id = n; + } + } + } + else + { + std::cout << "PRN " << gnss_synchro_vec.at(n).PRN << " has NO observations!\n"; + } + } - //Do not compare E5a with E5 RINEX due to the Doppler frequency discrepancy caused by the different center frequencies - //E5a_fc=1176.45e6, E5b_fc=1207.14e6, E5_fc=1191.795e6; - if (strcmp("5X\0", gnss_synchro_vec.at(n).Signal) != 0 or FLAGS_compare_with_5X) + arma::vec receiver_time_offset_ref_channel_s; + receiver_time_offset_ref_channel_s = (true_obs_vec.at(min_pr_ch_id).col(1)(0) - measured_obs_vec.at(min_pr_ch_id).col(4)(0)) / GPS_C_m_s; + std::cout << "Ref. channel initial Receiver time offset " << receiver_time_offset_ref_channel_s(0) * 1e3 << " [ms]" << std::endl; + + for (unsigned int n = 0; n < measured_obs_vec.size(); n++) + { + //debug save to .mat + std::vector tmp_vector_x(true_obs_vec.at(n).col(0).colptr(0), + true_obs_vec.at(n).col(0).colptr(0) + true_obs_vec.at(n).col(0).n_rows); + std::vector tmp_vector_y(true_obs_vec.at(n).col(1).colptr(0), + true_obs_vec.at(n).col(1).colptr(0) + true_obs_vec.at(n).col(1).n_rows); + save_mat_xy(tmp_vector_x, tmp_vector_y, std::string("true_pr_ch_" + std::to_string(n))); + + std::vector tmp_vector_x2(measured_obs_vec.at(n).col(0).colptr(0), + measured_obs_vec.at(n).col(0).colptr(0) + measured_obs_vec.at(n).col(0).n_rows); + std::vector tmp_vector_y2(measured_obs_vec.at(n).col(4).colptr(0), + measured_obs_vec.at(n).col(4).colptr(0) + measured_obs_vec.at(n).col(4).n_rows); + save_mat_xy(tmp_vector_x2, tmp_vector_y2, std::string("measured_pr_ch_" + std::to_string(n))); + + std::vector tmp_vector_x3(true_obs_vec.at(n).col(0).colptr(0), + true_obs_vec.at(n).col(0).colptr(0) + true_obs_vec.at(n).col(0).n_rows); + std::vector tmp_vector_y3(true_obs_vec.at(n).col(2).colptr(0), + true_obs_vec.at(n).col(2).colptr(0) + true_obs_vec.at(n).col(2).n_rows); + save_mat_xy(tmp_vector_x3, tmp_vector_y3, std::string("true_doppler_ch_" + std::to_string(n))); + + std::vector tmp_vector_x4(measured_obs_vec.at(n).col(0).colptr(0), + measured_obs_vec.at(n).col(0).colptr(0) + measured_obs_vec.at(n).col(0).n_rows); + std::vector tmp_vector_y4(measured_obs_vec.at(n).col(2).colptr(0), + measured_obs_vec.at(n).col(2).colptr(0) + measured_obs_vec.at(n).col(2).n_rows); + save_mat_xy(tmp_vector_x4, tmp_vector_y4, std::string("measured_doppler_ch_" + std::to_string(n))); + + + if (epoch_counters_vec.at(n) > 10) //discard non-valid channels + { + arma::vec true_TOW_ref_ch_s = true_obs_vec.at(min_pr_ch_id).col(0) - receiver_time_offset_ref_channel_s(0); + arma::vec true_TOW_ch_s = true_obs_vec.at(n).col(0) - receiver_time_offset_ref_channel_s(0); + //Compare measured observables + if (min_pr_ch_id != n) { - check_results_carrier_phase_double_diff(true_obs_vec.at(n), + check_results_code_pseudorange(true_obs_vec.at(n), true_obs_vec.at(min_pr_ch_id), true_TOW_ch_s, true_TOW_ref_ch_s, @@ -1651,34 +1899,47 @@ TEST_F(HybridObservablesTest, ValidationOfResults) measured_obs_vec.at(min_pr_ch_id), "[CH " + std::to_string(n) + "] PRN " + std::to_string(gnss_synchro_vec.at(n).PRN) + " "); - check_results_carrier_doppler_double_diff(true_obs_vec.at(n), - true_obs_vec.at(min_pr_ch_id), + //Do not compare E5a with E5 RINEX due to the Doppler frequency discrepancy caused by the different center frequencies + //E5a_fc=1176.45e6, E5b_fc=1207.14e6, E5_fc=1191.795e6; + if (strcmp("5X\0", gnss_synchro_vec.at(n).Signal) != 0 or FLAGS_compare_with_5X) + { + check_results_carrier_phase_double_diff(true_obs_vec.at(n), + true_obs_vec.at(min_pr_ch_id), + true_TOW_ch_s, + true_TOW_ref_ch_s, + measured_obs_vec.at(n), + measured_obs_vec.at(min_pr_ch_id), + "[CH " + std::to_string(n) + "] PRN " + std::to_string(gnss_synchro_vec.at(n).PRN) + " "); + + check_results_carrier_doppler_double_diff(true_obs_vec.at(n), + true_obs_vec.at(min_pr_ch_id), + true_TOW_ch_s, + true_TOW_ref_ch_s, + measured_obs_vec.at(n), + measured_obs_vec.at(min_pr_ch_id), + "[CH " + std::to_string(n) + "] PRN " + std::to_string(gnss_synchro_vec.at(n).PRN) + " "); + } + } + else + { + std::cout << "[CH " << std::to_string(n) << "] PRN " << std::to_string(gnss_synchro_vec.at(n).PRN) << " is the reference satellite" << std::endl; + } + if (FLAGS_compute_single_diffs) + { + check_results_carrier_phase(true_obs_vec.at(n), + true_TOW_ch_s, + measured_obs_vec.at(n), + "[CH " + std::to_string(n) + "] PRN " + std::to_string(gnss_synchro_vec.at(n).PRN) + " "); + check_results_carrier_doppler(true_obs_vec.at(n), true_TOW_ch_s, - true_TOW_ref_ch_s, measured_obs_vec.at(n), - measured_obs_vec.at(min_pr_ch_id), "[CH " + std::to_string(n) + "] PRN " + std::to_string(gnss_synchro_vec.at(n).PRN) + " "); } } else { - std::cout << "[CH " << std::to_string(n) << "] PRN " << std::to_string(gnss_synchro_vec.at(n).PRN) << " is the reference satellite" << std::endl; + std::cout << "PRN " << gnss_synchro_vec.at(n).PRN << " has NO observations!\n"; } - if (FLAGS_compute_single_diffs) - { - check_results_carrier_phase(true_obs_vec.at(n), - true_TOW_ch_s, - measured_obs_vec.at(n), - "[CH " + std::to_string(n) + "] PRN " + std::to_string(gnss_synchro_vec.at(n).PRN) + " "); - check_results_carrier_doppler(true_obs_vec.at(n), - true_TOW_ch_s, - measured_obs_vec.at(n), - "[CH " + std::to_string(n) + "] PRN " + std::to_string(gnss_synchro_vec.at(n).PRN) + " "); - } - } - else - { - std::cout << "PRN " << gnss_synchro_vec.at(n).PRN << " has NO observations!\n"; } } std::cout << "Test completed in " << elapsed_seconds.count() << " [s]" << std::endl;