| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | // This file is part of INSTINCT, the INS Toolkit for Integrated | ||
| 2 | // Navigation Concepts and Training by the Institute of Navigation of | ||
| 3 | // the University of Stuttgart, Germany. | ||
| 4 | // | ||
| 5 | // This Source Code Form is subject to the terms of the Mozilla Public | ||
| 6 | // License, v. 2.0. If a copy of the MPL was not distributed with this | ||
| 7 | // file, You can obtain one at https://mozilla.org/MPL/2.0/. | ||
| 8 | |||
| 9 | #include "AllanDeviation.hpp" | ||
| 10 | #include <algorithm> | ||
| 11 | #include <cstddef> | ||
| 12 | #include <cstdint> | ||
| 13 | #include <imgui.h> | ||
| 14 | #include <mutex> | ||
| 15 | |||
| 16 | #include "util/Logger.hpp" | ||
| 17 | |||
| 18 | #include <Eigen/src/Core/util/Meta.h> | ||
| 19 | #include <fmt/core.h> | ||
| 20 | #include "internal/FlowManager.hpp" | ||
| 21 | |||
| 22 | #include "internal/gui/widgets/HelpMarker.hpp" | ||
| 23 | #include "internal/gui/widgets/imgui_ex.hpp" | ||
| 24 | #include "implot.h" | ||
| 25 | |||
| 26 | #include "NodeData/IMU/ImuObs.hpp" | ||
| 27 | |||
| 28 | #include <chrono> | ||
| 29 | |||
| 30 | #include <unsupported/Eigen/MatrixFunctions> | ||
| 31 | |||
| 32 | 114 | NAV::AllanDeviation::AllanDeviation() | |
| 33 |
3/6✓ Branch 1 taken 114 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 114 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 114 times.
✗ Branch 9 not taken.
|
114 | : Node(typeStatic()) |
| 34 | { | ||
| 35 | LOG_TRACE("{}: called", name); | ||
| 36 | |||
| 37 | 114 | _hasConfig = true; | |
| 38 | 114 | _lockConfigDuringRun = false; | |
| 39 | 114 | _guiConfigDefaultWindowSize = { 630, 530 }; | |
| 40 | |||
| 41 |
4/8✓ Branch 1 taken 114 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 114 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 114 times.
✓ Branch 9 taken 114 times.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
|
342 | CreateInputPin("ImuObs", Pin::Type::Flow, { NAV::ImuObs::type() }, &AllanDeviation::receiveImuObs); |
| 42 | 228 | } | |
| 43 | |||
| 44 | 228 | NAV::AllanDeviation::~AllanDeviation() | |
| 45 | { | ||
| 46 | LOG_TRACE("{}: called", nameId()); | ||
| 47 | 228 | } | |
| 48 | |||
| 49 | 228 | std::string NAV::AllanDeviation::typeStatic() | |
| 50 | { | ||
| 51 |
1/2✓ Branch 1 taken 228 times.
✗ Branch 2 not taken.
|
456 | return "AllanDeviation"; |
| 52 | } | ||
| 53 | |||
| 54 | ✗ | std::string NAV::AllanDeviation::type() const | |
| 55 | { | ||
| 56 | ✗ | return typeStatic(); | |
| 57 | } | ||
| 58 | |||
| 59 | 114 | std::string NAV::AllanDeviation::category() | |
| 60 | { | ||
| 61 |
1/2✓ Branch 1 taken 114 times.
✗ Branch 2 not taken.
|
228 | return "Data Processor"; |
| 62 | } | ||
| 63 | |||
| 64 | ✗ | void NAV::AllanDeviation::guiConfig() | |
| 65 | { | ||
| 66 | ✗ | if (ImGui::Checkbox("Compute Allan Deviation last", &_updateLast)) { flow::ApplyChanges(); } | |
| 67 | ✗ | if (ImGui::Checkbox("Display Confidence Intervals", &_displayConfidence)) { flow::ApplyChanges(); } | |
| 68 | ✗ | if (_displayConfidence) | |
| 69 | { | ||
| 70 | ✗ | ImGui::SameLine(); | |
| 71 | ✗ | ImGui::SetNextItemWidth(200.0F); | |
| 72 | ✗ | if (ImGui::SliderFloat("Confidence Alpha Channel", &_confidenceFillAlpha, 0.0F, 1.0F, "%.2f")) | |
| 73 | { | ||
| 74 | ✗ | flow::ApplyChanges(); | |
| 75 | } | ||
| 76 | } | ||
| 77 | |||
| 78 | ✗ | const std::array<std::string, 3> legendEntries{ "x", "y", "z" }; | |
| 79 | ✗ | constexpr std::array<double, 5> slopeTicks = { -2, -1, 0, 1, 2 }; | |
| 80 | |||
| 81 | ✗ | if (ImGui::BeginTabBar("AllanDeviationTabBar", ImGuiTabBarFlags_None)) | |
| 82 | { | ||
| 83 | ✗ | for (size_t sensorIter = 0; sensorIter < SensorType_COUNT; sensorIter++) | |
| 84 | { | ||
| 85 | ✗ | auto sensorType = static_cast<SensorType>(sensorIter); | |
| 86 | ✗ | const char* sensorName = to_string(sensorType); | |
| 87 | ✗ | auto& sensor = _sensors.at(sensorIter); | |
| 88 | ✗ | std::scoped_lock<std::mutex> lg(sensor.mutex); | |
| 89 | |||
| 90 | ✗ | if (ImGui::BeginTabItem(fmt::format("{}##TabItem {}", sensorName, size_t(id)).c_str())) | |
| 91 | { | ||
| 92 | ✗ | if (ImPlot::BeginPlot(fmt::format("Allan Deviation of {}##Plot {}", sensorName, size_t(id)).c_str())) | |
| 93 | { | ||
| 94 | ✗ | ImPlot::SetupLegend(ImPlotLocation_SouthWest, ImPlotLegendFlags_None); | |
| 95 | ✗ | ImPlot::SetupAxes("τ [s]", fmt::format("σ(τ) [{}]", unitString(sensorType)).c_str(), | |
| 96 | ImPlotAxisFlags_AutoFit, ImPlotAxisFlags_AutoFit); | ||
| 97 | ✗ | ImPlot::SetupAxisScale(ImAxis_X1, ImPlotScale_Log10); | |
| 98 | ✗ | ImPlot::SetupAxisScale(ImAxis_Y1, ImPlotScale_Log10); | |
| 99 | ✗ | ImPlot::PushStyleVar(ImPlotStyleVar_FillAlpha, _confidenceFillAlpha); | |
| 100 | ✗ | for (size_t d = 0; d < 3; d++) | |
| 101 | { | ||
| 102 | ✗ | if (_displayConfidence & !_averagingTimes.empty()) | |
| 103 | { | ||
| 104 | ✗ | ImPlot::PlotShaded(legendEntries.at(d).data(), | |
| 105 | ✗ | _averagingTimes.data(), | |
| 106 | ✗ | sensor.allanDeviationConfidenceIntervals.at(d).at(0).data(), | |
| 107 | ✗ | sensor.allanDeviationConfidenceIntervals.at(d).at(1).data(), | |
| 108 | ✗ | static_cast<int>(_averagingTimes.size())); | |
| 109 | } | ||
| 110 | ✗ | ImPlot::PlotLine(legendEntries.at(d).data(), | |
| 111 | ✗ | _averagingTimes.data(), | |
| 112 | ✗ | sensor.allanDeviation.at(d).data(), | |
| 113 | ✗ | static_cast<int>(_averagingTimes.size())); | |
| 114 | } | ||
| 115 | ✗ | ImPlot::EndPlot(); | |
| 116 | } | ||
| 117 | ✗ | if (ImGui::TreeNode("Slopes")) | |
| 118 | { | ||
| 119 | ✗ | if (ImPlot::BeginPlot("Slopes of Allan Variance")) | |
| 120 | { | ||
| 121 | ✗ | ImPlot::SetupLegend(ImPlotLocation_SouthWest, ImPlotLegendFlags_None); | |
| 122 | ✗ | ImPlot::SetupAxis(ImAxis_X1, "τ [s]", ImPlotAxisFlags_AutoFit); | |
| 123 | ✗ | ImPlot::SetupAxis(ImAxis_Y1, "µ [ ]", ImPlotAxisFlags_Lock); | |
| 124 | ✗ | ImPlot::SetupAxisScale(ImAxis_X1, ImPlotScale_Log10); | |
| 125 | ✗ | ImPlot::SetupAxisLimits(ImAxis_Y1, -2.5, 2.5); | |
| 126 | ✗ | ImPlot::SetupAxisTicks(ImAxis_Y1, slopeTicks.data(), 5, nullptr, false); | |
| 127 | ✗ | for (size_t d = 0; d < 3; d++) | |
| 128 | { | ||
| 129 | ✗ | ImPlot::PlotLine(legendEntries.at(d).data(), | |
| 130 | ✗ | _averagingTimes.data(), | |
| 131 | ✗ | sensor.slope.at(d).data(), | |
| 132 | ✗ | static_cast<int>(_averagingTimes.size())); | |
| 133 | } | ||
| 134 | ✗ | ImPlot::EndPlot(); | |
| 135 | } | ||
| 136 | |||
| 137 | ✗ | ImGui::TreePop(); | |
| 138 | } | ||
| 139 | ✗ | ImGui::EndTabItem(); | |
| 140 | } | ||
| 141 | ✗ | } | |
| 142 | ✗ | ImGui::EndTabBar(); | |
| 143 | } | ||
| 144 | ✗ | } | |
| 145 | |||
| 146 | ✗ | [[nodiscard]] json NAV::AllanDeviation::save() const | |
| 147 | { | ||
| 148 | LOG_TRACE("{}: called", nameId()); | ||
| 149 | |||
| 150 | return { | ||
| 151 | ✗ | { "displayConfidence", _displayConfidence }, | |
| 152 | ✗ | { "confidenceFillAlpha", _confidenceFillAlpha }, | |
| 153 | ✗ | { "updateLast", _updateLast }, | |
| 154 | ✗ | }; | |
| 155 | ✗ | } | |
| 156 | |||
| 157 | ✗ | void NAV::AllanDeviation::restore(json const& j) | |
| 158 | { | ||
| 159 | LOG_TRACE("{}: called", nameId()); | ||
| 160 | |||
| 161 | ✗ | if (j.contains("displayConfidence")) { j.at("displayConfidence").get_to(_displayConfidence); } | |
| 162 | ✗ | if (j.contains("confidenceFillAlpha")) { j.at("confidenceFillAlpha").get_to(_confidenceFillAlpha); } | |
| 163 | ✗ | if (j.contains("updateLast")) { j.at("updateLast").get_to(_updateLast); } | |
| 164 | ✗ | } | |
| 165 | |||
| 166 | ✗ | bool NAV::AllanDeviation::initialize() | |
| 167 | { | ||
| 168 | LOG_TRACE("{}: called", nameId()); | ||
| 169 | |||
| 170 | ✗ | _startingInsTime.reset(); | |
| 171 | |||
| 172 | ✗ | for (auto& sensor : _sensors) | |
| 173 | { | ||
| 174 | ✗ | std::scoped_lock<std::mutex> lg(sensor.mutex); | |
| 175 | ✗ | sensor.cumSum = { Eigen::Vector3d::Zero() }; | |
| 176 | ✗ | sensor.allanSum = {}; | |
| 177 | ✗ | sensor.allanVariance = {}; | |
| 178 | ✗ | sensor.allanDeviation = {}; | |
| 179 | ✗ | sensor.slope = {}; | |
| 180 | ✗ | sensor.allanDeviationConfidenceIntervals = {}; | |
| 181 | ✗ | } | |
| 182 | |||
| 183 | ✗ | _averagingFactors.clear(); | |
| 184 | ✗ | _averagingTimes.clear(); | |
| 185 | ✗ | _observationCount.clear(); | |
| 186 | |||
| 187 | ✗ | _imuObsCount = 0; | |
| 188 | |||
| 189 | ✗ | _nextAveragingFactorExponent = 1; | |
| 190 | ✗ | _nextAveragingFactor = 1; | |
| 191 | ✗ | _samplingInterval = 0.0; | |
| 192 | |||
| 193 | ✗ | return true; | |
| 194 | } | ||
| 195 | |||
| 196 | ✗ | void NAV::AllanDeviation::deinitialize() | |
| 197 | { | ||
| 198 | LOG_TRACE("{}: called", nameId()); | ||
| 199 | ✗ | } | |
| 200 | |||
| 201 | ✗ | void NAV::AllanDeviation::receiveImuObs(NAV::InputPin::NodeDataQueue& queue, size_t /* pinIdx */) | |
| 202 | { | ||
| 203 | ✗ | auto obs = std::static_pointer_cast<const ImuObs>(queue.extract_front()); | |
| 204 | |||
| 205 | ✗ | bool lastMessage = false; | |
| 206 | ✗ | if (queue.empty() | |
| 207 | ✗ | && inputPins[INPUT_PORT_INDEX_IMU_OBS].isPinLinked() | |
| 208 | ✗ | && inputPins[INPUT_PORT_INDEX_IMU_OBS].link.getConnectedPin()->noMoreDataAvailable) | |
| 209 | { | ||
| 210 | ✗ | lastMessage = true; | |
| 211 | } | ||
| 212 | |||
| 213 | // save InsTime of first imuObs for sampling interval computation | ||
| 214 | ✗ | if (_startingInsTime.empty()) { _startingInsTime = obs->insTime; } | |
| 215 | |||
| 216 | ✗ | _imuObsCount++; | |
| 217 | { | ||
| 218 | ✗ | std::scoped_lock<std::mutex> lg(_sensors.at(Accel).mutex); | |
| 219 | ✗ | _sensors.at(Accel).cumSum.emplace_back(_sensors.at(Accel).cumSum.back() + obs->p_acceleration); | |
| 220 | ✗ | } | |
| 221 | { | ||
| 222 | ✗ | std::scoped_lock<std::mutex> lg(_sensors.at(Gyro).mutex); | |
| 223 | ✗ | _sensors.at(Gyro).cumSum.emplace_back(_sensors.at(Gyro).cumSum.back() + obs->p_angularRate); | |
| 224 | ✗ | } | |
| 225 | |||
| 226 | // extending _averagingFactors if necessary | ||
| 227 | ✗ | if (_imuObsCount == _nextAveragingFactor * 2) | |
| 228 | { | ||
| 229 | ✗ | _averagingFactors.push_back(_nextAveragingFactor); | |
| 230 | ✗ | _observationCount.push_back(0); | |
| 231 | |||
| 232 | ✗ | for (auto& sensor : _sensors) | |
| 233 | { | ||
| 234 | ✗ | for (size_t d = 0; d < 3; d++) | |
| 235 | { | ||
| 236 | ✗ | sensor.allanSum.at(d).push_back(0); | |
| 237 | } | ||
| 238 | } | ||
| 239 | |||
| 240 | // computation of next averaging factor | ||
| 241 | ✗ | while (static_cast<unsigned int>(round(pow(10., static_cast<double>(_nextAveragingFactorExponent) / _averagingFactorsPerDecade))) == _nextAveragingFactor) | |
| 242 | { | ||
| 243 | ✗ | _nextAveragingFactorExponent++; | |
| 244 | } | ||
| 245 | |||
| 246 | ✗ | _nextAveragingFactor = static_cast<unsigned int>(round(pow(10., static_cast<double>(_nextAveragingFactorExponent) / _averagingFactorsPerDecade))); | |
| 247 | } | ||
| 248 | |||
| 249 | // computation Allan sum | ||
| 250 | ✗ | for (size_t k = 0; k < _averagingFactors.size(); k++) | |
| 251 | { | ||
| 252 | ✗ | for (auto& sensor : _sensors) | |
| 253 | { | ||
| 254 | ✗ | std::scoped_lock<std::mutex> lg(sensor.mutex); | |
| 255 | ✗ | Eigen::Vector3d tempSum = sensor.cumSum.at(_imuObsCount) | |
| 256 | ✗ | - 2 * sensor.cumSum.at(_imuObsCount - static_cast<unsigned int>(_averagingFactors.at(k))) | |
| 257 | ✗ | + sensor.cumSum.at(_imuObsCount - 2 * static_cast<unsigned int>(_averagingFactors.at(k))); | |
| 258 | |||
| 259 | ✗ | for (size_t d = 0; d < 3; d++) | |
| 260 | { | ||
| 261 | ✗ | sensor.allanSum.at(d).at(k) += pow(tempSum(static_cast<Eigen::Index>(d)), 2); | |
| 262 | } | ||
| 263 | ✗ | } | |
| 264 | ✗ | _observationCount.at(k)++; | |
| 265 | } | ||
| 266 | |||
| 267 | // computation of Allan Variance and Deviation | ||
| 268 | ✗ | if (!_updateLast || lastMessage) | |
| 269 | { | ||
| 270 | ✗ | _samplingInterval = static_cast<double>((obs->insTime - _startingInsTime).count()) / _imuObsCount; | |
| 271 | ✗ | _averagingTimes.resize(_averagingFactors.size(), 0.); | |
| 272 | ✗ | _confidenceMultiplicationFactor.resize(_averagingFactors.size(), 0.); | |
| 273 | ✗ | for (size_t k = 0; k < _averagingFactors.size(); k++) | |
| 274 | { | ||
| 275 | ✗ | _averagingTimes.at(k) = _averagingFactors.at(k) * _samplingInterval; | |
| 276 | ✗ | _confidenceMultiplicationFactor.at(k) = sqrt(0.5 / (_imuObsCount / _averagingFactors.at(k) - 1)); | |
| 277 | } | ||
| 278 | |||
| 279 | ✗ | for (auto& sensor : _sensors) | |
| 280 | { | ||
| 281 | ✗ | std::scoped_lock<std::mutex> lg(sensor.mutex); | |
| 282 | ✗ | for (size_t d = 0; d < 3; d++) | |
| 283 | { | ||
| 284 | ✗ | sensor.allanVariance.at(d).resize(_averagingFactors.size(), 0.); | |
| 285 | ✗ | sensor.allanDeviation.at(d).resize(_averagingFactors.size(), 0.); | |
| 286 | |||
| 287 | ✗ | for (size_t j = 0; j < 2; j++) | |
| 288 | { | ||
| 289 | ✗ | sensor.allanDeviationConfidenceIntervals.at(d).at(j).resize(_averagingFactors.size(), 0.); | |
| 290 | } | ||
| 291 | |||
| 292 | ✗ | for (size_t k = 0; k < _averagingFactors.size(); k++) | |
| 293 | { | ||
| 294 | ✗ | sensor.allanVariance.at(d).at(k) = sensor.allanSum.at(d).at(k) / (2 * pow(_averagingFactors.at(k), 2) * _observationCount.at(k)); | |
| 295 | |||
| 296 | ✗ | sensor.allanDeviation.at(d).at(k) = sqrt(sensor.allanVariance.at(d).at(k)); | |
| 297 | |||
| 298 | ✗ | sensor.allanDeviationConfidenceIntervals.at(d).at(0).at(k) = sensor.allanDeviation.at(d).at(k) * (1 - _confidenceMultiplicationFactor.at(k)); | |
| 299 | ✗ | sensor.allanDeviationConfidenceIntervals.at(d).at(1).at(k) = sensor.allanDeviation.at(d).at(k) * (1 + _confidenceMultiplicationFactor.at(k)); | |
| 300 | } | ||
| 301 | } | ||
| 302 | |||
| 303 | // Compute Slopes | ||
| 304 | ✗ | for (size_t d = 0; d < 3; d++) { sensor.slope.at(d).resize(_averagingFactors.size(), 0.); } | |
| 305 | |||
| 306 | ✗ | for (size_t k = 0; k < _averagingFactors.size(); k++) | |
| 307 | { | ||
| 308 | ✗ | size_t lo = (k == 0 ? k : k - 1); | |
| 309 | ✗ | size_t hi = (k == _averagingFactors.size() - 1 ? k : k + 1); | |
| 310 | |||
| 311 | ✗ | double divisor = log(_averagingFactors.at(hi) / _averagingFactors.at(lo)); | |
| 312 | |||
| 313 | ✗ | for (size_t d = 0; d < 3; d++) | |
| 314 | { | ||
| 315 | ✗ | sensor.slope.at(d).at(k) = log(sensor.allanVariance.at(d).at(hi) / sensor.allanVariance.at(d).at(lo)) / divisor; | |
| 316 | } | ||
| 317 | } | ||
| 318 | ✗ | } | |
| 319 | } | ||
| 320 | |||
| 321 | // empty _cumSum for performance's sake | ||
| 322 | ✗ | if (lastMessage) | |
| 323 | { | ||
| 324 | ✗ | for (auto& sensor : _sensors) { sensor.cumSum.clear(); } | |
| 325 | } | ||
| 326 | ✗ | } | |
| 327 | |||
| 328 | ✗ | const char* NAV::AllanDeviation::to_string(SensorType sensorType) | |
| 329 | { | ||
| 330 | ✗ | switch (sensorType) | |
| 331 | { | ||
| 332 | ✗ | case SensorType::Accel: | |
| 333 | ✗ | return "Accelerometer"; | |
| 334 | ✗ | case SensorType::Gyro: | |
| 335 | ✗ | return "Gyroscope"; | |
| 336 | ✗ | case SensorType::SensorType_COUNT: | |
| 337 | ✗ | return "Error: Count"; | |
| 338 | } | ||
| 339 | ✗ | return ""; | |
| 340 | } | ||
| 341 | |||
| 342 | ✗ | const char* NAV::AllanDeviation::unitString(SensorType sensorType) | |
| 343 | { | ||
| 344 | ✗ | switch (sensorType) | |
| 345 | { | ||
| 346 | ✗ | case SensorType::Accel: | |
| 347 | ✗ | return "m/s²"; | |
| 348 | ✗ | case SensorType::Gyro: | |
| 349 | ✗ | return "rad/s"; | |
| 350 | ✗ | case SensorType::SensorType_COUNT: | |
| 351 | ✗ | return "Error: Count"; | |
| 352 | } | ||
| 353 | ✗ | return ""; | |
| 354 | } | ||
| 355 |