C implementation of the Matlab interp1 function (linear interpolation)

Question

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

Answer 1

I myself implemented this linear interpolation (some of them are written in Spanish, sorry). A function called encuentraValorMasProximo simply finds the closest value (elementoMasProximo) and index (indiceEnVector) to another (xx [i]) in the 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 may 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 result :

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 have ported Luis code to C ++. It seems to work, but I have not tested it, so be aware and 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

Experimental implementations of commonly used functions can be found in the book Numerical Recipes in C , which is available for viewing on the Internet. Chapter 3.1.2 presents a linear interpolation recipe; the rest of the chapter is devoted to more advanced ones.

I can highly recommend this book, it is very well written and covers a huge number of algorithms that are also implemented in a very efficient and understandable way.

Answer 4

There were problems with the C code provided by Louis. At first, encuentraValorMasProximo was missing. Secondly, array reservation is short. I also cleared 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 can look at the GSL (Digital Science Library). There are many Matlab-like functions, including one-dimensional interpolation.

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

Answer 6

Do you know what Matlab Coder is ? It automatically generates c / C ++ code from Matlab code. If you have this as part of your Matlab package, you can simply run the interp1 function through it and see what Matlab drinks.

Answer 7

@ user1097111, there is an error in your code, the findNearestNeighbourIndex function should have 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 error 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; }