C Implementation Of Matlab Interp1 Function (linear Interpolation)

Answer :

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

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

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.