【无标题】c++ 实现 interpolate.interp1d插值

提示：文章写完后，目录可以自动生成，如何生成可参考右边的帮助文档

---

使用c++实现 interpolate.interp1d插值

interpolate.interp1d 就是已知离散点拟合出拟合出平滑的曲线（有不同的方法拟合，线性，二次）比如一致x0 函数值y0,x2函数值为y2,x1对应的函数值为多少,；x0<x1<x2;就是求出（x0，y0）(x2,y2)之间的直线，然后就可以求出x1对应的函数值。代码：

#include <iostream>
#include <vector>
#include <limits>
#include <cmath>


template<typename T>
int nearestNeighbourIndex(std::vector<T>& x, T& value)
{
	T dist = std::numeric_limits<T>::max();
	T newDist = dist;
	size_t idx = 0;

	for (size_t i = 0; i < x.size(); ++i) {
		newDist = std::abs(value - x[i]);
		if (newDist <= dist) {
			dist = newDist;
			idx = i;
		}
	}

	return idx;
}



template<typename T>
std::vector<T> interp1(std::vector<T>& x, std::vector<T>& y, std::vector<T>& x_new)
{
	std::vector<T> y_new;
	T dx, dy, m, b;
	size_t x_max_idx = x.size() - 1;
	size_t x_new_size = x_new.size();

	y_new.reserve(x_new_size);

	for (size_t i = 0; i < x_new_size; ++i)
	{
		size_t idx = nearestNeighbourIndex(x, x_new[i]);

		if (x[idx] > x_new[i])
		{
			dx = idx > 0 ? (x[idx] - x[idx - 1]) : (x[idx + 1] - x[idx]);
			dy = idx > 0 ? (y[idx] - y[idx - 1]) : (y[idx + 1] - y[idx]);
		}
		else
		{
			dx = idx < x_max_idx ? (x[idx + 1] - x[idx]) : (x[idx] - x[idx - 1]);
			dy = idx < x_max_idx ? (y[idx + 1] - y[idx]) : (y[idx] - y[idx - 1]);
		}

		m = dy / dx;
		b = y[idx] - x[idx] * m;

		y_new.push_back(x_new[i] * m + b);
	}

	return y_new;
}

int main() {
	vector<float> x{ 0, 1, 2, 4, 7 };
	vector<float> y{ 37, 37, 37, 38, 39 };
	vector<float> newx{ 0, 1, 2,3, 4, 5, 6,7,8,9,10 };

	vector<float> res = interp1(x, y, newx);
	cout << res[9.6]<<endl;//vector值访问下表如果不是整数，会自动取整
	for (float i : res)
		cout <<i << " ";
	cout << endl;
}