Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

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

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

like image 445
Luis Andrés García Avatar asked Feb 22 '12 12:02

Luis Andrés García


People also ask

How does interp1 work in Matlab?

The interp1 command interpolates between data points. It finds values at intermediate points, of a one-dimensional function that underlies the data. This function is shown below, along with the relationship between vectors x , Y , xi , and yi . Interpolation is the same operation as table lookup.

How do you linear interpolate in Matlab?

The function to perform linear interpolation, MATLAB provides the interp1() function. Here, a sample point is a set of data points, which could be an array or a vector. The value of the unknown function on sample points is also a set that has the same size/length as the sample points.

Can you use interp1 to extrapolate outside the data range?

If the time at which we want an interpolated value is not within the range of times in the input time series, then interp1 will return NaN. That is, this function will not extrapolate outside of the range of input ( t ) values given.

Why does interp1 return NaN?

For the 'nearest' , 'linear' , and 'v5cubic' methods, interp1(x,Y,xi,method) returns NaN for any element of xi that is outside the interval spanned by x . For all other methods, interp1 performs extrapolation for out of range values.


2 Answers

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;
}
like image 67
Nikita Rokotyan Avatar answered Sep 20 '22 14:09

Nikita Rokotyan


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

like image 28
Luis Andrés García Avatar answered Sep 18 '22 14:09

Luis Andrés García