# Python pour le calcul scientifique/Interpolation, extrapolation et lissage

Rappelons que dorénavant les programmes commencent tous par :

## Interpolation à une dimension

Nous disposons de n points M de coordonnées (x$_{pi}$, y$_{pi}$)$_{0≤i≤n–1}$. Nous supposons que ces points décrivent une fonction continue ƒ et nous voulons connaître une approximation de ƒ(x) pour un point x quelconque. Les méthodes d'interpolation utilisent des fonctions qui passent exactement par les points (x$_{pi}$, y$_{pi}$). Entre ces points, c'est-à-dire pour une valeur de x comprise entre deux valeurs x$_{pi}$ et x$_{pi+1}$, on utilise typiquement une des solutions suivantes :
* on détermine le x$_{pj}$ le plus proche de x et l'on prend y = y$_{pj}$ ce qui donne une fonction en escalier ;
* on trace un segment de droite entre (x$_{pi}$, y$_{pi}$) et (x$_{pi+1}$, y$_{pi+1}$), le point (x, y) est situé sur ce segment de droite ; c'est la méthode d'interpolation linéaire, la fonction est affine par parties, son graphe est un polygone ouvert ;
* on utilise le polynôme de degré n – 1 passant pour tous les points, le point (x, y) est situé sur ce polynôme ;
* on utilise des splines (cerces), des polynômes par partie, en général de degré 2 ou 3, dont les dérivées sont raccordées.
En dehors de l'intervalle [x$_{p0}$ ; x$_{pn–1}$], on peut extrapoler les valeurs de la manière suivante :
* pas d'extrapolation, la fonction retourne des valeurs NaN ;
* extrapolation par une valeur constante, soit la valeur la plus proche (x$_{p0}$ à gauche, x$_{pn–1}$ à droite), soit une valeur fixée arbitrairement (éventuellement 0) ;
* soit en prolongeant la fonction avec la loi sur le segment [x$_{p0}$ ; x$_{p1}$] à gauche, [x$_{pn–2}$ ; x$_{pn–1}$] à droite ;
* soit en considérant que la fonction est périodique de période x$_{pn–1}$ – x$_{p0}$.

### La fonction numpy.interp

La fonction np.interp permet de faire des interpolations linéaires :
où
* xp et yp sont des séquences (vecteurs, listes, n-uplets) de réels, la liste des points décrivant la fonction ;
* x est la séquence de valeurs auxquelles on s'intéresse.
Par défaut, l'extrapolation se fait par des valeurs constantes, y$_{p0}$ et y$_{pn–1}$. On peut indiquer que la fonction a une période avec l'argument period. Par exemple :

### Le module scipy.interpolate

