From b66d8066c3c99737ea5da1257b27f26acbc08b94 Mon Sep 17 00:00:00 2001 From: Javier Date: Tue, 18 Feb 2020 16:45:56 +0100 Subject: [PATCH] Adding a new RINEX observables testing tool to command line utilities set --- src/utils/CMakeLists.txt | 1 + src/utils/rinex_tools/CMakeLists.txt | 190 +++++ src/utils/rinex_tools/obsdiff.cc | 1104 +++++++++++++++++++++++++ src/utils/rinex_tools/obsdiff_flags.h | 46 ++ 4 files changed, 1341 insertions(+) create mode 100644 src/utils/rinex_tools/CMakeLists.txt create mode 100644 src/utils/rinex_tools/obsdiff.cc create mode 100644 src/utils/rinex_tools/obsdiff_flags.h diff --git a/src/utils/CMakeLists.txt b/src/utils/CMakeLists.txt index aeef9ad1c..2b122dfe3 100644 --- a/src/utils/CMakeLists.txt +++ b/src/utils/CMakeLists.txt @@ -8,6 +8,7 @@ # add_subdirectory(front-end-cal) +add_subdirectory(rinex_tools) if(ENABLE_UNIT_TESTING_EXTRA OR ENABLE_SYSTEM_TESTING_EXTRA OR ENABLE_FPGA) add_subdirectory(rinex2assist) diff --git a/src/utils/rinex_tools/CMakeLists.txt b/src/utils/rinex_tools/CMakeLists.txt new file mode 100644 index 000000000..c5a975c04 --- /dev/null +++ b/src/utils/rinex_tools/CMakeLists.txt @@ -0,0 +1,190 @@ +# Copyright (C) 2012-2020 (see AUTHORS file for a list of contributors) +# +# GNSS-SDR is a software-defined Global Navigation Satellite Systems receiver +# +# This file is part of GNSS-SDR. +# +# SPDX-License-Identifier: GPL-3.0-or-later +# + +############################################################################ +# Local installation of GPSTk https://github.com/SGL-UT/GPSTk +############################################################################ +find_package(GPSTK) +set_package_properties(GPSTK PROPERTIES + PURPOSE "Used in some Extra Tests." + TYPE REQUIRED +) +if(NOT GPSTK_FOUND OR ENABLE_OWN_GPSTK) + set_package_properties(GPSTK PROPERTIES + PURPOSE "GPSTk v${GNSSSDR_GPSTK_LOCAL_VERSION} will be automatically downloaded and built when doing '${CMAKE_MAKE_PROGRAM_PRETTY_NAME}'." + ) + if("${TOOLCHAIN_ARG}" STREQUAL "") + set(TOOLCHAIN_ARG "-DCMAKE_CXX_FLAGS=-Wno-deprecated") + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wno-deprecated") + endif() + set(GPSTK_BUILD_COMMAND "${CMAKE_MAKE_PROGRAM}") + if(CMAKE_GENERATOR STREQUAL Xcode) + set(GPSTK_BUILD_COMMAND "xcodebuild" "-configuration" "${CMAKE_BUILD_TYPE}") + endif() + include(GNUInstallDirs) + if(CMAKE_VERSION VERSION_LESS 3.2) + ExternalProject_Add(gpstk-${GNSSSDR_GPSTK_LOCAL_VERSION} + GIT_REPOSITORY https://github.com/SGL-UT/GPSTk + GIT_TAG v${GNSSSDR_GPSTK_LOCAL_VERSION} + SOURCE_DIR ${CMAKE_CURRENT_SOURCE_DIR}/../../thirdparty/gpstk-${GNSSSDR_GPSTK_LOCAL_VERSION} + BINARY_DIR ${CMAKE_CURRENT_BINARY_DIR}/../../gpstk-${GNSSSDR_GPSTK_LOCAL_VERSION} + CMAKE_ARGS ${GTEST_COMPILER} ${TOOLCHAIN_ARG} -DCMAKE_INSTALL_PREFIX=${CMAKE_SOURCE_DIR}/thirdparty/gpstk-${GNSSSDR_GPSTK_LOCAL_VERSION}/install -DBUILD_EXT=OFF -DBUILD_PYTHON=OFF -DCMAKE_BUILD_TYPE=${CMAKE_BUILD_TYPE} + BUILD_COMMAND ${GPSTK_BUILD_COMMAND} + UPDATE_COMMAND "" + PATCH_COMMAND "" + ) + else() + ExternalProject_Add(gpstk-${GNSSSDR_GPSTK_LOCAL_VERSION} + GIT_REPOSITORY https://github.com/SGL-UT/GPSTk + GIT_TAG v${GNSSSDR_GPSTK_LOCAL_VERSION} + SOURCE_DIR ${CMAKE_CURRENT_SOURCE_DIR}/../../thirdparty/gpstk-${GNSSSDR_GPSTK_LOCAL_VERSION} + BINARY_DIR ${CMAKE_CURRENT_BINARY_DIR}/../../gpstk-${GNSSSDR_GPSTK_LOCAL_VERSION} + CMAKE_ARGS ${GTEST_COMPILER} ${TOOLCHAIN_ARG} -DCMAKE_INSTALL_PREFIX=${CMAKE_SOURCE_DIR}/thirdparty/gpstk-${GNSSSDR_GPSTK_LOCAL_VERSION}/install -DBUILD_EXT=OFF -DBUILD_PYTHON=OFF -DCMAKE_BUILD_TYPE=${CMAKE_BUILD_TYPE} + BUILD_COMMAND ${GPSTK_BUILD_COMMAND} + BUILD_BYPRODUCTS ${CMAKE_SOURCE_DIR}/thirdparty/gpstk-${GNSSSDR_GPSTK_LOCAL_VERSION}/install/${CMAKE_INSTALL_LIBDIR}/${CMAKE_FIND_LIBRARY_PREFIXES}gpstk${CMAKE_SHARED_LIBRARY_SUFFIX} + UPDATE_COMMAND "" + PATCH_COMMAND "" + ) + endif() + set(GPSTK_INCLUDE_DIRS ${CMAKE_CURRENT_SOURCE_DIR}/../../thirdparty/gpstk-${GNSSSDR_GPSTK_LOCAL_VERSION}/install/include CACHE PATH "Local GPSTK headers") + set(GPSTK_LIBRARY ${CMAKE_CURRENT_SOURCE_DIR}/../../thirdparty/gpstk-${GNSSSDR_GPSTK_LOCAL_VERSION}/install/${CMAKE_INSTALL_LIBDIR}/${CMAKE_FIND_LIBRARY_PREFIXES}gpstk${CMAKE_SHARED_LIBRARY_SUFFIX}) + set(GPSTK_BINDIR ${CMAKE_CURRENT_SOURCE_DIR}/../../thirdparty/gpstk-${GNSSSDR_GPSTK_LOCAL_VERSION}/install/bin/) + add_definitions(-DGPSTK_BINDIR="${GPSTK_BINDIR}") + add_library(Gpstk::gpstk SHARED IMPORTED) + add_dependencies(Gpstk::gpstk gpstk-${GNSSSDR_GPSTK_LOCAL_VERSION}) + file(MAKE_DIRECTORY ${GPSTK_INCLUDE_DIRS}) + file(MAKE_DIRECTORY ${GPSTK_INCLUDE_DIRS}/gpstk) + set_target_properties(Gpstk::gpstk PROPERTIES + IMPORTED_LINK_INTERFACE_LANGUAGES "CXX" + IMPORTED_LOCATION "${GPSTK_LIBRARY}" + INTERFACE_INCLUDE_DIRECTORIES "${GPSTK_INCLUDE_DIRS};${GPSTK_INCLUDE_DIRS}/gpstk" + INTERFACE_LINK_LIBRARIES "${GPSTK_LIBRARY}" + ) +else() + set(GPSTK_BINDIR ${GPSTK_INCLUDE_DIR}/../bin/) + add_definitions(-DGPSTK_BINDIR="${GPSTK_BINDIR}") +endif() + + +find_package(Gnuplot) +if(GNUPLOT_FOUND) + add_definitions(-DGNUPLOT_EXECUTABLE="${GNUPLOT_EXECUTABLE}") +endif() +set_package_properties(Gnuplot PROPERTIES + URL "http://www.gnuplot.info" + PURPOSE "Used to generate plots in some tests." + TYPE OPTIONAL +) +if(GNUPLOT_FOUND AND GNUPLOT_VERSION_STRING) + set_package_properties(Gnuplot PROPERTIES + DESCRIPTION "A portable command-line driven graphing utility (found: v${GNUPLOT_VERSION_STRING})" + ) +else() + set_package_properties(Gnuplot PROPERTIES + DESCRIPTION "A portable command-line driven graphing utility" + ) +endif() + + +set(LIST_INCLUDE_DIRS + ${CMAKE_SOURCE_DIR}/src/tests/common-files +) + +include_directories(${LIST_INCLUDE_DIRS}) + +#set(obsdiff_SOURCES obsdiff.cc) +set(obsdiff_HEADERS obsdiff_flags.h) +source_group(Headers FILES ${obsdiff_HEADERS}) + + +if(ENABLE_CLANG_TIDY) + if(CLANG_TIDY_EXE) + set_target_properties(obsdiff_lib + PROPERTIES + CXX_CLANG_TIDY "${DO_CLANG_TIDY}" + ) + endif() +endif() + +add_executable(obsdiff ${CMAKE_CURRENT_SOURCE_DIR}/obsdiff.cc) + +################################################################################ +# Detect availability of std::filesystem and set C++ standard accordingly +################################################################################ +set(FILESYSTEM_FOUND FALSE) +if(NOT (GNURADIO_VERSION VERSION_LESS 3.8) AND LOG4CPP_READY_FOR_CXX17) + # Check if we have std::filesystem + if(NOT (CMAKE_VERSION VERSION_LESS 3.8)) + find_package(FILESYSTEM COMPONENTS Final Experimental) + set_package_properties(FILESYSTEM PROPERTIES + URL "https://en.cppreference.com/w/cpp/filesystem" + DESCRIPTION "Provides facilities for performing operations on file systems and their components" + PURPOSE "Work with paths, regular files, and directories." + TYPE OPTIONAL + ) + if(${FILESYSTEM_FOUND}) + if(CMAKE_VERSION VERSION_LESS 3.13) + set(CMAKE_CXX_STANDARD 17) + else() + set(CMAKE_CXX_STANDARD 20) + endif() + set(CMAKE_CXX_STANDARD_REQUIRED ON) + endif() + endif() +endif() + + + + +target_link_libraries(obsdiff + PUBLIC + Boost::headers + PRIVATE + Armadillo::armadillo + Threads::Threads + channel_libs + algorithms_libs + core_receiver + core_libs + Gflags::gflags + Glog::glog + Gpstk::gpstk +) + +target_compile_definitions(obsdiff + PUBLIC -DGNSS_SDR_VERSION="${VERSION}" + PUBLIC -DGNSSSDR_INSTALL_DIR="${CMAKE_INSTALL_PREFIX}" +) + +if(ENABLE_CLANG_TIDY) + if(CLANG_TIDY_EXE) + set_target_properties(obsdiff + PROPERTIES + CXX_CLANG_TIDY "${DO_CLANG_TIDY}" + ) + endif() +endif() + +add_custom_command(TARGET obsdiff POST_BUILD + COMMAND ${CMAKE_COMMAND} -E copy $ + ${CMAKE_SOURCE_DIR}/install/$) + +install(TARGETS obsdiff + RUNTIME DESTINATION bin + COMPONENT "obsdiff" +) + +find_program(GZIP gzip + /bin + /usr/bin + /usr/local/bin + /opt/local/bin + /sbin +) + diff --git a/src/utils/rinex_tools/obsdiff.cc b/src/utils/rinex_tools/obsdiff.cc new file mode 100644 index 000000000..40927787f --- /dev/null +++ b/src/utils/rinex_tools/obsdiff.cc @@ -0,0 +1,1104 @@ +/*! + * \file obs-diff.cc + * \brief This class implements a single difference and double difference comparison algorithm to evaluate receiver's performance at observable level + * \authors
    + *
  • Javier Arribas, 2020. jarribas(at)cttc.es + *
