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 ("xi", "yi"), 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{"xi"}, max{"xi"}]. Une manière de faire consiste à utiliser une interpolation ; Scilab propose deux démarches :
Interpolation linéaire.
Le nuage de "n" points ("xi", "yi")1 ≤ "i" ≤ "n" est décrit par deux vecteurs : codice_1, qui contient les "xi" classés par ordre croissant, et codice_2, qui contient les "yi" correspondant.
Pour un interpolation linéaire, on considère une valeur codice_3 arbitraire, ou un vecteur de valeurs, les valeurs de codice_4 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 codice_5 : 
yp = interp1(x, y, xp, "nearest")
On peut aussi utiliser la commande d'interpolation à "n" dimensions codice_6 :
yp = linear_interpn(xp, x, y)
Si l'on veut faire une "extrapolation", c'est-à-dire avoir une valeur hors de l'intervalle ["x"1, "xn"], on peut ajouter une option :
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 codice_11 (les dérivées), soit trois vecteurs ( codice_1, codice_2, codice_11). Le vecteur codice_11 est obtenu par la commande codice_16 :
d = splin(x, y)
si l'on a maintenant une valeur codice_3 arbitraire, ou un vecteur de valeurs, les valeurs de codice_4 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 codice_1 et codice_2 sont « exactes ». Si le signal est bruité, il faut alors effectuer un lissage.
Scilab propose la commande codice_21 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 codice_22 (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 "zi" à des positions déterminées ("xi", "yi"), 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 :
L'interpolation linéaire se fait de la même manière que précédemment : on définit un nouveau couple de valeurs codice_3 et codice_4 (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 :
Et comme dans le plan, la fonction d'interpolation est décrite par (codice_1, codice_2, codice_33), les valeurs coordonnées des points du quadrillage et une matrice de coefficients codice_33 obtenue par
C = splin2d(x, y, z)
Si l'on considère un point du plan (codice_3, codice_4) (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 ("xi", "yi", "zi") sont décrits par une matrice M de trois colonnes et de "n" lignes
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 (codice_3, codice_4) (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"