From d9f4d540ef334013eb404ce91b3b446e5fc917ff Mon Sep 17 00:00:00 2001
From: Martin Braun <martin.braun@ettus.com>
Date: Tue, 31 Mar 2020 21:38:13 -0700
Subject: uhd: math: Add interpolation.hpp
- Moves linear_interp from cal to utils
- Moves the interp_mode enum class to interpolation.hpp
- Adds three interpolation methods for maps: at_interpolate_1d(),
at_nearest(), at_lin_interp()
- Adds unit tests
---
host/include/uhd/cal/CMakeLists.txt | 1 -
host/include/uhd/cal/interpolation.hpp | 16 ---
host/include/uhd/cal/iq_cal.hpp | 4 +-
host/include/uhd/utils/CMakeLists.txt | 1 +
host/include/uhd/utils/interpolation.hpp | 16 +++
host/include/uhd/utils/math.hpp | 25 +----
host/lib/cal/cal_python.hpp | 8 +-
host/lib/cal/iq_cal.cpp | 2 +
host/lib/include/uhdlib/utils/interpolation.hpp | 140 ++++++++++++++++++++++++
host/tests/CMakeLists.txt | 1 +
host/tests/cal_data_iq_test.cpp | 1 +
host/tests/interpolation_test.cpp | 35 ++++++
host/tests/math_test.cpp | 10 --
13 files changed, 208 insertions(+), 52 deletions(-)
delete mode 100644 host/include/uhd/cal/interpolation.hpp
create mode 100644 host/include/uhd/utils/interpolation.hpp
create mode 100644 host/lib/include/uhdlib/utils/interpolation.hpp
create mode 100644 host/tests/interpolation_test.cpp
diff --git a/host/include/uhd/cal/CMakeLists.txt b/host/include/uhd/cal/CMakeLists.txt
index 2aba1a91c..8d8dfa11e 100644
--- a/host/include/uhd/cal/CMakeLists.txt
+++ b/host/include/uhd/cal/CMakeLists.txt
@@ -7,7 +7,6 @@
UHD_INSTALL(FILES
database.hpp
container.hpp
- interpolation.hpp
iq_cal.hpp
iq_cal_generated.h
DESTINATION ${INCLUDE_DIR}/uhd/cal
diff --git a/host/include/uhd/cal/interpolation.hpp b/host/include/uhd/cal/interpolation.hpp
deleted file mode 100644
index 5cdbcfa54..000000000
--- a/host/include/uhd/cal/interpolation.hpp
+++ /dev/null
@@ -1,16 +0,0 @@
-//
-// Copyright 2020 Ettus Research, a National Instruments Brand
-//
-// SPDX-License-Identifier: GPL-3.0-or-later
-//
-
-#ifndef INCLUDED_LIBUHD_CAL_INTERP_HPP
-#define INCLUDED_LIBUHD_CAL_INTERP_HPP
-
-namespace uhd { namespace usrp { namespace cal {
-
-enum class interp_mode { NEAREST_NEIGHBOR, LINEAR };
-
-}}} // namespace uhd::usrp::cal
-
-#endif /* INCLUDED_LIBUHD_CAL_INTERP_HPP */
diff --git a/host/include/uhd/cal/iq_cal.hpp b/host/include/uhd/cal/iq_cal.hpp
index ded082698..78f1c9177 100644
--- a/host/include/uhd/cal/iq_cal.hpp
+++ b/host/include/uhd/cal/iq_cal.hpp
@@ -8,8 +8,8 @@
#define INCLUDED_LIBUHD_CAL_IQ_DATA_HPP
#include <uhd/cal/container.hpp>
-#include <uhd/cal/interpolation.hpp>
#include <uhd/config.hpp>
+#include <uhd/utils/interpolation.hpp>
#include <complex>
#include <memory>
#include <string>
@@ -35,7 +35,7 @@ public:
// \param interp The new interpolation mode
// \throws uhd::value_error if the given interpolation mode is not
// supported.
- virtual void set_interp_mode(const interp_mode interp) = 0;
+ virtual void set_interp_mode(const uhd::math::interp_mode interp) = 0;
//! Return a calibration coefficient for a given frequency
//
diff --git a/host/include/uhd/utils/CMakeLists.txt b/host/include/uhd/utils/CMakeLists.txt
index b70de86ba..17d5b2380 100644
--- a/host/include/uhd/utils/CMakeLists.txt
+++ b/host/include/uhd/utils/CMakeLists.txt
@@ -18,6 +18,7 @@ UHD_INSTALL(FILES
fp_compare_epsilon.ipp
gain_group.hpp
graph_utils.hpp
+ interpolation.hpp
log.hpp
log_add.hpp
math.hpp
diff --git a/host/include/uhd/utils/interpolation.hpp b/host/include/uhd/utils/interpolation.hpp
new file mode 100644
index 000000000..42a7d8fb1
--- /dev/null
+++ b/host/include/uhd/utils/interpolation.hpp
@@ -0,0 +1,16 @@
+//
+// Copyright 2020 Ettus Research, a National Instruments Brand
+//
+// SPDX-License-Identifier: GPL-3.0-or-later
+//
+
+#ifndef INCLUDED_UHD_INTERP_HPP
+#define INCLUDED_UHD_INTERP_HPP
+
+namespace uhd { namespace math {
+
+enum class interp_mode { NEAREST_NEIGHBOR, LINEAR };
+
+}} // namespace uhd::math
+
+#endif /* INCLUDED_UHD_INTERP_HPP */
diff --git a/host/include/uhd/utils/math.hpp b/host/include/uhd/utils/math.hpp
index bb721d3af..f5aef5ea1 100644
--- a/host/include/uhd/utils/math.hpp
+++ b/host/include/uhd/utils/math.hpp
@@ -1,4 +1,9 @@
-
+//
+// Copyright 2014-2015 Ettus Research LLC
+// Copyright 2018 Ettus Research, a National Instruments Company
+//
+// SPDX-License-Identifier: GPL-3.0-or-later
+//
#ifndef INCLUDED_UHD_UTILS_MATH_HPP
#define INCLUDED_UHD_UTILS_MATH_HPP
@@ -246,24 +251,6 @@ inline IntegerType gcd(IntegerType x, IntegerType y)
return _bmint::gcd<IntegerType>(x, y);
}
-//! Linearly interpolate f(x) given f(x0) = y0 and f(x1) = y1
-//
-// This draws a line through the coordinates x0/y0 and x1/y1, and then returns
-// the y-value for the given x-value on said line.
-//
-// \throws uhd::runtime_error if x0 == x1, since that doesn't allow us to
-// interpolate.
-template <typename InterpType>
-inline InterpType linear_interp(
- InterpType x, InterpType x0, InterpType y0, InterpType x1, InterpType y1)
-{
- if (x0 == x1) {
- throw uhd::runtime_error("linear_interp(): x0 and x1 must differ!");
- }
- return y0 + (x - x0) * (y1 - y0) / (x1 - x0);
-}
-
-
} // namespace math
} // namespace uhd
diff --git a/host/lib/cal/cal_python.hpp b/host/lib/cal/cal_python.hpp
index 39a13d94c..e8f19eef9 100644
--- a/host/lib/cal/cal_python.hpp
+++ b/host/lib/cal/cal_python.hpp
@@ -8,8 +8,8 @@
#define INCLUDED_UHD_CAL_PYTHON_HPP
#include <uhd/cal/database.hpp>
-#include <uhd/cal/interpolation.hpp>
#include <uhd/cal/iq_cal.hpp>
+#include <uhd/utils/interpolation.hpp>
std::vector<uint8_t> pybytes_to_vector(const py::bytes& data)
{
@@ -59,9 +59,9 @@ void export_cal(py::module& m)
database::write_cal_data(key, serial, pybytes_to_vector(data));
});
- py::enum_<interp_mode>(m, "interp_mode")
- .value("NEAREST_NEIGHBOR", interp_mode::NEAREST_NEIGHBOR)
- .value("LINEAR", interp_mode::LINEAR);
+ py::enum_<uhd::math::interp_mode>(m, "interp_mode")
+ .value("NEAREST_NEIGHBOR", uhd::math::interp_mode::NEAREST_NEIGHBOR)
+ .value("LINEAR", uhd::math::interp_mode::LINEAR);
py::class_<container, std::shared_ptr<container>>(m, "container")
.def("get_name", &container::get_name)
diff --git a/host/lib/cal/iq_cal.cpp b/host/lib/cal/iq_cal.cpp
index e1ed8c9cb..f5640b01e 100644
--- a/host/lib/cal/iq_cal.cpp
+++ b/host/lib/cal/iq_cal.cpp
@@ -8,10 +8,12 @@
#include <uhd/cal/iq_cal_generated.h>
#include <uhd/exception.hpp>
#include <uhd/utils/math.hpp>
+#include <uhdlib/utils/interpolation.hpp>
#include <map>
#include <string>
using namespace uhd::usrp::cal;
+using namespace uhd::math;
constexpr int VERSION_MAJOR = 1;
constexpr int VERSION_MINOR = 0;
diff --git a/host/lib/include/uhdlib/utils/interpolation.hpp b/host/lib/include/uhdlib/utils/interpolation.hpp
new file mode 100644
index 000000000..10fc06559
--- /dev/null
+++ b/host/lib/include/uhdlib/utils/interpolation.hpp
@@ -0,0 +1,140 @@
+//
+// Copyright 2020 Ettus Research, a National Instruments Brand
+//
+// SPDX-License-Identifier: GPL-3.0-or-later
+//
+
+// Various interpolation functions used within UHD
+
+#ifndef INCLUDED_UHD_UTILS_INTERP_HPP
+#define INCLUDED_UHD_UTILS_INTERP_HPP
+
+#include <uhd/utils/math.hpp>
+#include <uhd/utils/interpolation.hpp>
+#include <map>
+
+namespace uhd { namespace math {
+
+//! Linearly interpolate f(x) given f(x0) = y0 and f(x1) = y1
+//
+// This draws a line through the coordinates x0/y0 and x1/y1, and then returns
+// the y-value for the given x-value on said line.
+//
+// \throws uhd::runtime_error if x0 == x1, since that doesn't allow us to
+// interpolate.
+template <typename InterpType>
+inline InterpType linear_interp(
+ InterpType x, InterpType x0, InterpType y0, InterpType x1, InterpType y1)
+{
+ if (x0 == x1) {
+ throw uhd::runtime_error("linear_interp(): x0 and x1 must differ!");
+ }
+ return y0 + (x - x0) * (y1 - y0) / (x1 - x0);
+}
+
+/*! Like std::at, but will interpolate the data between two points
+ *
+ * Most often, this is used with a std::map (keys need to be sorted!). Unlike a
+ * regular map, this function allows looking up values in the map even if they're
+ * not a data point, by interpolating between the two data points that are closest
+ * to the given key.
+ *
+ * Example:
+ * ~~~{.cpp}
+ * std::map<double, double> data{{1.0, 2.0}, {2.0, 3.0}};
+ * std::cout << data[1.5] << std::endl; // This will fail!
+ * // This will print something depending on my_interpolation_func:
+ * std::cout << at_interpolate_1d(data, 1.5, &my_interpolation_func) << std::endl;
+ * ~~~
+ *
+ * Note: When \p key exceeds the max key of \p data, or is lower than the min key
+ * of \p data, it will return the value for the max or min key, respectively.
+ *
+ * \param data A map that stores x/y values. It must implemented lower_bound().
+ * \param key The x-value to look up
+ * \param interp_func A function that takes 5 inputs (x, x0, y0, x1, y1) and
+ * returns f(x) for f(x0) == y0 and f(x1) == y1. The inputs
+ * x, x0, and x1 must be of the key type of \p data, and the
+ * return value must be of the mapped type of \p data.
+ * \returns An interpolation of the data at key \p key. The return value type
+ * is the mapped type of \p data.
+ */
+template <typename map_type, typename interp_func_type>
+typename map_type::mapped_type at_interpolate_1d(const map_type& data,
+ const typename map_type::key_type& key,
+ interp_func_type&& interp_func)
+{
+ // Get an iterator to the next item
+ const auto next_it = data.lower_bound(key);
+ if (next_it == data.cend()) {
+ // This means key is larger than our biggest key in data, and thus we
+ // return the value of the largest key.
+ return data.crbegin()->second;
+ }
+ if (next_it == data.cbegin()) {
+ // This means freq is smaller than our smallest key, and thus we
+ // return the value of the smallest key.
+ return data.cbegin()->second;
+ }
+
+ // Get an iterator to the previous item
+ auto prev_it = next_it;
+ prev_it--;
+ const auto hi_key = next_it->first;
+ const auto hi_value = next_it->second;
+ const auto lo_key = prev_it->first;
+ const auto lo_value = prev_it->second;
+
+ return interp_func(key, lo_key, lo_value, hi_key, hi_value);
+}
+
+//! Like std::map::at, except with an approximate index
+//
+// Example:
+// ~~~{.cpp}
+// std::map<double, double> data{{1.0, 2.0}, {2.0, 3.0}};
+// std::cout << at_nearest(data, 1.72) << std::endl; // prints 3.0
+// ~~~
+//
+// This is in fact a shorthand for at_interpolate_1d(). It will look up the
+// value in \p data with the key that most closely matches \p key, i.e.,
+// at_nearest(data, key) == data[key'] if key' == argmin abs(key' - key).
+template <typename map_type>
+typename map_type::mapped_type at_nearest(
+ const map_type& data, const typename map_type::key_type& key)
+{
+ return at_interpolate_1d(
+ data,
+ key,
+ [&](const typename map_type::key_type x,
+ const typename map_type::key_type x0,
+ const typename map_type::mapped_type y0,
+ const typename map_type::key_type x1,
+ const typename map_type::mapped_type y1) ->
+ typename map_type::mapped_type { return (x1 - x < x - x0) ? y1 : y0; });
+}
+
+//! Like std::map::at, except it will linearly interpolate in one dimension
+//
+// Example:
+// ~~~{.cpp}
+// std::map<double, double> data{{1.0, 2.0}, {2.0, 3.0}};
+// std::cout << at_lin_interp(data, 1.5) << std::endl; // prints 2.5
+// ~~~
+//
+// This treats the map as a set of x/y coordinates, and returns the value from
+// the map that corresponds to a linear interpolation on those coordinates.
+//
+// For x-values greater than the maximum key, or smaller than the minimum key
+// of \p data, we return the value for the closest available key.
+template <typename map_type>
+typename map_type::mapped_type at_lin_interp(
+ const map_type& data, const typename map_type::key_type& key)
+{
+ return at_interpolate_1d(
+ data, key, &uhd::math::linear_interp<typename map_type::mapped_type>);
+}
+
+}} // namespace uhd::math
+
+#endif /* INCLUDED_UHD_UTILS_INTERP_HPP */
diff --git a/host/tests/CMakeLists.txt b/host/tests/CMakeLists.txt
index 3665a9d17..33bcdb734 100644
--- a/host/tests/CMakeLists.txt
+++ b/host/tests/CMakeLists.txt
@@ -36,6 +36,7 @@ set(test_sources
fp_compare_delta_test.cpp
fp_compare_epsilon_test.cpp
gain_group_test.cpp
+ interpolation_test.cpp
isatty_test.cpp
log_test.cpp
math_test.cpp
diff --git a/host/tests/cal_data_iq_test.cpp b/host/tests/cal_data_iq_test.cpp
index 97b35a331..523a4310e 100644
--- a/host/tests/cal_data_iq_test.cpp
+++ b/host/tests/cal_data_iq_test.cpp
@@ -10,6 +10,7 @@
#include <iostream>
using namespace uhd::usrp::cal;
+using namespace uhd::math;
BOOST_AUTO_TEST_CASE(test_iq_cal_api)
{
diff --git a/host/tests/interpolation_test.cpp b/host/tests/interpolation_test.cpp
new file mode 100644
index 000000000..b3314c34a
--- /dev/null
+++ b/host/tests/interpolation_test.cpp
@@ -0,0 +1,35 @@
+//
+// Copyright 2020 Ettus Research, a National Instruments Brand
+//
+// SPDX-License-Identifier: GPL-3.0-or-later
+//
+
+#include <uhdlib/utils/interpolation.hpp>
+#include <boost/test/unit_test.hpp>
+
+BOOST_AUTO_TEST_CASE(test_interp)
+{
+ const double x0 = 1.0, x1 = 2.0;
+ const double y0 = 2.0, y1 = 4.0;
+ const double x = 1.5;
+ const double y_exp = 3.0;
+
+ BOOST_CHECK_EQUAL(uhd::math::linear_interp<double>(x, x0, y0, x1, y1), y_exp);
+}
+
+BOOST_AUTO_TEST_CASE(test_map_interp)
+{
+ using namespace uhd::math;
+ std::map<double, double> test_data{{1.0, 1.0}, {2.0, 2.0}, {3.0, 3.0}};
+
+ BOOST_CHECK_EQUAL(at_nearest(test_data, 1.1), 1.0);
+ BOOST_CHECK_EQUAL(at_nearest(test_data, 1.9), 2.0);
+ BOOST_CHECK_EQUAL(at_nearest(test_data, 2.1), 2.0);
+ BOOST_CHECK_EQUAL(at_nearest(test_data, 1e9), 3.0);
+
+ BOOST_CHECK_CLOSE(at_lin_interp(test_data, 1.5), 1.5, 1e-6);
+ BOOST_CHECK_EQUAL(at_lin_interp(test_data, 0.1), 1.0);
+ BOOST_CHECK_EQUAL(at_lin_interp(test_data, 2.0), 2.0);
+ BOOST_CHECK_EQUAL(at_lin_interp(test_data, 137.0), 3.0);
+}
+
diff --git a/host/tests/math_test.cpp b/host/tests/math_test.cpp
index bff6fc16f..c1b6793bb 100644
--- a/host/tests/math_test.cpp
+++ b/host/tests/math_test.cpp
@@ -17,13 +17,3 @@ BOOST_AUTO_TEST_CASE(test_gcd)
{
BOOST_CHECK_EQUAL(uhd::math::gcd<int>(6, 15), 3);
}
-
-BOOST_AUTO_TEST_CASE(test_interp)
-{
- const double x0 = 1.0, x1 = 2.0;
- const double y0 = 2.0, y1 = 4.0;
- const double x = 1.5;
- const double y_exp = 3.0;
-
- BOOST_CHECK_EQUAL(uhd::math::linear_interp<double>(x, x0, y0, x1, y1), y_exp);
-}
--
cgit v1.2.3