您知道 Matlab interp1 函数的任何 C 实现（只是“线性”函数）吗？我知道一个用于Java的。
问问题
22327 次
9 回答
13
我已将 Luis 的代码移植到 C++ 中。它似乎有效，但我没有检查过很多，所以请注意并重新检查你的结果。
#include <vector>
#include <cfloat>
#include <math.h>
vector< float > interp1( vector< float > &x, vector< float > &y, vector< float > &x_new )
{
vector< float > y_new;
y_new.reserve( x_new.size() );
std::vector< float > dx, dy, slope, intercept;
dx.reserve( x.size() );
dy.reserve( x.size() );
slope.reserve( x.size() );
intercept.reserve( x.size() );
for( int i = 0; i < x.size(); ++i ){
if( i < x.size()-1 )
{
dx.push_back( x[i+1] - x[i] );
dy.push_back( y[i+1] - y[i] );
slope.push_back( dy[i] / dx[i] );
intercept.push_back( y[i] - x[i] * slope[i] );
}
else
{
dx.push_back( dx[i-1] );
dy.push_back( dy[i-1] );
slope.push_back( slope[i-1] );
intercept.push_back( intercept[i-1] );
}
}
for ( int i = 0; i < x_new.size(); ++i )
{
int idx = findNearestNeighbourIndex( x_new[i], x );
y_new.push_back( slope[idx] * x_new[i] + intercept[idx] );
}
}
int findNearestNeighbourIndex( float value, vector< float > &x )
{
float dist = FLT_MAX;
int idx = -1;
for ( int i = 0; i < x.size(); ++i ) {
float newDist = value - x[i];
if ( newDist > 0 && newDist < dist ) {
dist = newDist;
idx = i;
}
}
return idx;
}
于 2013-02-04T16:38:49.283 回答
8
我自己实现了这个线性插值（有些是用西班牙语写的，抱歉）。名为encuentraValorMasProximo的函数只是在数组 (xD) 中找到与另一个 (xx[i]) 最接近的值 (elementoMasProximo) 和索引 (indiceEnVector)。
void interp1(int *x, int x_tam, double *y, int *xx, int xx_tam, double *yy)
{
double *dx, *dy, *slope, *intercept, *elementoMasProximo, *xD;
int i, *indiceEnVector;
dx=(double *)calloc(x_tam-1,sizeof(double));
dy=(double *)calloc(x_tam-1,sizeof(double));
slope=(double *)calloc(x_tam-1,sizeof(double));
intercept=(double *)calloc(x_tam-1,sizeof(double));
indiceEnVector=(int *) malloc(sizeof(int));
elementoMasProximo=(double *) malloc(sizeof(double));
xD=(double *)calloc(x_tam,sizeof(double));
for(i=0;i<x_tam;i++){
xD[i]=x[i];
}
for(i = 0; i < x_tam; i++){
if(i<x_tam-1){
dx[i] = x[i + 1] - x[i];
dy[i] = y[i + 1] - y[i];
slope[i] = dy[i] / dx[i];
intercept[i] = y[i] - x[i] * slope[i];
}else{
dx[i]=dx[i-1];
dy[i]=dy[i-1];
slope[i]=slope[i-1];
intercept[i]=intercept[i-1];
}
}
for (i = 0; i < xx_tam; i++) {
encuentraValorMasProximo(xx[i], xD, x_tam, x_tam, elementoMasProximo, indiceEnVector);
yy[i]=slope[*indiceEnVector] * xx[i] + intercept[*indiceEnVector];
}
}
测试功能可以是：
void main(){
int x_tam, xx_tam, i;
double *yy;
int x[]={3,6,9};
double y[]={6,12,18};
int xx[]={1,2,3,4,5,6,7,8,9,10};
x_tam=3;
xx_tam=10;
yy=(double *) calloc(xx_tam,sizeof(double));
interp1(x, x_tam, y, xx, xx_tam, yy);
for(i=0;i<xx_tam;i++){
printf("%d\t%f\n",xx[i],yy[i]);
}
}
及其结果：
1 2.000000
2 4.000000
3 6.000000
4 8.000000
5 10.000000
6 12.000000
7 14.000000
8 16.000000
9 18.000000
10 20.000000
于 2012-02-22T21:52:17.297 回答
6
常用函数的出色实现可以在《C 中的数值食谱》一书中找到，该书可在线免费查看。第 3.1.2 章有一个线性插值方法，本章的其余部分介绍更高级的方法。
我强烈推荐这本书，它写得很好，涵盖了大量的算法，这些算法也以一种非常有效且仍然可以理解的方式实现。
于 2012-02-22T13:39:30.773 回答
5
Luis 提交的 C 代码存在问题。首先缺少 encuentraValorMasProximo。第二，数组预留少一个数字。我也清理了这个功能。下面是具有 encuentraValorMasProximo 函数（重命名为 findNearestNeighbourIndex）的功能 C 代码。
#include <float.h>
int findNearestNeighbourIndex( double value, double *x, int len )
{
double dist;
int idx;
int i;
idx = -1;
dist = DBL_MAX;
for ( i = 0; i < len; i++ ) {
double newDist = value - x[i];
if ( newDist > 0 && newDist < dist ) {
dist = newDist;
idx = i;
}
}
return idx;
}
void interp1(double *x, int x_tam, double *y, double *xx, int xx_tam, double *yy)
{
double dx, dy, *slope, *intercept;
int i, indiceEnVector;
slope=(double *)calloc(x_tam,sizeof(double));
intercept=(double *)calloc(x_tam,sizeof(double));
for(i = 0; i < x_tam; i++){
if(i<x_tam-1){
dx = x[i + 1] - x[i];
dy = y[i + 1] - y[i];
slope[i] = dy / dx;
intercept[i] = y[i] - x[i] * slope[i];
}else{
slope[i]=slope[i-1];
intercept[i]=intercept[i-1];
}
}
for (i = 0; i < xx_tam; i++) {
indiceEnVector = findNearestNeighbourIndex( xx[i], x, x_tam);
if (indiceEnVector != -1)
yy[i]=slope[indiceEnVector] * xx[i] + intercept[indiceEnVector];
else
yy[i]=DBL_MAX;
}
free(slope);
free(intercept);
}
于 2015-01-15T20:30:00.293 回答
3
@user1097111，您的代码存在一个错误，在函数 findNearestNeighbourIndex 中，它应该是 if ( newDist >= 0 && newDist < dist )，而不是 if ( newDist > 0 && newDist < dist )。
于 2015-06-18T09:32:00.020 回答
2
您可以查看 GSL（数值科学图书馆）。有很多类似 Matlab 的函数，其中还有一维插值。
我现在在手机上，抱歉，无法提供链接。
于 2012-02-22T12:51:18.680 回答
2
你知道Matlab 编码器吗？它自动从 Matlab 代码生成 c/c++ 代码。如果你把它作为你的 Matlab 包的一部分，你可以通过它运行 interp1 函数，看看 Matlab 吐出什么。
于 2012-02-22T13:28:28.243 回答
2
如果这对将来的某人有所帮助，这是一个没有临时数组且没有 0 错误的版本。
#include <iostream>
#include <vector>
#include <limits>
#include <cmath>
template<typename Real>
int nearestNeighbourIndex(std::vector<Real> &x, Real &value)
{
Real dist = std::numeric_limits<Real>::max();
Real 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 Real>
std::vector<Real> interp1(std::vector<Real> &x, std::vector<Real> &y, std::vector<Real> &x_new)
{
std::vector<Real> y_new;
Real 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{1, 2, 3, 4, 5};
vector<float> y{5, 6, 7, 8, 9};
vector<float> newx{0, 5, 6.123, 12.123, 2, 4};
auto res = interp1(x, y, newx);
for (auto i: res)
cout << i << " ";
cout << endl;
}
于 2016-09-23T18:22:52.413 回答
1
请参阅文件交换中的 lininterp1f。
于 2013-10-22T23:49:15.587 回答