Rappelons que dorénavant les programmes commencent tous par :
import numpy as np
import matplotlib.pyplot as plt
Interpolation à une dimension.
Nous disposons de "n" points M de coordonnées ("x"p"i", "y"p"i")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"p"i", "y"p"i"). Entre ces points, c'est-à-dire pour une valeur de "x" comprise entre deux valeurs "x"p"i" et "x"p"i" + 1, on utilise typiquement une des solutions suivantes :
En dehors de l'intervalle ["x"p 0 ; "x"p "n" – 1], on peut extrapoler les valeurs de la manière suivante :
La fonction codice_1.
La fonction codice_2 permet de faire des interpolations linéaires :
y = np.interp(x, xp, yp)
où
Par défaut, l'extrapolation se fait par des valeurs constantes, "y"p 0 et "y"p "n" – 1. On peut indiquer que la fonction a une période avec l'argument codice_6. Par exemple :
import numpy as np
import matplotlib.pyplot as plt
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 codice_7.
Le module SciPy fournit d'autres méthodes. Pour cela, il faut charger le module :
from scipy import interpolate
La fonction codice_8 (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
Pour extrapoler avec une valeur constante, par exemple des zéros, il faut utiliser deux paramètres :
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 :
Cette fonction propose d'autre méthodes d'interpolation : avec l'option codice_12 :
Par exemple :
La classe codice_20 correspond à une fonction polynôme passant par un ensemble de points ("x"p, "y"p) :
f = interpolate.BarycentricInterpolator(xp, yp)
y = f.__call__(x)
On a le même résultat avec la fonction codice_21 :
y = interpolate.barycentric_interpolate(xp, yp, x)
La classe codice_22 et la fonction codice_23 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"p"i" dans le vecteur "x"p, et de mettre "y"p"i", "y"’p"i" dans le vecteur "y"p. Si une valeur "x"p"i" 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 "()" avec la classe codice_24 ou la fonction codice_25 :
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 codice_26 définit une fonction d'interpolation selon la méthode Akima : elle construit des sous-"" à 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.
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)
finterpa = interpolate.CubicSpline(xp, yp)
xa = np.linspace(-0.5*np.pi, 2.5*np.pi, 57)
ya = finterpa(xa)
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)
finterpf = interpolate.PchipInterpolator(xp, yp)
xf = np.linspace(-0.5*np.pi, 2.5*np.pi, 49)
yf = finterpf.__call__(xf)
finterpg = interpolate.Akima1DInterpolator(xp, yp)
finterpg.extrapolate = True
xg = np.linspace(-0.9*np.pi, 2.9*np.pi, 97)
yg = finterpg.__call__(xg)
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.
Avec le module codice_27.
La bibliothèque SciPy dispose de fonctions de traitement du signal, dont des fonctions de lissage.
La fonction la plus simple est codice_28 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
y = 0.1*np.exp(-(10*(x - 1)*x + 0.125)) + epsilon*np.random.normal(x)
ylisse = signal.medfilt(y, 7)
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 : codice_29 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 :
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 codice_30
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 codice_31 ; puis déterminer les conditions initiales du filtre avec codice_32 ; enfin appliquer ce filtre au signal avec la fonction codice_33. Un filtre de Butterworth a deux paramètres :
Par exemple :
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 "()" ; 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 formula_1.
Dans l'exemple précédent, nous avons codice_34 = 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 codice_35.
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, ""). La fonction codice_36 mise en œuvre dans le module codice_37 de SciPy place cependant les fréquences élevées au centre du spectre ; pour faciliter le filtrage, il faut décaler les coefficients "()" 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 "()". Cela utilise les fonctions codice_38 pour réordonner les coefficients, puis codice_39 pour les remettre dans l'ordre initial.
Le code de base est donc :
import scipy.fft as fft
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 codice_40 permet de trouver la longueur de signal optimale la plus proche. Le code optimisé devient donc :
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
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,
Dans le cas présent, le signal est à valeurs réelles et non pas complexes. Nous pouvons donc utiliser les fonctions codice_41 et codice_42. 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 :
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
ff[-25:] = 0 # atténuation du spectre à la fin
ylissef = fft.irfft(ff)[:L] # calcule la transformée de Fourier inverse,
Notes et références.
Statistiques < ↑ > Régression et optimisation