Découvrir Scilab/Interpolation, extrapolation et lissage

5. Interpolation, extrapolation et lissage

Interpolation dans le plan

Présentation du problème

On dispose de points expérimentaux (x_i, y_i), qui sont censés suivre une loi y = ƒ(x) inconnue, et l'on désire avoir une valeur y correspondant à une valeur x arbitraire, mais contenue dans l'intervalle [min{x_i}, max{x_i}]. Une manière de faire consiste à utiliser une interpolation ; Scilab propose deux démarches :

interpolation linéaire : si x_i ≤ x ≤ x_{i + 1}, alors on prend y sur le segment de droite [(x_i, y_i) ; ((x_{i + 1}, y_{i + 1})] :
$y=y_{i}+{\frac {y_{i+1}-y_{i}}{x_{i+1}-x_{i}}}(x-x_{i})$ ;
interpolation par un polynôme de degré trois, en utilisant des splines cubiques d'Hermite :
y = p(x)
la fonction p étant un polynôme de degré trois par parties, les portions de polynôme étant raccordées aux points expérimentaux (même valeur, même pente).

Pour plus de détails voir : wikipedia:fr:Interpolation linéaire et wikipedia:fr:Spline cubique d'Hermite.

Interpolation linéaire

Le nuage de n points (x_i, y_i)_{1 ≤ i ≤ n} est décrit par deux vecteurs : x, qui contient les x_i classés par ordre croissant, et y, qui contient les y_i correspondant.

Pour un interpolation linéaire, on considère une valeur xp arbitraire, ou un vecteur de valeurs, les valeurs de yp correspondantes sont obtenues par

yp = interp1(x, y, xp)

Plutôt que d'avoir les valeurs interpolées, on peut préférer les valeurs les plus proches ; pour cela, on utilise l'option "nearest" :

yp = interp1(x, y, xp, "nearest")

On peut aussi utiliser la commande d'interpolation à n dimensions linear_interpn() :

yp = linear_interpn(xp, x, y)

Si l'on veut faire une extrapolation, c'est-à-dire avoir une valeur hors de l'intervalle [x₁, x_n], on peut ajouter une option :

"by_zero" : y = 0 hors intervalle ;
"C0" : y = y₁ si x < x₁, y = y_n si x > x_n ;
"natural" : si la fonction est supposée monotone hors de l'intervalle de définition, on utilise le segment de droite le plus proche ;
"periodic" : si la fonction est supposée périodique, on répète le motif obtenu.

Par exemple :

yp = linear_interpn(xp, x, y, "natural")

Interpolation par une spline cubique

Pour une interpolation par spline cubique d'Hermite, la fonction d'interpolation est décrite par les coordonnées des points et un vecteur de coefficients d (les dérivées), soit trois vecteurs ( x, y, d). Le vecteur d est obtenu par la commande splin() :

d = splin(x, y)

si l'on a maintenant une valeur xp arbitraire, ou un vecteur de valeurs, les valeurs de yp correspondantes sont obtenues par

yp = interp(xp, x, y, d)

Lissage

L'interpolation considère que les données sont peu bruitées ; les points mesurés sont « fiables », les valeurs des x et y sont « exactes ». Si le signal est bruité, il faut alors effectuer un lissage.

Scilab propose la commande lsq_splin() qui utilise des splines cubique comme modèle, déterminées par la méthode des moindres carrés (least square). Si l'on a deux vecteurs x et y de m valeurs chacun, on choisit n abscisses xs pour déterminer la spline, avec dans l'idéal n << m. Typiquement :

// ***** Génération des données *****

m = 200; // nombre de points mesurés
bruit = 0.2; // amplitude du bruit
x = linspace(0, 2*%pi, m);
y = sin(x) + 0.2*rand(x, "normal"); // valeurs bruitées

// ***** Détermination de la spline *****

n = 6; // nombre de points de contrôle de la spline
xs = linspace(x(1), x($), n); // abscisse des points de contrôle
// répartition uniforme
[ys, d] = lsq_splin(x, y, xs); // détermination de la spline

// ***** Calcul des valeurs lissées de y *****

yliss = interp(x, xs, ys, d); // calcul des données lissées

// ***** Tracé *****

plot(x, y) // données originales
plot(x, yliss, "r") // données lissées

Ainsi, Scilab a construit une spline cubique passant par les abscisses xs (en l'occurrence 0, 2π/5, 4π/5, 6π/5, 8π/5 et 2π), et a déterminé les dérivées d en ces n points.

Si l'on veut éviter des effets de bord, on peut choisir une répartition différente des points de contrôle. On utilise fréquemment une répartition sinusoïdale pour avoir plus de points aux bords qu'au centre (comme pour les polynômes de Tchebychev).

theta = linspace(-%pi, 0, n); // angles uniformément répartis
xs = mean([x(1), x($)]) + 0.5*(x($) - x(1))*cos(theta);
// abscisse des points de contrôle

Ceci est en particulier important si l'on veut extrapoler, c'est-à-dire calculer des valeurs de y pour des x hors des points de mesure.

Cela ne présente d'intérêt que pour n ≥ 4.

L'utilisateur peut évidemment coder un autre algorithme. L'image ci-contre montre le calcul de la dérivée et de la dérivée seconde par l'algorithme de Savitzky-Golay ; cliquer sur l'image pour avoir accès à la page de description indiquant le code utilisé.

On peut aussi effectuer un simple lissage par la méthode des moyennes glissantes en utilisant la fonction convol() (voir ci-après). Si l'on appelle p le nombre de points de l'intervalle glissant :

p = 10; // largeur de l'intervalle

h = 1/p*ones(1, p); // fonction porte pour le lissage

yliss = conv(h, y);

plot(x, y) // données originales
plot(x, yliss(1+p/2:$-p/2+1), "r") // données lissées

Interpolation d'une surface

Présentation du problème

On dispose de valeurs z_i à des positions déterminées (x_i, y_i), qui sont censés suivre une loi z = ƒ(x, y) inconnue, ce qui définit donc une surface dans l'espace, et l'on désire avoir une valeur z correspondant à un point (x, y) arbitraire.

Interpolation linéaire sur un quadrillage

Les données sont décrites par :

un vecteur x de dimension m et un vecteur y de dimension n, strictement croissants, qui définissent un quadrillage du plan, donc un ensemble de m×n points {(x₁, y₁) ; (x₁, y₂) ; … ; (x₁, y_n) ; (x₂, y₁) ; … ; (x_m, y_n) ;
une matrice z de dimension m×n, contenant les valeurs correspondante : z_{i, j} = ƒ(x_i, y_j).

L'interpolation linéaire se fait de la même manière que précédemment : on définit un nouveau couple de valeurs xp et yp (ou un couple de vecteurs défnissant un nouveau quadrillage), et l'on a :

zp = linear_interp(xp, yp, x, y, z)

c'est-à-dire qu'à partir des points (x, y, z), on définit les points interpolés (xp, yp, zp).

Interpolation cubique sur un quadrillage

Comme précédemment, les données sont décrites par :

un vecteur x de dimension m et un vecteur y de dimension n, strictement croissants, qui définissent un quadrillage du plan, donc un ensemble de m×n points {(x₁, y₁) ; (x₁, y₂) ; … ; (x₁, y_n) ; (x₂, y₁) ; … ; (x_m, y_n) ;
une matrice z de dimension m×n, contenant les valeurs correspondante : z_{i, j} = ƒ(x_i, y_j).

Et comme dans le plan, la fonction d'interpolation est décrite par (x, y, C), les valeurs coordonnées des points du quadrillage et une matrice de coefficients C obtenue par

C = splin2d(x, y, z)

Si l'on considère un point du plan (xp, yp) (ou un couple de vecteurs défnissant un nouveau quadrillage), les valeurs de z sur ces nouveaux points sont :

zp = interp2d(xp, yp, x, y, C)

Interpolation cubique sur un nuage de points dispersés

On ne dispose pas toujours de relevés de valeurs sur des coordonnées régulières formant un quadrillage. Dans ce cas-là, on utilise la méthode de Shepard (pondération des plus proches voisins). Les n points expérimentaux (x_i, y_i, z_i) sont décrits par une matrice M de trois colonnes et de n lignes

\mathrm {M} ={\begin{pmatrix}x_{1}&y_{1}&z_{1}\\\vdots &\vdots &\vdots \\x_{n}&y_{n}&z_{n}\end{pmatrix}}

La fonction d'interpolation est définie par une liste typée T obtenue par

T = cshep2d(M)

Si l'on considère un point du plan (xp, yp) (ou un couple de vecteurs défnissant un nouveau quadrillage), les valeurs de z sur ces nouveaux points sont :

zp = eval_cshep2d(xp, yp, T)

Notes

Statistiques < ↑ > Régression