Implementazione C della funzione interp1 di Matlab (interpolazione lineare)

Question

Conoscete qualche implementazione C della funzione interp1 di Matlab (solo una 'lineare')? Ne conosco uno per Java.

Answer 1

Ho eseguito personalmente questa interpolazione lineare (in parte è scritta in spagnolo, mi spiace). La funzione chiamata encuentraValorMasProximo trova solo il valore più vicino (elementoMasProximo) e l'indice (indiceEnVector) in un altro (xx [i]), in un 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]; 
} 
} 

La funzione prova potrebbe essere:

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

} 

E sua esito:

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

Si può guardare la GSL (biblioteca scientifica numerica). Ci sono molte funzioni simili a Matlab, tra queste e l'interpolazione unidimensionale.

Sono sul mio telefono ora, mi dispiace, non posso fornire il collegamento.

Answer 3

Sei al corrente dello Matlab Coder? Genera automaticamente codice c/C++ dal codice Matlab. Se lo hai come parte del tuo pacchetto Matlab, puoi semplicemente eseguire la funzione interp1 attraverso di esso e vedere ciò che Matlab sputa.

Answer 4

Implementazioni eccellenti delle funzioni di uso comune sono disponibili nel libro Numerical Recipes in C, che è possibile visualizzare gratuitamente online. Il capitolo 3.1.2 ha una ricetta di interpolazione lineare, il resto del capitolo riguarda quelli più avanzati.

Posso fortemente raccomandare questo libro, è scritto molto bene e copre una grande quantità di algoritmi, che sono anche implementati in modo molto efficiente e ancora comprensibile.

Answer 5

Ho portato il codice di Luis in C++. Sembra che funzioni, ma non l'ho controllato molto, quindi fai attenzione e ricontrolla i risultati.

#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 6

vedere lininterp1f in fileexchange.

Answer 7

Ci sono stati problemi con il codice C inviato da Luis. Prima mancava l'encuentraValorMasProximo. In secondo luogo, la prenotazione dell'array è un numero breve. Ho anche ripulito la funzione.Di seguito è riportato il codice C funzionale con la funzione encuentraValorMasProximo (rinominato 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 8

@ user1097111, il codice esiste un bug, in funzione di findNearestNeighbourIndex, dovrebbe essere se (newDist> = 0 & & newDist < dist), se non (newDist> 0 & & newDist < dist).

Answer 9

Se questo aiuta qualcuno in futuro, questa è una versione senza array temporanei e senza il bug 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; 
}