Interpolation 1D avec Matlab (2024)

• Parler des jointures SQL en français
• Limites de suites et précision limitée - Partie 2

Si vous souhaitez, notamment,

• rééchantillonner des signaux,
• estimer une fonction coûteuse en un point à partir des valeurs en d’autres points,
• tracer des courbes lissées à partir d’un faible nombre de points,

alors vous aller avoir besoin d’interpoler, voire d’extrapoler !

La fonction interp1 de Matlab permet d’effectuer des interpolations et extrapolations à une dimension. Nous verrons dans ce tutoriel comment l’utiliser, et surtout comment choisir ses paramètres en fonction de vos besoins.

• Syntaxe
• Exemple travaillé
• Choisir les options adéquates

Syntaxe

Nous nous concentrerons sur l’usage avec tous les arguments, y compris optionnels.¹ Cet usage se résume par la ligne suivante :

1

vq = interp1(x, v, xq, method, extrapolation);

où :

• x est le vecteur des abscisses des points d’échantillonnage,
• v est le vecteur des valeurs d’échantillonnage, correspondant à x.
• xq est le vecteur des abscisses des points dont on souhaite interpoler la valeur,
• method est la méthode d’interpolation (à choisir notamment entre 'nearest', 'linear', 'pchip' et 'spline'),
• extrapolation est la stratégie d’extrapolation (soit 'extrap', soit une constante numérique),
• vq est le vecteur des valeurs interpolées.

Autrement dit, vq = interp1(x, v, xq, method, extrapolation) retourne la valeur vq aux points xq de la fonction valant v aux points x, et dont la valeur en dehors de ces points est calculée en utilisant la méthode d’interpolation method et la stratégie d’extrapolation extrapolation.

Exemple travaillé

Pour clarifier les explications de la section précédente, nous allons travailler sur un exemple.

Imaginons que l’on dispose de relevés de l’élévation du soleil dans le ciel toutes les heures entre 9 heures et 19 heures.

12

temps =[9 10 11 12 13 14 15 16 17 18 19];elevation =[28.5 40 50 60 67.5 70 65 56 45.5 35 23];

On souhaite estimer la position du soleil de 8 heures à 20 heures avec une résolution d’un quart d’heure, de sorte que le vecteur de temps que l’on souhaite est le suivant.

1

tempsInterp = 8:0.25:20;

Pour estimer l’élévation pour chaque valeur du vecteur précédent, il est nécessaire d’interpoler et d’extrapoler car les données ont seulement une résolution d’une heure et que les relevés n’ont été effectués que de 9 heures à 19 heures.

Pour guider notre choix de méthode d’interpolation et une stratégie d’extrapolation, on examine les informations dont on dispose sur la courbe. On sait que la fonction qui représente l’élévation en fonction du temps est:

• continue (le soleil ne saute pas de position),
• à dérivée continue (la vitesse d’élévation ne change pas brutalement),
• à dérivée seconde continue (si la vitesse d’élévation varie, alors elle varie continûment).

Les informations ci-dessous signifient que la fonction est au moins de classe $C^2$. La méthode d’interpolation 'spline' et la méthode d’extrapolation 'extrap' semblent des choix judicieux (voir la partie suivante pour plus de détails sur le choix d’une méthode d’interpolation).

1

elevationInterp = interp1(temps, elevation, tempsInterp, 'spline', 'extrap');

Visualisons le résultat obtenu !

123456789

figure;scatter(temps, elevation, 'o');hold all;scatter(tempsInterp, elevationInterp, 'x');xlabel('Temps (heures)');ylabel('Élévation (°)');title('Élévation du soleil au cours d''une journée');legend('mesures','interpolation');ylim([0 80]);

Choisir les options adéquates

Choisir les bonnes méthodes d’interpolation et d’extrapolation n’est pas trivial dans la plupart des cas. Cette section vise à vous guider pour choisir en connaissance de cause.

Méthode d’interpolation

