一维线性插值算法C++详细实现

算法概述：

算法实现：

下面提供一套健壮的一维线性插值 C++ 实现，重点考虑了数据预处理、边界处理、重复点合并及浮点容差等鲁棒性细节。

#include <vector> #include <algorithm> #include <stdexcept> #include <cmath> #include <limits> #include <utility> /** * 鲁棒的一维线性插值器 * @tparam T 浮点类型（float/double/long double） */ template<typename T> class LinearInterpolator { public: using value_type = T; /** * 构造函数：接收原始数据点，自动进行排序、去重（相同 x 取 y 均值） * @param x 自变量数组 * @param y 因变量数组 * @throw std::invalid_argument 如果 x 和 y 长度不相等 */ LinearInterpolator(const std::vector<T>& x, const std::vector<T>& y) { if (x.size() != y.size()) { throw std::invalid_argument("x and y must have the same size"); } if (x.empty()) { return; // 空数据，后续插值返回 NaN } // 1. 合并为点对并排序 std::vector<std::pair<T, T>> points; points.reserve(x.size()); for (size_t i = 0; i < x.size(); ++i) { points.emplace_back(x[i], y[i]); } std::sort(points.begin(), points.end(), [](const auto& a, const auto& b) { return a.first < b.first; }); // 2. 合并相同 x 的点（考虑浮点容差） const T eps = std::numeric_limits<T>::epsilon() * 100; std::vector<std::pair<T, T>> merged; merged.reserve(points.size()); for (const auto& p : points) { if (merged.empty() || std::abs(p.first - merged.back().first) > eps * (std::abs(p.first) + T(1))) { merged.push_back(p); } else { // 对相同 x 的 y 取平均值（可根据需求改为其他策略） merged.back().second = (merged.back().second + p.second) / T(2); } } // 3. 存储处理后的数据 x_.reserve(merged.size()); y_.reserve(merged.size()); for (auto& p : merged) { x_.push_back(p.first); y_.push_back(p.second); } } /** * 单点插值（外推时返回最近点的 y 值，即限幅 clamp） * @param x_query 待插值的自变量 * @return 插值结果（无数据点时返回 NaN） */ T interpolate(T x_query) const { if (x_.empty()) { return std::numeric_limits<T>::quiet_NaN(); } // 处理边界：直接返回最近点的 y 值（限幅） if (x_query <= x_.front()) { return y_.front(); } if (x_query >= x_.back()) { return y_.back(); } // 二分查找第一个大于 x_query 的索引 auto it = std::upper_bound(x_.begin(), x_.end(), x_query); size_t right_idx = it - x_.begin(); // 右侧点索引 size_t left_idx = right_idx - 1; // 左侧点索引 T x_left = x_[left_idx]; T x_right = x_[right_idx]; T y_left = y_[left_idx]; T y_right = y_[right_idx]; // 线性插值（区间宽度已保证 > 0，因为重复 x 已合并） T t = (x_query - x_left) / (x_right - x_left); return y_left + t * (y_right - y_left); } /** * 多点插值 * @param x_queries 待插值的自变量数组 * @return 对应的插值结果数组 */ std::vector<T> interpolate(const std::vector<T>& x_queries) const { std::vector<T> results; results.reserve(x_queries.size()); for (T xq : x_queries) { results.push_back(interpolate(xq)); } return results; } private: std::vector<T> x_; ///< 排序且去重后的自变量 std::vector<T> y_; ///< 对应的因变量 };

鲁棒性说明：

1. 输入校验

   • 检查x和y长度是否一致，不一致则抛异常。

   • 空数据集时插值返回NaN，避免未定义行为。

2. 数据预处理

   • 自动排序：无论原始数据顺序如何，均按x升序排列，保证后续二分查找正确。

   • 合并重复x：使用相对容差判断x是否相等（考虑浮点精度），对相同x的y取平均值，避免除零或歧义。

   • 经过预处理后，内部存储的x严格递增且无重复，插值区间宽度恒为正。

3. 边界处理

   • 查询点超出数据范围时采用限幅（clamp）策略，直接返回最近数据点的y值，避免不合理外推导致结果失控。

   • 若数据仅有一个点，任何查询均返回该点的y。

4. 数值稳定性

   • 合并重复点时使用相对容差，避免因浮点误差导致微小区间被保留，从而防止后续除零风险。

   • 插值计算t = (x_query - x_left) / (x_right - x_left)在区间宽度 > 0 时安全。

