Linear interpolation in C++

[RezaAb]

Hi Guys,

I am trying to linearly interpolate a 2D vector. I tried different methods offered by the people, but I could not find an appropriate answer.
Here is an example:

I have the 2D vector:
A[2][5]={{1,5,10,15,20}}{{0.3,0.5,0.8,0.10,0.14}}
and also I have the following 1D vector:
B[20]={{1,2,3,4,5,6,...,17,18,19,20}}

So since I am not very familiar with C++ programming, my question is how can I write an inline function for solving this issue?

Thanks,

[jonnin]

not sure what your data is meaning..

if the user types in (???) you want to get (???) back out, from which (or both) of the vectors?

LI is not too hard -- you treat an array as a function/lookup for data points, you get the point above and below the desired point, and compute the answer. If its outside your bounds, you can either reject it or continue the last line from the last 2 points.

[RezaAb]

Thanks for the reply, the final answer which I am trying to evaluate is a vector similar to B including interpolated values based on the second row of the A. It might be a 2D vector:

B'[2][20]={{1,2,3,4,5,6,...,17,18,19,20}}{{0.3,0.35,0.4,0.5,....0.8,0.10,0.13,0.14}}

I was wondering if we have a library for this. I am trying to be efficient in time!

Thanks.

[lastchance]

#include <iostream>
#include <iomanip>
#include <vector>
using namespace std;

//======================================================================

// Returns interpolated value at x from parallel arrays ( xData, yData )
//   Assumes that xData has at least two elements, is sorted and is strictly monotonic increasing
//   boolean argument extrapolate determines behaviour beyond ends of array (if needed)

double interpolate( vector<double> &xData, vector<double> &yData, double x, bool extrapolate )
{
   int size = xData.size();

   int i = 0;                                                                  // find left end of interval for interpolation
   if ( x >= xData[size - 2] )                                                 // special case: beyond right end
   {
      i = size - 2;
   }
   else
   {
      while ( x > xData[i+1] ) i++;
   }
   double xL = xData[i], yL = yData[i], xR = xData[i+1], yR = yData[i+1];      // points on either side (unless beyond ends)
   if ( !extrapolate )                                                         // if beyond ends of array and not extrapolating
   {
      if ( x < xL ) yR = yL;
      if ( x > xR ) yL = yR;
   }

   double dydx = ( yR - yL ) / ( xR - xL );                                    // gradient

   return yL + dydx * ( x - xL );                                              // linear interpolation
}

//======================================================================

int main()
{
   // Original data
   vector<double> xData = { 1, 5, 10, 15, 20 };
   vector<double> yData = { 0.3, 0.5, 0.8, 0.1, 0.14 };

   // Set up some points for interpolation in xVals
   const int NPTS = 20;
   vector<double> xVals, yVals;
   for ( int i = 1; i <= NPTS; i++ ) xVals.push_back( (double)i );

   // Interpolate
   for ( double x : xVals )
   {
      double y = interpolate( xData, yData, x, true );
      yVals.push_back( y );
   }

   // Output
   #define SP << fixed << setw( 15 ) << setprecision( 6 ) <<
   #define NL << '\n'
   cout << "Original data:\n";
   for ( int i = 0; i < xData.size(); i++ ) cout SP xData[i] SP yData[i] NL;
   cout << "\nInterpolated data:\n";
   for ( int i = 0; i < xVals.size(); i++ ) cout SP xVals[i] SP yVals[i] NL;
}

Original data:
       1.000000       0.300000
       5.000000       0.500000
      10.000000       0.800000
      15.000000       0.100000
      20.000000       0.140000

Interpolated data:
       1.000000       0.300000
       2.000000       0.350000
       3.000000       0.400000
       4.000000       0.450000
       5.000000       0.500000
       6.000000       0.560000
       7.000000       0.620000
       8.000000       0.680000
       9.000000       0.740000
      10.000000       0.800000
      11.000000       0.660000
      12.000000       0.520000
      13.000000       0.380000
      14.000000       0.240000
      15.000000       0.100000
      16.000000       0.108000
      17.000000       0.116000
      18.000000       0.124000
      19.000000       0.132000
      20.000000       0.140000

[jonnin]

if I understood this

you have
y = f(x)
where X is currently 1,5,10,15 (which is annoying, because 1-5 is not the same as 5-10 .. its +4, +5, +5, +5, ...?!)
and you want the Y data for
X = 1,2,3.... N

Is that correct?

I don't know of anything built into the language to do it, but what I don't know can fill a lot of books. Im sure there is a library out there that does this, its a common task. I couldn't say which one to pick. If you want to DIY, I can probably help there.

[RezaAb]

Thanks, Jonnin. Lastchance's code solved my problem.