nearest

Si vous souhaitez une très grande performance au détriment de la continuité, alors 'nearest' est fait pour vous ! Vous aurez la plus grande rapidité possible, mais la précision peut en pâtir, puisque la fonction d’interpolation n’est pas continue. Vous aurez des marches d’escalier. Elle nécessite au minimum deux points.

Vous retrouverez ci-dessous le script ayant généré le graphique ci-dessus.

 1 2 3 4 5 6 7 8 91011121314151617

f = @(x) abs(exp(-(x-2.5).^2).*(x-2).^3);nSamples = 17;xSamples = linspace(0, 10, nSamples);ySamples = f(xSamples);xInterp = linspace(-1, 11, 10000);yNearest = interp1(xSamples, ySamples, xInterp, 'nearest', NaN);figure;plot(xInterp, f(xInterp), '-m');hold all;stairs(xInterp, yNearest, 'b');scatter(xSamples, ySamples, '.k');xlim([0 10]);grid on;legend('original', 'interpolation','échantillons');xlabel('x');ylabel('y');

linear

La méthode 'linear' fournit une fonction d’interpolation continue. C’est un bon compromis entre performance et précision. Elle fournit une fonction linéaire par morceau, avec des cassures aux points de données. C’est la méthode à choisir si vous n’avez pas de contrainte particulière de continuité ou de performance. Elle nécessite au minimum deux points.

Vous retrouverez ci-dessous le script ayant généré le graphique ci-dessus.

 1 2 3 4 5 6 7 8 91011121314151617

f = @(x) abs(exp(-(x-2.5).^2).*(x-2).^3);nSamples = 17;xSamples = linspace(0, 10, nSamples);ySamples = f(xSamples);xInterp = linspace(-1, 11, 10000);yNearest = interp1(xSamples, ySamples, xInterp, 'linear', NaN);figure;plot(xInterp, f(xInterp), '-m');hold all;stairs(xInterp, yNearest, 'b');scatter(xSamples, ySamples, '.k');xlim([0 10]);grid on;legend('original', 'interpolation','échantillons');xlabel('x');ylabel('y');

pchip

Pour avoir une fonction d’interpolation à dérivée continue, il faut passer à la méthode 'pchip', qui se base sur des splines cubiques de Hermite. Cette méthode présente l’avantage de préserver la forme : sur chaque intervalle, la fonction d’interpolation est monotone, et donc ne dépasse pas les points de données. Par contre, les crêtes peuvent être coupées. Elle nécessite au minimum quatre points.

 1 2 3 4 5 6 7 8 91011121314151617

f = @(x) abs(exp(-(x-2.5).^2).*(x-2).^3);nSamples = 17;xSamples = linspace(0, 10, nSamples);ySamples = f(xSamples);xInterp = linspace(-1, 11, 10000);yNearest = interp1(xSamples, ySamples, xInterp, 'pchip', NaN);figure;plot(xInterp, f(xInterp), '-m');hold all;stairs(xInterp, yNearest, 'b');scatter(xSamples, ySamples, '.k');xlim([0 10]);grid on;legend('original', 'interpolation','échantillons');xlabel('x');ylabel('y');

spline

La méthode 'spline' est la plus lisse fournie par Matlab. Il s’agit d’un interpolateur à dérivée seconde continue. La fonction d’interpolation est ainsi très lisse, et permet de représenter assez bien les pics. Elle a malheureusem*nt tendance à créer des ondulations¹ si les données forment des plateaux. Elle nécessite au minimum quatre points.

Vous retrouverez ci-dessous le script ayant généré le graphique ci-dessus.

 1 2 3 4 5 6 7 8 91011121314151617