Le module SciPy fournit d'autres méthodes. Pour cela, il faut charger le module :
La fonction interpolate.interp1d (interpolation à une dimension, l'avant-dernière lettre est le chiffre 1) crée une fonction d'interpolation. Pour avoir une interpolation linéaire comme précédemment, on écrit :
On peut autoriser l'extrapolation linéaire avec l'option
    fill_value="extrapolate"
Pour extrapoler avec une valeur constante, par exemple des zéros, il faut utiliser deux paramètres :
    fill_value=0, bounds_error=False
En absence de ces paramètres, une valeur de x en dehors de l'intervalle des x$_p$ génère une erreur. Si l'on veut extrapoler avec les valeurs aux extrémités, on écrit :
    fill_value=(yp[0], yp[-1]), bounds_error=False
Cette fonction propose d'autre méthodes d'interpolation : avec l'option kind= :
* "linear" : méthode par défaut ;
* "nearest" : renvoie la valeur de y$_p$ correspondant au x$_p$ le plus proche ;
* "zero" : fait une interpolation par une spline d'ordre 0 ;
* "slinear" : " par une spline d'ordre 1 ;
* "quadratic" : " par une spline d'ordre 2 ;
* "cubic" : " par une spline d'ordre 3.
Par exemple :
    interpolate.interp1d(xp, yp, fill_value="extrapolate", kind="nearest")
La classe interpolate.BarycentricInterpolator correspond à une fonction polynôme passant par un ensemble de points (x$_p$, y$_p$) :
On a le même résultat avec la fonction interpolate.barycentric_interpolate :
La classe interpolate.KroghInterpolator et la fonction interpolate.krogh_interpolate permettent de fixer les dérivées successives du polynôme : pour définir la dérivée au point i, il suffit de mettre deux fois la valeur x$_{pi}$ dans le vecteur x$_p$, et de mettre y$_{pi}$, y’$_{pi}$ dans le vecteur y$_p$. Si une valeur x$_{pi}$ figure n fois, on définit les dérivées jusqu'à l'ordre n – 1.
On peut avoir la dérivée d'ordre n avec la syntaxe :
On peut avoir une interpolation par un polynôme cubique par parties PCHIP (piecewise cubic hermite interpolating polynomial) avec la classe interpolate.PchipInterpolator ou la fonction interpolate.pchip_interpolate :
La classe interpolate.Akima1DInterpolator définit une fonction d'interpolation selon la méthode Akima : elle construit des sous-splines à partir de polynômes cubiques par partie :
Voici à titre d'exemple le code source de l'image ci-contre.
Code source

## Lissage
### Avec le module scipy.signal

La bibliothèque SciPy dispose de fonctions de traitement du signal, dont des fonctions de lissage.
La fonction la plus simple est scipy.signal.medfilt(y, h) qui, pour chaque point du tenseur y, calcule la médiane sur une fenêtre de largeur h points ; la valeur par défaut pour la fenêtre est h = 3. Par exemple, pour une courbe (filtrage d'un tenseur d'ordre 1) :
Ce module dispose également d'un filtre de Savitzky-Golay : scipy.signal.savgol_filter(y, h, p) effectue un lissage avec une fenêtre de largeur h et un polynôme de degré p. Dans l'exemple ci-dessus, pour filtrer avec un polynôme de degré 3 et une fenêtre de largeur 7 points, on écrit :
L'algorithme permet également d'obtenir une dérivée n-ième (n inférieur à p) avec la syntaxe scipy.signal.savgol_filter(y, h, p,n )
Le lissage peut être vu comme un filtre passe-bas : le bruit correspond à des fréquences élevées. On peut donc utiliser d'autres fonctions de filtre.
Un des plus simples est le filtre de Butterworth. Il faut dans un premier temps générer les caractéristiques du filtre, avec la fonction scipy.signal.butter ; puis déterminer les conditions initiales du filtre avec scipy.signal.lfilter_zi ; enfin appliquer ce filtre au signal avec la fonction scipy.signal.lfilter. Un filtre de Butterworth a deux paramètres :
* l'ordre, qui détermine la vitesse d'atténuation lorsque la fréquence s'élève ; un filtre d'ordre 1 décroît à –6 dB/octave, un filtre d'ordre 2 à –12 dB/octave, etc.
* la fréquence de coupure w$_s$ ; on doit avoir 0 < w$_s$ < 1 (1 est la fréquence de Nyquist).
Par exemple :
Si l'on veut donner un sens physique à la fréquence de coupure, il faut préciser la fréquence d'échantillonnage ƒ$_s$ (sample frequency) ; la fréquence de coupure w$_s$ est alors exprimée dans la même unité et l'on doit avoir 0 < w$_s$ < ƒ$_s$. La fréquence w$_s$ correspond à une atténuation d'un facteur ${\sqrt {2}}$.
Dans l'exemple précédent, nous avons len(y) = 50 points ; en supposant que la fenêtre représente une seconde, nous avons une fréquence d'échantillonnage de ƒ$_s$ = 50/1 = 50. Si nous voulons filtrer des phénomènes qui se répètent cinq fois dans la fenêtre, nous prenons w$_s$ = 5 et donc nous générons le filtre de Butterworth d'ordre 1 par :

### Avec le module scipy.fft

Puisque le lissage est essentiellement un filtre passe-bas, nous pouvons envisager d'utiliser la transformée de Fourier et de mettre à 0 les fréquences élevées du spectre avant de faire la transformée de Fourier inverse. Le plus simple et le plus rapide consiste à utiliser la transformée de Fourier rapide (FFT, fast Fourier transform). La fonction scipy.fft.fft mise en œuvre dans le module fft de SciPy place cependant les fréquences élevées au centre du spectre ; pour faciliter le filtrage, il faut décaler les coefficients (shift) pour mettre les fréquences élevées à l'extérieur (au début et à la fin du vecteur), on peut ainsi facilement les mettre à 0 par tranchage (slicing). Cela utilise les fonctions scipy.fft.fftshift pour réordonner les coefficients, puis scipy.fft.ifftshift pour les remettre dans l'ordre initial.
Le code de base est donc :
Par ailleurs, la transformée de Fourier rapide est plus efficace pour certaines longueurs de signal. La fonction scipy.fft.next_fast_len permet de trouver la longueur de signal optimale la plus proche. Le code optimisé devient donc :
Dans le cas présent, le signal est à valeurs réelles et non pas complexes. Nous pouvons donc utiliser les fonctions scipy.fft.rfft et scipy.fft.irfft. Dans ce cas-là, la transformée de Fourier n'a pas de coefficients à indice négatif, il suffit donc de mettre à zéro les derniers coefficients, sans avoir besoin de changer leur ordre. Le code est alors :

## Notes et références

Statistiques < ↑ > Régression et optimisation