C mise en œuvre de Matlab interp1 fonction (interpolation linéaire)

J'ai porté Luis du code c++. Il semble fonctionner, mais je n'ai pas vérifié, donc, être conscient et re-vérifier vos résultats.

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

OriginalL'auteur Nikita Rokotyan

Excelent implémentations de fonctions fréquemment utilisées peuvent être trouvés dans le livre Numerical Recipes in C, qui est visible gratuitement en ligne. Chapitre 3.1.2 a une interpolation linéaire de la recette, le reste de ce chapitre couvre les plus avancés.

Je peux que recommander fortement ce livre, c'est très bien écrit et couvre une grande quantité d'algorithmes, qui sont également mis en œuvre dans une très efficace et toujours compréhensible.

OriginalL'auteur jupp0r

Vous pouvez regarder la GSL (bibliothèque scientifique numérique). Il y a beaucoup de Matlab-like fonctions, parmi eux, et à une dimension de l'interpolation.

Je suis sur mon téléphone maintenant, désolé, ne peut pas fournir de lien.

OriginalL'auteur Ivan Angelov

@user1097111, votre code existe un bug, en fonction findNearestNeighbourIndex, il devrait l'être, si ( newDist >= 0 && newDist < dist ), pas si ( newDist > 0 && newDist < dist ).

OriginalL'auteur Haiwei Liu

Si cela aide quelqu'un dans le futur, c'est une version sans temporaires tableaux et sans le 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;
}

OriginalL'auteur canesin