f = @(x) abs(exp(-(x-2.5).^2).*(x-2).^3);nSamples = 17;xSamples = linspace(0, 10, nSamples);ySamples = f(xSamples);xInterp = linspace(-1, 11, 10000);yNearest = interp1(xSamples, ySamples, xInterp, 'spline', NaN);figure;plot(xInterp, f(xInterp), '-m');hold all;stairs(xInterp, yNearest, 'b');scatter(xSamples, ySamples, '.k');xlim([0 10]);grid on;legend('original', 'interpolation','échantillons');xlabel('x');ylabel('y');

Synthèse

Le tableau ci-dessous résume les caractéristiques des différentes méthodes d’interpolation.

Autres méthodes

Matlab propose d’autres méthodes d’interpolation, je vous laisse regarder en détail la documentation.

Extrapolation

L’extrapolation est un sujet délicat. Dans la majorité des cas, je conseille d’éviter les extrapolation. Si vous êtes forcé à le faire, interp1 propose deux méthodes.

A partir de la méthode d’interpolation ('extrap')

La première consiste à extrapoler en se basant sur la méthode d’interpolation. Pour cela, il faut utiliser le paramètre 'extrap'.

Vous retrouverez ci-dessous le script ayant généré le graphique ci-dessus.

 1 2 3 4 5 6 7 8 91011121314151617

f = @(x) abs(exp(-(x-2.5).^2).*(x-2).^3);nSamples = 17;xSamples = linspace(0, 10, nSamples);ySamples = f(xSamples);xInterp = linspace(-1, 11, 10000);yNearest = interp1(xSamples, ySamples, xInterp, 'spline', 'extrap');figure;plot(xInterp, f(xInterp), '-m');hold all;stairs(xInterp, yNearest, 'b');scatter(xSamples, ySamples, '.k');xlim([-1 11]);grid on;legend('original', 'interpolation','échantillons');xlabel('x');ylabel('y');

Avec une valeur numérique constante

Autrement, il est possible de spécifier une valeur qui sera utilisée pour tous les points tombant en dehors de l’étendue des $x_i$.

Vous retrouverez ci-dessous le script ayant généré le graphique ci-dessus.

 1 2 3 4 5 6 7 8 91011121314151617

f = @(x) abs(exp(-(x-2.5).^2).*(x-2).^3);nSamples = 17;xSamples = linspace(0, 10, nSamples);ySamples = f(xSamples);xInterp = linspace(-1, 11, 10000);yNearest = interp1(xSamples, ySamples, xInterp, 'spline', 0);figure;plot(xInterp, f(xInterp), '-m');hold all;stairs(xInterp, yNearest, 'b');scatter(xSamples, ySamples, '.k');xlim([-1 11]);grid on;legend('original', 'interpolation','échantillons');xlabel('x');ylabel('y');

Cette utilisation correspond parfois à un cas où l’on ne contrôle pas parfaitement les données d’entrées (paramètres de l’utilisateur, sortie d’un autre script, etc). Il est dans ce cas intéressant de marquer ces valeurs grâce à la valeur spéciale NaN.

Vous savez désormais comment interpoler à une dimension avec Matlab, ainsi que choisir au mieux les paramètres en fonction de votre problème et de vos données.

Si êtes curieux et souhaitez en apprendre plus, les quelques liens suivant devraient élargir vos horizons.

Références

Sur le fonctionnement de Matlab

• Documentation de interp1 sur le site de The Mathworks.
• Article détaillant les différences entre pchips et splines, par Cleve Moler de The Mathworks.

Sur l’interpolation au plus proche

• Une implémentation simpliste en Python de l’interpolation au plus proche.

Sur l’interpolation linéaire

• Une implémentation simpliste en Python de l’interpolation linéaire.
• Plus de détails sur l’interpolation linéaire sur Wikipédia

Sur l’interpolation cubique de type pchip

• Une implémentation simpliste en Python de l’interpolation cubique de type pchip.
• Plus de détails sur les interpolations cubiques de type pchip sur Wikipédia.

Sur l’interpolation cubique de type spline

• Plus de détails sur les splines sur Wikipédia.

• Parler des jointures SQL en français
• Limites de suites et précision limitée - Partie 2