+ * + * + * ------------------------------------------------------------------------- + * + * Copyright (C) 2010-2020 (see AUTHORS file for a list of contributors) + * + * GNSS-SDR is a software defined Global Navigation + * Satellite Systems receiver + * + * This file is part of GNSS-SDR. + * + * GNSS-SDR is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * GNSS-SDR is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with GNSS-SDR. If not, see . + * + * ------------------------------------------------------------------------- + */ + +#include "MATH_CONSTANTS.h" +#include "gnuplot_i.h" +#include "obsdiff_flags.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#if HAS_STD_FILESYSTEM +#include +namespace fs = std::filesystem; +#else +#include +namespace fs = boost::filesystem; +#endif + +// Create the lists of GNSS satellites +std::set available_gps_prn = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, + 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, + 29, 30, 31, 32}; + +std::set available_sbas_prn = {123, 131, 135, 136, 138}; + +std::set available_galileo_prn = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, + 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, + 29, 30, 31, 32, 33, 34, 35, 36}; + + +bool file_exist(const char* fileName) +{ + std::ifstream infile(fileName); + return infile.good(); +} + +std::map ReadRinexObs(std::string rinex_file, char system, std::string signal) +{ + std::map obs_map; + if (not file_exist(rinex_file.c_str())) + { + std::cout << "Warning: RINEX Obs file " << rinex_file << " does not exist\n"; + return obs_map; + } + // Open and read reference RINEX observables file + try + { + gpstk::Rinex3ObsStream r_ref(rinex_file); + r_ref.exceptions(std::ios::failbit); + gpstk::Rinex3ObsData r_ref_data; + gpstk::Rinex3ObsHeader r_ref_header; + + gpstk::RinexDatum dataobj; + + r_ref >> r_ref_header; + + std::set PRN_set; + gpstk::SatID prn; + switch (system) + { + case 'G': + prn.system = gpstk::SatID::systemGPS; + PRN_set = available_gps_prn; + break; + case 'E': + prn.system = gpstk::SatID::systemGalileo; + PRN_set = available_galileo_prn; + break; + default: + prn.system = gpstk::SatID::systemGPS; + PRN_set = available_gps_prn; + } + + std::cout << "Reading RINEX OBS file " << rinex_file << " ...\n"; + while (r_ref >> r_ref_data) + { + for (std::set::iterator prn_it = PRN_set.begin(); prn_it != PRN_set.end(); ++prn_it) + { + prn.id = *prn_it; + gpstk::CommonTime time = r_ref_data.time; + double sow(static_cast(time).sow); + + auto pointer = r_ref_data.obs.find(prn); + + if (pointer != r_ref_data.obs.end()) + { + //insert next column + try + { + obs_map.at(prn.id); + } + catch (const std::out_of_range& oor) + { + obs_map[prn.id] = arma::mat(); + } + + arma::mat& obs_mat = obs_map[prn.id]; + obs_mat.insert_rows(obs_mat.n_rows, arma::zeros(1, 4)); + + if (strcmp("1C\0", signal.c_str()) == 0) + { + obs_mat.at(obs_mat.n_rows - 1, 0) = sow; + dataobj = r_ref_data.getObs(prn, "C1C", r_ref_header); + obs_mat.at(obs_mat.n_rows - 1, 1) = dataobj.data; //C1C P1 (psudorange L1) + dataobj = r_ref_data.getObs(prn, "D1C", r_ref_header); + obs_mat.at(obs_mat.n_rows - 1, 2) = dataobj.data; //D1C Carrier Doppler + dataobj = r_ref_data.getObs(prn, "L1C", r_ref_header); + obs_mat.at(obs_mat.n_rows - 1, 3) = dataobj.data; //L1C Carrier Phase + } + else if (strcmp("1B\0", signal.c_str()) == 0) + { + obs_mat.at(obs_mat.n_rows - 1, 0) = sow; + dataobj = r_ref_data.getObs(prn, "C1B", r_ref_header); + obs_mat.at(obs_mat.n_rows - 1, 1) = dataobj.data; + dataobj = r_ref_data.getObs(prn, "D1B", r_ref_header); + obs_mat.at(obs_mat.n_rows - 1, 2) = dataobj.data; + dataobj = r_ref_data.getObs(prn, "L1B", r_ref_header); + obs_mat.at(obs_mat.n_rows - 1, 3) = dataobj.data; + } + else if (strcmp("2S\0", signal.c_str()) == 0) //L2M + { + obs_mat.at(obs_mat.n_rows - 1, 0) = sow; + dataobj = r_ref_data.getObs(prn, "C2S", r_ref_header); + obs_mat.at(obs_mat.n_rows - 1, 1) = dataobj.data; + dataobj = r_ref_data.getObs(prn, "D2S", r_ref_header); + obs_mat.at(obs_mat.n_rows - 1, 2) = dataobj.data; + dataobj = r_ref_data.getObs(prn, "L2S", r_ref_header); + obs_mat.at(obs_mat.n_rows - 1, 3) = dataobj.data; + } + else if (strcmp("L5\0", signal.c_str()) == 0) + { + obs_mat.at(obs_mat.n_rows - 1, 0) = sow; + dataobj = r_ref_data.getObs(prn, "C5I", r_ref_header); + obs_mat.at(obs_mat.n_rows - 1, 1) = dataobj.data; + dataobj = r_ref_data.getObs(prn, "D5I", r_ref_header); + obs_mat.at(obs_mat.n_rows - 1, 2) = dataobj.data; + dataobj = r_ref_data.getObs(prn, "L5I", r_ref_header); + obs_mat.at(obs_mat.n_rows - 1, 3) = dataobj.data; + } + else if (strcmp("5X\0", signal.c_str()) == 0) //Simulator gives RINEX with E5a+E5b. Doppler and accumulated Carrier phase WILL differ + { + obs_mat.at(obs_mat.n_rows - 1, 0) = sow; + dataobj = r_ref_data.getObs(prn, "C8I", r_ref_header); + obs_mat.at(obs_mat.n_rows - 1, 1) = dataobj.data; + dataobj = r_ref_data.getObs(prn, "D8I", r_ref_header); + obs_mat.at(obs_mat.n_rows - 1, 2) = dataobj.data; + dataobj = r_ref_data.getObs(prn, "L8I", r_ref_header); + obs_mat.at(obs_mat.n_rows - 1, 3) = dataobj.data; + } + else + { + std::cout << "ReadRinexObs unknown signal requested: " << signal << std::endl; + return obs_map; + } + } + } + } // end while + } // End of 'try' block + catch (const gpstk::FFStreamError& e) + { + std::cout << e; + return obs_map; + } + catch (const gpstk::Exception& e) + { + std::cout << e; + return obs_map; + } + catch (const std::exception& e) + { + std::cout << "Exception: " << e.what(); + std::cout << "unknown error. I don't feel so well..." << std::endl; + return obs_map; + } + + return obs_map; +} + +bool save_mat_xy(std::vector& x, std::vector& y, std::string filename) +{ + try + { + // WRITE MAT FILE + mat_t* matfp; + matvar_t* matvar; + filename.append(".mat"); + //std::cout << "save_mat_xy write " << filename << std::endl; + matfp = Mat_CreateVer(filename.c_str(), nullptr, MAT_FT_MAT5); + if (reinterpret_cast(matfp) != nullptr) + { + size_t dims[2] = {1, x.size()}; + matvar = Mat_VarCreate("x", MAT_C_DOUBLE, MAT_T_DOUBLE, 2, dims, &x[0], 0); + Mat_VarWrite(matfp, matvar, MAT_COMPRESSION_ZLIB); // or MAT_COMPRESSION_NONE + Mat_VarFree(matvar); + + matvar = Mat_VarCreate("y", MAT_C_DOUBLE, MAT_T_DOUBLE, 2, dims, &y[0], 0); + Mat_VarWrite(matfp, matvar, MAT_COMPRESSION_ZLIB); // or MAT_COMPRESSION_NONE + Mat_VarFree(matvar); + } + else + { + std::cout << "save_mat_xy: error creating file" << std::endl; + } + Mat_Close(matfp); + return true; + } + catch (const std::exception& ex) + { + std::cout << "save_mat_xy: " << ex.what() << std::endl; + return false; + } +} + + +bool save_mat_x(std::vector* x, std::string filename) +{ + try + { + // WRITE MAT FILE + mat_t* matfp; + matvar_t* matvar; + filename.append(".mat"); + std::cout << "save_mat_x write " << filename << std::endl; + matfp = Mat_CreateVer(filename.c_str(), nullptr, MAT_FT_MAT5); + if (reinterpret_cast(matfp) != nullptr) + { + std::array dims{1, x->size()}; + matvar = Mat_VarCreate("x", MAT_C_DOUBLE, MAT_T_DOUBLE, 2, dims.data(), &x[0], 0); + Mat_VarWrite(matfp, matvar, MAT_COMPRESSION_ZLIB); // or MAT_COMPRESSION_NONE + Mat_VarFree(matvar); + } + else + { + std::cout << "save_mat_x: error creating file" << std::endl; + } + Mat_Close(matfp); + return true; + } + catch (const std::exception& ex) + { + std::cout << "save_mat_x: " << ex.what() << std::endl; + return false; + } +} + +void check_results_carrier_phase_double_diff( + arma::mat& true_ch0, + arma::mat& true_ch1, + arma::mat& measured_ch0, + arma::mat& measured_ch1, + const std::string& data_title) +{ + // 1. True value interpolation to match the measurement times + arma::vec measurement_time = measured_ch0.col(0); + + arma::vec true_ch0_carrier_phase_interp; + arma::vec true_ch1_carrier_phase_interp; + arma::interp1(true_ch0.col(0), true_ch0.col(3), measurement_time, true_ch0_carrier_phase_interp); + arma::interp1(true_ch1.col(0), true_ch1.col(3), measurement_time, true_ch1_carrier_phase_interp); + + arma::vec meas_ch1_carrier_phase_interp; + arma::interp1(measured_ch1.col(0), measured_ch1.col(3), measurement_time, meas_ch1_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_true_carrier_phase_cycles = (true_ch0_carrier_phase_interp - true_ch0_carrier_phase_interp(0)) - (true_ch1_carrier_phase_interp - true_ch1_carrier_phase_interp(0)); + arma::vec delta_measured_carrier_phase_cycles = (measured_ch0.col(3) - measured_ch0.col(3)(0)) - (meas_ch1_carrier_phase_interp - meas_ch1_carrier_phase_interp(0)); + + // remove NaN + arma::uvec NaN_in_true_data = arma::find_nonfinite(delta_true_carrier_phase_cycles); + arma::uvec NaN_in_measured_data = arma::find_nonfinite(delta_measured_carrier_phase_cycles); + + arma::mat tmp_mat = arma::conv_to::from(delta_true_carrier_phase_cycles); + tmp_mat.shed_rows(arma::join_cols(NaN_in_true_data, NaN_in_measured_data)); + delta_true_carrier_phase_cycles = tmp_mat.col(0); + + tmp_mat = arma::conv_to::from(delta_measured_carrier_phase_cycles); + tmp_mat.shed_rows(arma::join_cols(NaN_in_true_data, NaN_in_measured_data)); + delta_measured_carrier_phase_cycles = tmp_mat.col(0); + + tmp_mat = arma::conv_to::from(measurement_time); + tmp_mat.shed_rows(arma::join_cols(NaN_in_true_data, NaN_in_measured_data)); + measurement_time = tmp_mat.col(0); + + std::vector + time_vector(measurement_time.colptr(0), measurement_time.colptr(0) + measurement_time.n_rows); + + if (measurement_time.size() > 0) + { + //debug + // std::vector tmp_time_vec(measurement_time.colptr(0), + // measurement_time.colptr(0) + measurement_time.n_rows); + // std::vector tmp_vector_y6(delta_true_carrier_doppler_cycles.colptr(0), + // delta_true_carrier_doppler_cycles.colptr(0) + delta_true_carrier_doppler_cycles.n_rows); + // save_mat_xy(tmp_time_vec, tmp_vector_y6, std::string("true_delta_doppler")); + // std::vector tmp_vector_y7(delta_measured_carrier_doppler_cycles.colptr(0), + // delta_measured_carrier_doppler_cycles.colptr(0) + delta_measured_carrier_doppler_cycles.n_rows); + // save_mat_xy(tmp_time_vec, tmp_vector_y7, std::string("measured_delta_doppler")); + + // 2. RMSE + arma::vec err; + err = delta_measured_carrier_phase_cycles - delta_true_carrier_phase_cycles; + arma::vec err2 = arma::square(err); + double rmse = sqrt(arma::mean(err2)); + + // 3. Mean err and variance + double error_mean = arma::mean(err); + double error_var = arma::var(err); + + // 4. Peaks + double max_error = arma::max(err); + double min_error = arma::min(err); + + // 5. report + std::streamsize ss = std::cout.precision(); + std::cout << std::setprecision(10) << data_title << "Double diff Carrier Phase RMSE = " + << rmse << ", mean = " << error_mean + << ", stdev = " << sqrt(error_var) + << " (max,min) = " << max_error + << "," << min_error + << " [Cycles]" << std::endl; + std::cout.precision(ss); + + // plots + if (FLAGS_show_plots) + { + Gnuplot g3("linespoints"); + g3.set_title(data_title + "Double diff Carrier Phase error [Cycles]"); + g3.set_grid(); + g3.set_xlabel("Time [s]"); + g3.set_ylabel("Double diff Carrier Phase error [Cycles]"); + // conversion between arma::vec and std:vector + std::vector range_error_m(err.colptr(0), err.colptr(0) + err.n_rows); + g3.cmd("set key box opaque"); + g3.plot_xy(time_vector, range_error_m, + "Double diff Carrier Phase error"); + g3.set_legend(); + g3.savetops(data_title + "double_diff_carrier_phase_error"); + + g3.showonscreen(); // window output + } + // check results against the test tolerance + // ASSERT_LT(rmse, 0.25); + // ASSERT_LT(error_mean, 0.2); + // ASSERT_GT(error_mean, -0.2); + // ASSERT_LT(error_var, 0.5); + // ASSERT_LT(max_error, 0.5); + // ASSERT_GT(min_error, -0.5); + } + else + { + std::cout << "No valid data\n"; + } +} + +void check_results_carrier_phase_single_diff( + arma::mat& measured_ch0, + arma::mat& measured_ch1, + const std::string& data_title) +{ + // 1. True value interpolation to match the measurement times + arma::vec measurement_time = measured_ch0.col(0); + + arma::vec meas_ch1_carrier_phase_interp; + arma::interp1(measured_ch1.col(0), measured_ch1.col(3), measurement_time, meas_ch1_carrier_phase_interp); + // generate single 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 = (measured_ch0.col(3) - measured_ch0.col(3)(0)) - (meas_ch1_carrier_phase_interp - meas_ch1_carrier_phase_interp(0)); + + // remove NaN + arma::uvec NaN_in_measured_data = arma::find_nonfinite(delta_measured_carrier_phase_cycles); + + arma::mat tmp_mat = arma::conv_to::from(delta_measured_carrier_phase_cycles); + tmp_mat.shed_rows(NaN_in_measured_data); + delta_measured_carrier_phase_cycles = tmp_mat.col(0); + + tmp_mat = arma::conv_to::from(measurement_time); + tmp_mat.shed_rows(NaN_in_measured_data); + measurement_time = tmp_mat.col(0); + + std::vector + time_vector(measurement_time.colptr(0), measurement_time.colptr(0) + measurement_time.n_rows); + + if (measurement_time.size() > 0) + { + // 2. RMSE + arma::vec err; + err = delta_measured_carrier_phase_cycles; + arma::vec err2 = arma::square(err); + double rmse = sqrt(arma::mean(err2)); + + // 3. Mean err and variance + double error_mean = arma::mean(err); + double error_var = arma::var(err); + + // 4. Peaks + double max_error = arma::max(err); + double min_error = arma::min(err); + + // 5. report + std::streamsize ss = std::cout.precision(); + std::cout << std::setprecision(10) << data_title << "Single diff Carrier Phase RMSE = " + << rmse << ", mean = " << error_mean + << ", stdev = " << sqrt(error_var) + << " (max,min) = " << max_error + << "," << min_error + << " [Cycles]" << std::endl; + std::cout.precision(ss); + + // plots + if (FLAGS_show_plots) + { + Gnuplot g3("linespoints"); + g3.set_title(data_title + "Single diff Carrier Phase error [Cycles]"); + g3.set_grid(); + g3.set_xlabel("Time [s]"); + g3.set_ylabel("Single diff Carrier Phase error [Cycles]"); + // conversion between arma::vec and std:vector + std::vector range_error_m(err.colptr(0), err.colptr(0) + err.n_rows); + g3.cmd("set key box opaque"); + g3.plot_xy(time_vector, range_error_m, + "Single diff Carrier Phase error"); + g3.set_legend(); + g3.savetops(data_title + "single_diff_carrier_phase_error"); + + g3.showonscreen(); // window output + } + } + else + { + std::cout << "No valid data\n"; + } +} + +void check_results_carrier_doppler_double_diff( + arma::mat& true_ch0, + arma::mat& true_ch1, + arma::mat& measured_ch0, + arma::mat& measured_ch1, + const std::string& data_title) +{ + // 1. True value interpolation to match the measurement times + arma::vec measurement_time = measured_ch0.col(0); + + arma::vec true_ch0_carrier_doppler_interp; + arma::vec true_ch1_carrier_doppler_interp; + arma::interp1(true_ch0.col(0), true_ch0.col(2), measurement_time, true_ch0_carrier_doppler_interp); + arma::interp1(true_ch1.col(0), true_ch1.col(2), measurement_time, true_ch1_carrier_doppler_interp); + + arma::vec meas_ch1_carrier_doppler_interp; + arma::interp1(measured_ch1.col(0), measured_ch1.col(2), measurement_time, meas_ch1_carrier_doppler_interp); + // generate double difference carrier Doppler + arma::vec delta_true_carrier_doppler_cycles = true_ch0_carrier_doppler_interp - true_ch1_carrier_doppler_interp; + arma::vec delta_measured_carrier_doppler_cycles = measured_ch0.col(2) - meas_ch1_carrier_doppler_interp; + + // remove NaN + arma::uvec NaN_in_true_data = arma::find_nonfinite(delta_true_carrier_doppler_cycles); + arma::uvec NaN_in_measured_data = arma::find_nonfinite(delta_measured_carrier_doppler_cycles); + + arma::mat tmp_mat = arma::conv_to::from(delta_true_carrier_doppler_cycles); + tmp_mat.shed_rows(arma::join_cols(NaN_in_true_data, NaN_in_measured_data)); + delta_true_carrier_doppler_cycles = tmp_mat.col(0); + + tmp_mat = arma::conv_to::from(delta_measured_carrier_doppler_cycles); + tmp_mat.shed_rows(arma::join_cols(NaN_in_true_data, NaN_in_measured_data)); + delta_measured_carrier_doppler_cycles = tmp_mat.col(0); + + tmp_mat = arma::conv_to::from(measurement_time); + tmp_mat.shed_rows(arma::join_cols(NaN_in_true_data, NaN_in_measured_data)); + measurement_time = tmp_mat.col(0); + + std::vector + time_vector(measurement_time.colptr(0), measurement_time.colptr(0) + measurement_time.n_rows); + + if (measurement_time.size() > 0) + { + //debug + // std::vector tmp_time_vec(measurement_time.colptr(0), + // measurement_time.colptr(0) + measurement_time.n_rows); + // std::vector tmp_vector_y6(delta_true_carrier_doppler_cycles.colptr(0), + // delta_true_carrier_doppler_cycles.colptr(0) + delta_true_carrier_doppler_cycles.n_rows); + // save_mat_xy(tmp_time_vec, tmp_vector_y6, std::string("true_delta_doppler")); + // std::vector tmp_vector_y7(delta_measured_carrier_doppler_cycles.colptr(0), + // delta_measured_carrier_doppler_cycles.colptr(0) + delta_measured_carrier_doppler_cycles.n_rows); + // save_mat_xy(tmp_time_vec, tmp_vector_y7, std::string("measured_delta_doppler")); + + // 2. RMSE + arma::vec err; + err = delta_measured_carrier_doppler_cycles - delta_true_carrier_doppler_cycles; + arma::vec err2 = arma::square(err); + double rmse = sqrt(arma::mean(err2)); + + // 3. Mean err and variance + double error_mean = arma::mean(err); + double error_var = arma::var(err); + + // 4. Peaks + double max_error = arma::max(err); + double min_error = arma::min(err); + + // 5. report + std::streamsize ss = std::cout.precision(); + std::cout << std::setprecision(10) << data_title << "Double diff Carrier Doppler RMSE = " + << rmse << ", mean = " << error_mean + << ", stdev = " << sqrt(error_var) + << " (max,min) = " << max_error + << "," << min_error + << " [Hz]" << std::endl; + std::cout.precision(ss); + + // plots + if (FLAGS_show_plots) + { + Gnuplot g3("linespoints"); + g3.set_title(data_title + "Double diff Carrier Doppler error [Hz]"); + g3.set_grid(); + g3.set_xlabel("Time [s]"); + g3.set_ylabel("Double diff Carrier Doppler error [Hz]"); + // conversion between arma::vec and std:vector + std::vector range_error_m(err.colptr(0), err.colptr(0) + err.n_rows); + g3.cmd("set key box opaque"); + g3.plot_xy(time_vector, range_error_m, + "Double diff Carrier Doppler error"); + g3.set_legend(); + g3.savetops(data_title + "double_diff_carrier_doppler_error"); + + g3.showonscreen(); // window output + } + + // check results against the test tolerance + // ASSERT_LT(error_mean, 5); + // ASSERT_GT(error_mean, -5); + // // assuming PLL BW=35 + // ASSERT_LT(error_var, 250); + // ASSERT_LT(max_error, 100); + // ASSERT_GT(min_error, -100); + // ASSERT_LT(rmse, 30); + } + else + { + std::cout << "No valid data\n"; + } +} + +void check_results_carrier_doppler_single_diff( + arma::mat& measured_ch0, + arma::mat& measured_ch1, + const std::string& data_title) +{ + // 1. True value interpolation to match the measurement times + arma::vec measurement_time = measured_ch0.col(0); + + arma::vec meas_ch1_carrier_doppler_interp; + arma::interp1(measured_ch1.col(0), measured_ch1.col(2), measurement_time, meas_ch1_carrier_doppler_interp); + // generate single difference carrier Doppler + arma::vec delta_measured_carrier_doppler_cycles = measured_ch0.col(2) - meas_ch1_carrier_doppler_interp; + + // remove NaN + arma::uvec NaN_in_measured_data = arma::find_nonfinite(delta_measured_carrier_doppler_cycles); + + arma::mat tmp_mat = arma::conv_to::from(delta_measured_carrier_doppler_cycles); + tmp_mat.shed_rows(NaN_in_measured_data); + delta_measured_carrier_doppler_cycles = tmp_mat.col(0); + + tmp_mat = arma::conv_to::from(measurement_time); + tmp_mat.shed_rows(NaN_in_measured_data); + measurement_time = tmp_mat.col(0); + + std::vector + time_vector(measurement_time.colptr(0), measurement_time.colptr(0) + measurement_time.n_rows); + + if (measurement_time.size() > 0) + { + // 2. RMSE + arma::vec err; + err = delta_measured_carrier_doppler_cycles; + arma::vec err2 = arma::square(err); + double rmse = sqrt(arma::mean(err2)); + + // 3. Mean err and variance + double error_mean = arma::mean(err); + double error_var = arma::var(err); + + // 4. Peaks + double max_error = arma::max(err); + double min_error = arma::min(err); + + // 5. report + std::streamsize ss = std::cout.precision(); + std::cout << std::setprecision(10) << data_title << "Single diff Carrier Doppler RMSE = " + << rmse << ", mean = " << error_mean + << ", stdev = " << sqrt(error_var) + << " (max,min) = " << max_error + << "," << min_error + << " [Hz]" << std::endl; + std::cout.precision(ss); + + // plots + if (FLAGS_show_plots) + { + Gnuplot g3("linespoints"); + g3.set_title(data_title + "Single diff Carrier Doppler error [Hz]"); + g3.set_grid(); + g3.set_xlabel("Time [s]"); + g3.set_ylabel("Single diff Carrier Doppler error [Hz]"); + // conversion between arma::vec and std:vector + std::vector range_error_m(err.colptr(0), err.colptr(0) + err.n_rows); + g3.cmd("set key box opaque"); + g3.plot_xy(time_vector, range_error_m, + "Single diff Carrier Doppler error"); + g3.set_legend(); + g3.savetops(data_title + "single_diff_carrier_doppler_error"); + + g3.showonscreen(); // window output + } + } + else + { + std::cout << "No valid data\n"; + } +} + +void check_results_code_pseudorange_double_diff( + arma::mat& true_ch0, + arma::mat& true_ch1, + arma::mat& measured_ch0, + arma::mat& measured_ch1, + const std::string& data_title) +{ + // 1. True value interpolation to match the measurement times + arma::vec measurement_time = measured_ch0.col(0); + arma::vec true_ch0_obs_interp; + arma::vec true_ch1_obs_interp; + arma::interp1(true_ch0.col(0), true_ch0.col(1), measurement_time, true_ch0_obs_interp); + arma::interp1(true_ch1.col(0), true_ch1.col(1), measurement_time, true_ch1_obs_interp); + + arma::vec meas_ch1_obs_interp; + arma::interp1(measured_ch1.col(0), measured_ch1.col(1), measurement_time, meas_ch1_obs_interp); + // generate double difference carrier Doppler + arma::vec delta_true_obs = true_ch0_obs_interp - true_ch1_obs_interp; + arma::vec delta_measured_obs = measured_ch0.col(1) - meas_ch1_obs_interp; + + // remove NaN + arma::uvec NaN_in_true_data = arma::find_nonfinite(delta_true_obs); + arma::uvec NaN_in_measured_data = arma::find_nonfinite(delta_measured_obs); + + arma::mat tmp_mat = arma::conv_to::from(delta_true_obs); + tmp_mat.shed_rows(arma::join_cols(NaN_in_true_data, NaN_in_measured_data)); + delta_true_obs = tmp_mat.col(0); + + tmp_mat = arma::conv_to::from(delta_measured_obs); + tmp_mat.shed_rows(arma::join_cols(NaN_in_true_data, NaN_in_measured_data)); + delta_measured_obs = tmp_mat.col(0); + + tmp_mat = arma::conv_to::from(measurement_time); + tmp_mat.shed_rows(arma::join_cols(NaN_in_true_data, NaN_in_measured_data)); + measurement_time = tmp_mat.col(0); + + std::vector + time_vector(measurement_time.colptr(0), measurement_time.colptr(0) + measurement_time.n_rows); + + if (measurement_time.size() > 0) + { + //debug + // std::vector tmp_time_vec(measurement_time.colptr(0), + // measurement_time.colptr(0) + measurement_time.n_rows); + // std::vector tmp_vector_y6(delta_true_carrier_doppler_cycles.colptr(0), + // delta_true_carrier_doppler_cycles.colptr(0) + delta_true_carrier_doppler_cycles.n_rows); + // save_mat_xy(tmp_time_vec, tmp_vector_y6, std::string("true_delta_doppler")); + // std::vector tmp_vector_y7(delta_measured_carrier_doppler_cycles.colptr(0), + // delta_measured_carrier_doppler_cycles.colptr(0) + delta_measured_carrier_doppler_cycles.n_rows); + // save_mat_xy(tmp_time_vec, tmp_vector_y7, std::string("measured_delta_doppler")); + + // 2. RMSE + arma::vec err; + + err = delta_measured_obs - delta_true_obs; + + arma::vec err2 = arma::square(err); + double rmse = sqrt(arma::mean(err2)); + + // 3. Mean err and variance + double error_mean = arma::mean(err); + double error_var = arma::var(err); + + // 4. Peaks + double max_error = arma::max(err); + double min_error = arma::min(err); + + // 5. report + std::streamsize ss = std::cout.precision(); + std::cout << std::setprecision(10) << data_title << "Double diff Pseudorange RMSE = " + << rmse << ", mean = " << error_mean + << ", stdev = " << sqrt(error_var) + << " (max,min) = " << max_error + << "," << min_error + << " [meters]" << std::endl; + std::cout.precision(ss); + + // plots + if (FLAGS_show_plots) + { + Gnuplot g3("linespoints"); + g3.set_title(data_title + "Double diff Pseudorange error [m]"); + g3.set_grid(); + g3.set_xlabel("Time [s]"); + g3.set_ylabel("Double diff Pseudorange error [m]"); + // conversion between arma::vec and std:vector + std::vector range_error_m(err.colptr(0), err.colptr(0) + err.n_rows); + g3.cmd("set key box opaque"); + g3.plot_xy(time_vector, range_error_m, + "Double diff Pseudorrange error"); + g3.set_legend(); + g3.savetops(data_title + "double_diff_pseudorrange_error"); + + g3.showonscreen(); // window output + } + + // check results against the test tolerance + // ASSERT_LT(rmse, 3.0); + // ASSERT_LT(error_mean, 1.0); + // ASSERT_GT(error_mean, -1.0); + // ASSERT_LT(error_var, 10.0); + // ASSERT_LT(max_error, 10.0); + // ASSERT_GT(min_error, -10.0); + } + else + { + std::cout << "No valid data\n"; + } +} + +void check_results_code_pseudorange_single_diff( + arma::mat& measured_ch0, + arma::mat& measured_ch1, + const std::string& data_title) +{ + // 1. True value interpolation to match the measurement times + arma::vec measurement_time = measured_ch0.col(0); + + arma::vec meas_ch1_obs_interp; + arma::interp1(measured_ch1.col(0), measured_ch1.col(1), measurement_time, meas_ch1_obs_interp); + // generate single difference carrier Doppler + arma::vec delta_measured_obs = measured_ch0.col(1) - meas_ch1_obs_interp; + + // remove NaN + arma::uvec NaN_in_measured_data = arma::find_nonfinite(delta_measured_obs); + + arma::mat tmp_mat = arma::conv_to::from(delta_measured_obs); + tmp_mat.shed_rows(NaN_in_measured_data); + delta_measured_obs = tmp_mat.col(0); + + tmp_mat = arma::conv_to::from(measurement_time); + tmp_mat.shed_rows(NaN_in_measured_data); + measurement_time = tmp_mat.col(0); + + std::vector + time_vector(measurement_time.colptr(0), measurement_time.colptr(0) + measurement_time.n_rows); + + if (measurement_time.size() > 0) + { + // 2. RMSE + arma::vec err; + + err = delta_measured_obs; + + arma::vec err2 = arma::square(err); + double rmse = sqrt(arma::mean(err2)); + + // 3. Mean err and variance + double error_mean = arma::mean(err); + double error_var = arma::var(err); + + // 4. Peaks + double max_error = arma::max(err); + double min_error = arma::min(err); + + // 5. report + std::streamsize ss = std::cout.precision(); + std::cout << std::setprecision(10) << data_title << "Single diff Pseudorange RMSE = " + << rmse << ", mean = " << error_mean + << ", stdev = " << sqrt(error_var) + << " (max,min) = " << max_error + << "," << min_error + << " [meters]" << std::endl; + std::cout.precision(ss); + + // plots + if (FLAGS_show_plots) + { + Gnuplot g3("linespoints"); + g3.set_title(data_title + "Single diff Pseudorange error [m]"); + g3.set_grid(); + g3.set_xlabel("Time [s]"); + g3.set_ylabel("Single diff Pseudorange error [m]"); + // conversion between arma::vec and std:vector + std::vector range_error_m(err.colptr(0), err.colptr(0) + err.n_rows); + g3.cmd("set key box opaque"); + g3.plot_xy(time_vector, range_error_m, + "Single diff Pseudorrange error"); + g3.set_legend(); + g3.savetops(data_title + "Single_diff_pseudorrange_error"); + + g3.showonscreen(); // window output + } + } + else + { + std::cout << "No valid data\n"; + } +} + + +void RINEX_doublediff_dupli_sat() +{ + // special test mode for duplicated satellites + // read rinex receiver-under-test observations + std::map test_obs = ReadRinexObs(FLAGS_test_rinex_obs, 'G', std::string("1C\0")); + if (test_obs.size() == 0) + { + return; + } + // Cut measurement initial transitory of the measurements + double initial_transitory_s = FLAGS_skip_obs_transitory_s; + std::cout << "Skipping initial transitory of " << initial_transitory_s << " [s]" << std::endl; + arma::uvec index; + for (std::map::iterator it = test_obs.begin(); it != test_obs.end(); ++it) + { + index = arma::find(it->second.col(0) >= (it->second.col(0)(0) + initial_transitory_s), 1, "first"); + if ((!index.empty()) and (index(0) > 0)) + { + it->second.shed_rows(0, index(0)); + } + } + + + std::vector prn_pairs; + std::stringstream ss(FLAGS_dupli_sat_prns); + unsigned int i; + while (ss >> i) + { + 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 + { + for (unsigned int n = 0; n < prn_pairs.size(); n = n + 2) + { + //compute double differences + if (test_obs.find(prn_pairs.at(n)) != test_obs.end() and test_obs.find(prn_pairs.at(n + 1)) != test_obs.end()) + { + std::cout << "Computing single difference observables for duplicated SV pairs..." << std::endl; + std::cout << "SD=OBS_ROVER(SV" << prn_pairs.at(n) << ")-OBS_ROVER(SV" << prn_pairs.at(n + 1) << ")" << std::endl; + + check_results_code_pseudorange_single_diff(test_obs.at(prn_pairs.at(n)), + test_obs.at(prn_pairs.at(n + 1)), + "SD=OBS(SV" + std::to_string(prn_pairs.at(n)) + ")-OBS(SV" + std::to_string(prn_pairs.at(n + 1)) + ") "); + + check_results_carrier_phase_single_diff(test_obs.at(prn_pairs.at(n)), + test_obs.at(prn_pairs.at(n + 1)), + "SD=OBS(SV" + std::to_string(prn_pairs.at(n)) + ")-OBS(SV" + std::to_string(prn_pairs.at(n + 1)) + ") "); + + check_results_carrier_doppler_single_diff(test_obs.at(prn_pairs.at(n)), + test_obs.at(prn_pairs.at(n + 1)), + "SD=OBS(SV" + std::to_string(prn_pairs.at(n)) + ")-OBS(SV" + std::to_string(prn_pairs.at(n + 1)) + ") "); + } + else + { + std::cout << "Satellite ID " << prn_pairs.at(n) << " and/or " << prn_pairs.at(n + 1) << " not found in RINEX file\n"; + } + } + } +} +void RINEX_doublediff() +{ + // read rinex reference observations + std::map ref_obs = ReadRinexObs(FLAGS_ref_rinex_obs, 'G', std::string("1C\0")); + // read rinex receiver-under-test observations + std::map test_obs = ReadRinexObs(FLAGS_test_rinex_obs, 'G', std::string("1C\0")); + + if (ref_obs.size() == 0 or test_obs.size() == 0) + { + return; + } + // Cut measurement initial transitory of the measurements + double initial_transitory_s = FLAGS_skip_obs_transitory_s; + std::cout << "Skipping initial transitory of " << initial_transitory_s << " [s]" << std::endl; + arma::uvec index; + for (std::map::iterator it = test_obs.begin(); it != test_obs.end(); ++it) + { + index = arma::find(it->second.col(0) >= (it->second.col(0)(0) + initial_transitory_s), 1, "first"); + if ((!index.empty()) and (index(0) > 0)) + { + it->second.shed_rows(0, index(0)); + } + } + + // Cut observation vectors ends to the shortest one (base or rover) + arma::colvec ref_obs_time = ref_obs.begin()->second.col(0); + arma::colvec test_obs_time = test_obs.begin()->second.col(0); + + if (ref_obs_time.back() < test_obs_time.back()) + { + //there are more rover observations than base observations + //cut rover vector + std::cout << "Cutting rover observations vector end.." << std::endl; + arma::uvec index2; + for (std::map::iterator it = test_obs.begin(); it != test_obs.end(); ++it) + { + index = arma::find(it->second.col(0) >= ref_obs_time.back(), 1, "first"); + if ((!index.empty()) and (index(0) > 0)) + { + it->second.shed_rows(index(0), it->second.n_rows - 1); + } + } + } + else + { + //there are more base observations than rover observations + //cut base vector + std::cout << "Cutting base observations vector end.." << std::endl; + for (std::map::iterator it = ref_obs.begin(); it != ref_obs.end(); ++it) + { + index = arma::find(it->second.col(0) >= test_obs_time.back(), 1, "first"); + if ((!index.empty()) and (index(0) > 0)) + { + it->second.shed_rows(index(0), it->second.n_rows - 1); + } + } + } + + // also skip last seconds of the observations (some artifacts are present in some RINEX endings) + ref_obs_time = ref_obs.begin()->second.col(0); + test_obs_time = test_obs.begin()->second.col(0); + + double skip_ends_s = FLAGS_skip_obs_ends_s; + std::cout << "Skipping last " << skip_ends_s << " [s] of observations" << std::endl; + for (std::map::iterator it = test_obs.begin(); it != test_obs.end(); ++it) + { + index = arma::find(it->second.col(0) >= (test_obs_time.back() - skip_ends_s), 1, "first"); + if ((!index.empty()) and (index(0) > 0)) + { + it->second.shed_rows(index(0), it->second.n_rows - 1); + } + } + for (std::map::iterator it = ref_obs.begin(); it != ref_obs.end(); ++it) + { + index = arma::find(it->second.col(0) >= (ref_obs_time.back() - skip_ends_s), 1, "first"); + if ((!index.empty()) and (index(0) > 0)) + { + it->second.shed_rows(index(0), it->second.n_rows - 1); + } + } + + // Save observations in .mat files + std::cout << "Saving RAW observables inputs to .mat files...\n"; + for (std::map::iterator it = ref_obs.begin(); it != ref_obs.end(); ++it) + { + // std::cout << it->first << " => " << it->second.n_rows << '\n'; + // std::cout << it->first << " has NaN values: " << it->second.has_nan() << '\n'; + //debug + std::vector tmp_time_vec(it->second.col(0).colptr(0), + it->second.col(0).colptr(0) + it->second.n_rows); + std::vector tmp_vector_y6(it->second.col(2).colptr(0), + it->second.col(2).colptr(0) + it->second.n_rows); + save_mat_xy(tmp_time_vec, tmp_vector_y6, std::string("ref_doppler_sat" + std::to_string(it->first))); + } + for (std::map::iterator it = test_obs.begin(); it != test_obs.end(); ++it) + { + // std::cout << it->first << " => " << it->second.n_rows << '\n'; + // std::cout << it->first << " has NaN values: " << it->second.has_nan() << '\n'; + //debug + std::vector tmp_time_vec(it->second.col(0).colptr(0), + it->second.col(0).colptr(0) + it->second.n_rows); + std::vector tmp_vector_y6(it->second.col(2).colptr(0), + it->second.col(2).colptr(0) + it->second.n_rows); + save_mat_xy(tmp_time_vec, tmp_vector_y6, std::string("measured_doppler_sat" + std::to_string(it->first))); + } + + // select reference satellite + std::set PRN_set = available_gps_prn; + double abs_min_doppler = 1e6; + int ref_sat_id = 1; + for (std::set::iterator ref_prn_it = PRN_set.begin(); ref_prn_it != PRN_set.end(); ++ref_prn_it) + { + if (ref_obs.find(*ref_prn_it) != ref_obs.end() and test_obs.find(*ref_prn_it) != test_obs.end()) + { + if (fabs(test_obs.at(*ref_prn_it).at(0, 2)) < abs_min_doppler) + { + abs_min_doppler = fabs(test_obs.at(*ref_prn_it).at(0, 2)); + ref_sat_id = *ref_prn_it; + } + } + } + + //compute double differences + if (ref_obs.find(ref_sat_id) != ref_obs.end() and test_obs.find(ref_sat_id) != test_obs.end()) + { + std::cout << "Using reference satellite SV " << ref_sat_id << " with minimum abs Doppler of " << abs_min_doppler << " [Hz]" << std::endl; + for (std::set::iterator current_prn_it = PRN_set.begin(); current_prn_it != PRN_set.end(); ++current_prn_it) + { + int current_sat_id = *current_prn_it; + if (current_sat_id != ref_sat_id) + { + if (ref_obs.find(current_sat_id) != ref_obs.end() and test_obs.find(current_sat_id) != test_obs.end()) + { + std::cout << "Computing double difference observables for SV " << current_sat_id << std::endl; + std::cout << "DD=(OBS_ROVER(SV" << current_sat_id << ")-OBS_ROVER(SV" << ref_sat_id << "))" + << "-(OBS_BASE(SV" << current_sat_id << ")-OBS_BASE(SV" << ref_sat_id << "))" << std::endl; + + check_results_code_pseudorange_double_diff(ref_obs.at(ref_sat_id), + ref_obs.at(current_sat_id), + test_obs.at(ref_sat_id), + test_obs.at(current_sat_id), + "PRN " + std::to_string(current_sat_id) + " "); + + check_results_carrier_phase_double_diff(ref_obs.at(ref_sat_id), + ref_obs.at(current_sat_id), + test_obs.at(ref_sat_id), + test_obs.at(current_sat_id), + "PRN " + std::to_string(current_sat_id) + " "); + + check_results_carrier_doppler_double_diff(ref_obs.at(ref_sat_id), + ref_obs.at(current_sat_id), + test_obs.at(ref_sat_id), + test_obs.at(current_sat_id), + "PRN " + std::to_string(current_sat_id) + " "); + } + } + } + } + else + { + std::cout << "Satellite ID " << ref_sat_id << " not found in both RINEX files\n"; + } +} + +int main(int argc, char** argv) +{ + std::cout << "Running RINEX observables difference tool..." << std::endl; + google::ParseCommandLineFlags(&argc, &argv, true); + google::InitGoogleLogging(argv[0]); + if (FLAGS_dupli_sat) + { + RINEX_doublediff_dupli_sat(); + } + else + { + RINEX_doublediff(); + } + google::ShutDownCommandLineFlags(); + return 0; +} diff --git a/src/utils/rinex_tools/obsdiff_flags.h b/src/utils/rinex_tools/obsdiff_flags.h new file mode 100644 index 000000000..a659350e2 --- /dev/null +++ b/src/utils/rinex_tools/obsdiff_flags.h @@ -0,0 +1,46 @@ +/*! + * \file rinex_observable_tests_flags.h + * \brief Helper file for unit testing + * \author Javier Arribas, 2019. jarribas(at)cttc.es + * + * ------------------------------------------------------------------------- + * + * Copyright (C) 2010-2019 (see AUTHORS file for a list of contributors) + * + * GNSS-SDR is a software defined Global Navigation + * Satellite Systems receiver + * + * This file is part of GNSS-SDR. + * + * GNSS-SDR is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * GNSS-SDR is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with GNSS-SDR. If not, see . + * + * ------------------------------------------------------------------------- + */ + +#ifndef GNSS_SDR_OBS_DIFF_FLAGS_H_ +#define GNSS_SDR_OBS_DIFF_FLAGS_H_ + +#include +#include + +DEFINE_double(skip_obs_transitory_s, 30.0, "Skip the initial observable outputs to avoid transitory results [s]"); +DEFINE_double(skip_obs_ends_s, 5.0, "Skip the lasts observable outputs to avoid transitory results [s]"); +DEFINE_bool(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(dupli_sat, false, "Enable special observable test mode where the scenario contains duplicated satellite orbits"); +DEFINE_string(dupli_sat_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)"); +DEFINE_string(ref_rinex_obs, "reference.obs", "Filename of reference RINEX navigation file"); +DEFINE_string(test_rinex_obs, "test.obs", "Filename of test RINEX navigation file"); + +#endif