C implementation of Matlab interp1 function (linear interpolation) [closed]

Question

Do you know any C implementation of the Matlab interp1 function (just 'linear' one)? I know one for Java.

Answer 1

I have implemented this linear interpolation myself (some of it is written in Spanish, sorry). The function called encuentraValorMasProximo just finds the nearest value (elementoMasProximo) and index (indiceEnVector) to another one (xx[i]), in an array (xD).

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

The test function could be:

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

}

And its outcome:

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

Answer 2

I've ported Luis's code to c++. It seems to be working but I haven't checked it a lot, so be aware and re-check your results.

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

Answer 3

Excelent implementations of commonly used functions can be found in the book Numerical Recipes in C, which is viewable for free online. Chapter 3.1.2 has a linear interpolation recipe, the rest of the chapter covers more advanced ones.

I can strongly recommend this book, it's very well written and covers a vast amount of algorithms, which are also implemented in a very efficient and still understandable fashion.

Answer 4

There were problems with the C-Code submitted by Luis. First the encuentraValorMasProximo was missing. Second, the array reservation is one number short. I have also cleaned up the function. Below is the functional C-Code with the encuentraValorMasProximo function (renamed findNearestNeighbourIndex).

#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);
}

Answer 5

You may look at the GSL (numerical scientific library). There are many Matlab-like functions, among them and one-dimensional interpolation.

I am on my phone now, sorry, can not provide link.

Answer 6

Are you aware of the Matlab Coder? It automatically generates c/c++ code from Matlab code. If you have that as part of your Matlab package, you could just run the interp1 function through it and see what Matlab spits out.

Answer 7

@user1097111, your code exists a bug, in function findNearestNeighbourIndex, it should be if ( newDist >= 0 && newDist < dist ), not if ( newDist > 0 && newDist < dist ).

Answer 8

see lininterp1f in fileexchange.

Answer 9

If this helps someone in the future, this is a version without temporary arrays and without the 0 bug.

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