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

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_p 0 ; x_{p 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 (x_p 0 à gauche, x_{p n – 1} à droite), soit une valeur fixée arbitrairement (éventuellement 0) ;
soit en prolongeant la fonction avec la loi sur le segment [x_p 0 ; x_p 1] à gauche, [x_{p n – 2} ; x_{p n – 1}] à droite ;
soit en considérant que la fonction est périodique de période x_{p n – 1} – x_p 0.

La fonction `numpy.interp()`

La fonction np.interp() permet de faire des interpolations linéaires :

y = np.interp(x, xp, yp)

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_p 0 et y_{p 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`

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 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) :

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 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.

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

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) :

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 w_s ; on doit avoir 0 < w_s < 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 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 :

# Lissage
(b, a) = signal.butter(1, 5, fs=len(y))
...

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 :

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

La bibliothèque 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()).

Notes et références

Statistiques < ↑ > Régression et optimisation

Interpolation à une dimension

La fonction numpy.interp()

Le module scipy.interpolate