5. 通用性

   • 模板支持float、double、long double等浮点类型。

   • 可轻松扩展外推模式（如线性外推），只需修改边界处理部分。

使用示例：

#include <iostream> #include <vector> int main() { // 原始数据（包含重复 x 和乱序） std::vector<double> x = {3.0, 1.0, 2.0, 2.0, 4.0}; std::vector<double> y = {30.0, 10.0, 20.0, 25.0, 40.0}; // 构造插值器（自动排序并合并 x=2.0 的点，取 y 均值 (20+25)/2 = 22.5） LinearInterpolator<double> interp(x, y); // 插值 std::vector<double> queries = {0.0, 1.5, 2.0, 3.5, 5.0}; auto results = interp.interpolate(queries); // 输出结果 for (size_t i = 0; i < queries.size(); ++i) { std::cout << "f(" << queries[i] << ") = " << results[i] << '\n'; } return 0; } 输出： f(0) = 10 // 限幅到最近点 x=1, y=10 f(1.5) = 15 // 插值于 (1,10) 和 (2,22.5) f(2) = 22.5 // 合并后的点 f(3.5) = 35 // 插值于 (3,30) 和 (4,40) f(5) = 40 // 限幅到 x=4, y=40

单点插值 int main() { // 原始数据（包含重复 x 和乱序） std::vector<double> x = {0.02, 0.04, 0.06, 0.08, 0.10}; std::vector<double> y = {13.0, 10.0, 8.0, 6.0, 6.0}; // 构造插值器（自动排序并合并 x=2.0 的点，取 y 均值 (20+25)/2 = 22.5） LinearInterpolator<double> interp(x, y); // 插值 auto results = interp.interpolate(0.075); std::cout<< "results = " << results << std::endl; return 0; }

该实现可在大多数工程场景下稳定运行，用户也可根据具体需求调整合并策略（如取第一个值）或外推方式。

语法解读：

std::sort(points.begin(), points.end(), [](const auto& a, const auto& b) { return a.first < b.first; });

1. ‌代码功能解析‌

• ‌points.begin(),points.end()‌：指定排序范围（整个容器）。
• ‌Lambda表达式‌：[](const auto& a, const auto& b) { return a.first < b.first; }
  定义比较规则：比较两个元素的first成员，若a.first < b.first，则a应排在b前面（升序）。

---

2. ‌关键点‌

• ‌适用容器‌：points需支持随机访问迭代器（如std::vector、std::deque等）。
• ‌排序稳定性‌：std::sort不是稳定排序（相等元素的相对顺序可能改变）。若需稳定排序，改用std::stable_sort。
• ‌性能‌：时间复杂度为O(n log n)。

3. ‌常见问题‌

• ‌降序排序‌：将Lambda中的<改为>即可。

const T eps = std::numeric_limits<T>::epsilon() * 100;

‌1. 核心功能‌

• ‌std::numeric_limits<T>::epsilon()‌
  返回类型T的机器精度（即 1 与大于 1 的最小可表示值之间的差值）。例如：
  • float：约1.19209e-07
  • double：约2.22045e-16
• ‌乘以 100‌：
  将机器精度放大 100 倍，通常用于避免浮点数比较时的精度误差（例如比较两个浮点数是否“近似相等”）。

‌2. 典型应用场景‌

• ‌浮点数比较‌：
  用于判断两个浮点数是否在允许的误差范围内相等：

  template<typename T> bool isApproxEqual(T a, T b, T eps) { return std::abs(a - b) <= eps; }

auto it = std::upper_bound(x_.begin(), x_.end(), x_value);

这段代码使用std::upper_bound在有序范围[x_.begin(), x_.end())中查找第一个‌大于‌x_value的元素，并返回指向该元素的迭代器it。以下是详细解析：

‌1. 核心功能‌

• ‌std::upper_bound‌：
  二分查找算法，要求输入范围‌有序‌（默认按升序），返回第一个满足*it > x_value的迭代器。
• ‌返回值it‌：
  • 若找到，it指向第一个大于x_value的元素。
  • 若所有元素均不大于x_value，则返回x_.end()。

‌自定义比较规则‌：
若元素为自定义类型或需降序排序，需传入比较函数：

auto it = std::upper_bound(x_.begin(), x_.end(), x_value, [](const auto& a, const auto& b) { return a > b; }); // 降序规则