Python pour le calcul scientifique/Interpolation, extrapolation et lissage
Rappelons que dorénavant les programmes commencent tous par :
#!/usr/bin/python3
import numpy as np
import matplotlib.pyplot as plt
Interpolation à une dimension
modifier
- avec une valeur constante ;
- extrapolation périodique ;
- extrapolation linéaire.
Nous disposons de n points M de coordonnées (xpi, ypi)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 (xpi, ypi). Entre ces points, c'est-à-dire pour une valeur de x comprise entre deux valeurs xpi et xpi + 1, on utilise typiquement une des solutions suivantes :
- on détermine le xpj le plus proche de x et l'on prend y = ypj ce qui donne une fonction en escalier ;
- on trace un segment de droite entre (xpi, ypi) et (xpi + 1, ypi + 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 [xp 0 ; xp n – 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 (xp 0 à gauche, xp n – 1 à droite), soit une valeur fixée arbitrairement (éventuellement 0) ;
- soit en prolongeant la fonction avec la loi sur le segment [xp 0 ; xp 1] à gauche, [xp n – 2 ; xp n – 1] à droite ;
- soit en considérant que la fonction est périodique de période xp n – 1 – xp 0.
La fonction numpy.interp()
modifier
La fonction np.interp() permet de faire des interpolations linéaires :
y = np.interp(x, xp, yp)
où
xpetypsont des séquences (vecteurs, listes, n-uplets) de réels, la liste des points décrivant la fonction ;xest la séquence de valeurs auxquelles on s'intéresse.
Par défaut, l'extrapolation se fait par des valeurs constantes, yp 0 et yp n – 1. On peut indiquer que la fonction a une période avec l'argument period. Par exemple :
import numpy as np
import matplotlib.pyplot as plt
# **********************
# * numpy.interp() *
# * avec option period *
# **********************
xp = np.linspace(0, 2*np.pi, 5) # Cinq points
yp = np.sin(xp)
x = np.linspace(-2*np.pi, 4*np.pi, 25)
y = np.interp(x, xp, yp, period=2*np.pi) # interpolation linéaire
plt.plot(x, y, "y-")
plt.plot(xp, yp, "r+")
Le module scipy.interpolate
modifier
Le module SciPy fournit d'autres méthodes. Pour cela, il faut charger le module :
from scipy import interpolate
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 :
f = interpolate.interp1d(xp, yp)
y = f(x)
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 xp 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 yp correspondant au xp 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 (xp, yp) :
f = interpolate.BarycentricInterpolator(xp, yp)
y = f.__call__(x)
On a le même résultat avec la fonction interpolate.barycentric_interpolate() :
y = interpolate.barycentric_interpolate(xp, yp, x)
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 xpi dans le vecteur xp, et de mettre ypi, y’pi dans le vecteur yp. Si une valeur xpi figure n fois, on définit les dérivées jusqu'à l'ordre n – 1.
f = interpolate.KroghInterpolator(xp, yp)
y = f.__call__(x)
y = interpolate.Krogh_interpolate(xp, yp, x)
On peut avoir la dérivée d'ordre n avec la syntaxe :
f = interpolate.KroghInterpolator(xp, yp)
y_n = f.derivative(x, n) # dérivée n-ième
y_N = f.derivatives(x, N) # dérivées première à n-ième ;
# notez le « s » au nom de la fonction
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() :
f = interpolate.PchipInterpolator(xp, yp) # renvoie NaN en cas d'extrapolation
f = interpolate.PchipInterpolator(xp, yp, extrapolate=True) # autorise l'extrapolation
y = f.__call__(x) # on peut ajouter extrapolate=True lors de l'évaluation
y_n = f.__call__(x, n) # valeurs de la dérivée n-ième
f_n = f.derivative(n) # construit la fonction dérivée n-ième
F_n = f.antiderivative(n) # construit la n-ième primitive
racines = f.roots() # racines du polynôme
y = interpolate.pchip_interpolate(xi, yi, x)
y_n = interpolate.pchip_interpolate(xi, yi, x, n) # dérivée n-ième
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 :
f = interpolate.Akima1DInterpolator(xp, yp)
y = f.__call__(x)
f_n = f.derivative(n) # construit la fonction dérivée n-ième
F_n = f.antiderivative(n) # construit la n-ième primitive
racines = f.roots() # racines du polynôme
Voici à titre d'exemple le code source de l'image ci-contre.
Code source
import numpy as np
from scipy import interpolate
import matplotlib.pyplot as plt
xp = np.linspace(0, 2*np.pi, 5)
yp = np.sin(xp)
# ***************
# * CubicSpline *
# ***************
finterpa = interpolate.CubicSpline(xp, yp)
xa = np.linspace(-0.5*np.pi, 2.5*np.pi, 57)
ya = finterpa(xa)
# **********************
# * Interp1d quadratic *
# **********************
finterpb = interpolate.interp1d(xp, yp, kind="quadratic", fill_value="extrapolate")
xb = np.linspace(-0.5*np.pi, 2.5*np.pi, 49)
yb = finterpb(xb)
# *********
# * Pchip *
# *********
finterpf = interpolate.PchipInterpolator(xp, yp)
xf = np.linspace(-0.5*np.pi, 2.5*np.pi, 49)
yf = finterpf.__call__(xf)
# ***********
# * Akima1D *
# ***********
finterpg = interpolate.Akima1DInterpolator(xp, yp)
finterpg.extrapolate = True
xg = np.linspace(-0.9*np.pi, 2.9*np.pi, 97)
yg = finterpg.__call__(xg)
# *********
# * Tracé *
# *********
fig = plt.figure(figsize=[20, 12])
plt.plot(xb, yb, label="interp1d quadratic")
plt.plot(xa, ya, label="CubicSpline")
plt.plot(xf, yf, label="PchipInterpolator")
plt.plot(xg, yg, label="Akima1DInterpolator")
plt.plot(xp, yp, "k+", markersize=16, markeredgewidth=1.5)
plt.legend()
plt.savefig("comparaison_interpolation_scipy_interpolate_polynomes.svg", format="svg")
print("CubicSpline == interp1d cubic == BarycentricInterpolator == KroghInterpolator == splrep PPoly")
Lissage
modifierAvec le module scipy.signal
modifier
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) :
import numpy as np
import scipy.signal as signal
import matplotlib.pyplot as plt
epsilon = 0.05 # amplitude du bruit
x = np.arange(0, 1, 0.02) # points d'abscisse
# Génération d'un signal gaussien bruité
y = 0.1*np.exp(-(10*(x - 1)*x + 0.125)) + epsilon*np.random.normal(x)
# Lissage
ylisse = signal.medfilt(y, 7)
# Affichage
fig = plt.figure(figsize = [10, 6])
plt.plot(x, y, "b.")
plt.plot(x, ylisse, "k-", linewidth="0.5")
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 :
# Lissage
ylisse = signal.savgol_filter(y, 7, 3)
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 ws ; on doit avoir 0 < ws < 1 (1 est la fréquence de Nyquist).
Par exemple :
# Lissage
(b, a) = signal.butter(2, 0.3)
zi = signal.lfilter_zi(b, a)
ylisse = signal.lfilter(b, a, y, zi=zi*y[0])
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 ws est alors exprimée dans la même unité et l'on doit avoir 0 < ws < ƒs. La fréquence ws correspond à une atténuation d'un facteur .
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 ws = 5 et donc nous générons le filtre de Butterworth d'ordre 1 par :
# Lissage
(b, a) = signal.butter(1, 5, fs=len(y))
...
Avec le module scipy.fft
modifier
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 :
import scipy.fft as fft
# Lissage
ff = fft.fftshift(fft.fft(y))) # calcule la transformée de Fourier rapide et ordonne les coefficients
ff[-25:] = 0 # atténuation du spectre aux extrémités
ff[:25] = 0
ylissef = fft.ifft(fft.ifftshift(ff)) # réordonne les coefficients et calcule transformée de Fourier inverse
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 :
# Lissage
L = len(y) # calcule la longueur du signal
ff = fft.fftshift(fft.fft(y, fft.next_fast_len(L, real=True))) # calcule la transformée de Fourier rapide de manière optimisée
# en ajoutant des 0 en queue de signal, et ordonne les coefficients
ff[-25:] = 0 # atténuation du spectre aux extrémités
ff[:25] = 0
ylissef = fft.ifft(fft.ifftshift(ff))[:L] # réordonne les coefficients et calcule la transformée de Fourier inverse,
# et élimine les 0 ajoutés artificiellement pour optimiser
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 :
# Lissage
L = len(y) # calcule la longueur du signal
ff = fft.rfft(y, fft.next_fast_len(L, real=True)) # calcule la transformée de Fourier rapide de manière optimisée
# en ajoutant des 0 en queue de signal
ff[-25:] = 0 # atténuation du spectre à la fin
ylissef = fft.irfft(ff)[:L] # calcule la transformée de Fourier inverse,
# et élimine les 0 ajoutés artificiellement pour optimiser
numpy dispose également d'un moule fft qui reprend les fonctions de scipy.fft. Cependant, Il ne dispose pas de toutes les fonctions, il manque en particulier les fonctions d'optimisation (fft.next_fast_len()) et de réordonnancement (fft.fftshift() et fft.ifftshift()).