Matlab interp1函数的C实现（线性插值）

Question

您知道任何C语言实现Matlab interp1函数（仅“linear”函数）的方法吗？ 我知道有一个适用于Java。

- Luis Andrés García

---

为什么这个问题被关闭了？我发现最佳答案（纠正索引错误的答案）非常有用，即使我不是在寻找C++实现。我可以轻松地将其翻译成其他语言。 - Frank Hileman

Answer 1

我已经将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;
}

- Nikita Rokotyan

---

他特别要求使用C语言实现，因此这并不是一个合适的答案。 - antonijn

这很好，因为我正在寻找一个C++实现。谢谢！ - rayryeng

Answer 2

我已经自己实现了这个线性插值（其中一部分是用西班牙语编写的，抱歉）。名为encuentraValorMasProximo的函数只是在数组(xD)中寻找最近的值(elementoMasProximo)和索引(indiceEnVector)，以匹配另一个值(xx[i])。

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];
}
}

test函数可能是这样的：

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

- Luis Andrés García

Answer 3

优秀的常用函数实现可以在书籍《Numerical Recipes in C》中找到，该书可在线免费阅读。第3.1.2章节讲解了线性插值方法，其余章节涵盖了更高级的内容。

我强烈推荐这本书，它写得非常好，并涵盖了大量算法，这些算法的实现方式也非常高效而易于理解。

- jupp0r

Answer 4

路易斯提交的C代码存在问题。首先，找不到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);
}

- user1097111

Answer 5

如果有人在未来需要帮助，这是一个没有临时数组和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;
}

- canesin

---

非常有帮助且易于翻译成其他语言。 - Frank Hileman

Answer 6

@user1097111，你的代码存在一个错误，在函数findNearestNeighbourIndex中应该是 if ( newDist >= 0 && newDist < dist )，而不是 if ( newDist > 0 && newDist < dist )。

- Haiwei Liu

Answer 7

您可以查看 GSL（数值科学库）。其中有许多类似于Matlab的函数，包括一维插值。

很抱歉，我现在在手机上，无法提供链接。

- Ivan Angelov

---

因为GSL插值，所以给它加上+1。然而，仅仅为了一个函数而使用它太过沉重，但作为替代方案还是可以的。 - Luis Andrés García

Answer 8

你是否了解Matlab Coder？它可以自动生成c/c++代码，从Matlab代码中生成。如果你的Matlab包中有这个功能，你只需要通过它运行interp1函数，看看Matlab会输出什么。

- prototoast

---

我已经尝试过了，但是它不可读。 - Luis Andrés García

Answer 9

在fileexchange中查看lininterp1f。

- Felipe G. Nievinski