TP1 traitement numérique du signal

Indications : les TPs ne serviront dans l’évaluation de ce cours qu’à des points bonus. Il faudra à la fin des 3 heures de TP, envoyer une feuille Jupyter (.ipynb) à chupin@ceremade.dauphine.fr. Cette feuille doit être un compte rendu de votre travail, il vous faudra donc à partir de la feuille de sujet, rédiger vos réponses pour expliquer ce que vous avez produit. Ne laissez dans cette feuille de compte rendu que ce qui est utile !

Remarque : L'utilisation de Jupyter est proposée pour simplifier l'utilisation de Python par l'utilisation de feuilles de calculs : *les notebooks*. Bien évidemment, celles et ceux qui voudront utiliser un autre outil, coder simplement avec Python dans un fichier texte avec son éditeur favoris, etc., peuvent tout à fait le faire.

En première cellule de nos notebooks, on chargera les librairies dont nous avons besoin, en particulier numpy que l'on va présenter très rapidement dans la suite.

import numpy as np
import matplotlib
import matplotlib.pyplot as plt

Transformée de Fourier et théorème de Shannon

Nous allons commencer sur un exemple simple. Pour cela, générons un signal sinusoïdal, représentons-le, puis étudions son spectre dans le domaine de Fourier.

Remarque : Ici, puisqu'il s'agit de travaux pratiques informatiques, nous ne pouvons travailler qu'avec des signaux déjà échantillonnés et quantifiés.

Exercice : On veut générer la deux listes $t_n$ et $s_n=\sin(2 \pi f_0 t_n)$ pour $n\in[0,\dots,N-1]$ ($t_n$ est la discrétisation de l'intervalle $[0,1]$).

1. Quelle est la fréquence d'échantillonnage ?
2. En utilisant la fonction np.arange générez la suite (de type numpy.ndarray) $(t_n)$ et gérérez la suite $(s_n)$ en utilisant la fonction np.sin (qui peut prendre en argument une liste). On définira une variable fo égale à 3 et une variable N égale à 1000.

Pour comprendre un peu mieux ce que l'on fait, nous allons représenter notre signal. Pour cela on utilisera la librairie matplotlib.

On peut alors faire appel aux outils de tracés. Pour tracer une liste y en fonction d'une autre de même taille x, avec légende, il suffit d'exécuter :

plt.figure() # pour générer la figure
plt.plot(x,y)
plt.title('$y$ en fonction de $x$')
plt.xlabel("Abscisse")
plt.ylabel("Ordonnées")

Exercice : Tracer le signal défini précedemment.

Exercice : En utilisant la fonction de transformée de Fourier rapide (FFT) de numpy (np.fft.fft), calculer puis représentez le module de la transformée de Fourier du signal discrétisé $(s_n)$. On utilisera la fonction np.fft.fftfreq(<nbr points>,<time step>) pour générer la liste des fréquences pour pouvoir tracer la FFT. Voir https://docs.scipy.org/doc/numpy/reference/generated/numpy.fft.fftfreq.html#numpy.fft.fftfreq et https://docs.scipy.org/doc/numpy/reference/generated/numpy.fft.fft.html#numpy.fft.fft. Retrouve-t-on la fréquence du signal ?

Exercice : On dispose maintenant d'un signal sauvé dans le fichier signal.npz à télécharger ici, que l'on chargera sous forme d'un np.ndarray avec les lignes suivantes :

signal=np.load('signal.npz') # 

Ce fichier contient les valeurs du signal temporel nommé en variable x que l'on obtiendra grâce à

s=signal['x'].flatten() # 

Ce signal est échantillonné à la fréquence $F_e=32$.

1. Quel est le pas de temps $T_e$ entre deux échantillons ?
2. Construire la suite $(t_n)$ de la même taille que le signal correspondant au temps. Tracer le signal.
3. Faite l'analyse spectrale grâce à la FFT, on utilisera un nombre 8 fois supérieur à la taille du signal temporel pour le calcul de la FFT.

Exercice : Créez un nouveau signal en répertant le signal lu 8 fois (fonction np.tile). Faites-en ensuite l'analyse spectrale. On n'hésitera pas à zomer sur le spectre pour bien voir se qu'il se passe (plt.xlim), on le comparera au spectre du signal non périodisé en le superposant lors du tracé : il suffit de faire plusieurs plt.plot sans ouvrir de nouvelle plt.figure. On pourra utiliser la fonction plt.stem au lieu de plt.plot pour le spectre du signal périodisé.

Remarque : il y a plusieurs définition de la TFD, et l'amplitude peut être normalisé ou non. Il est possible qu'il faille pour bien comparer les deux spectres, renormaliser un des deux spectre.

Nous allons maintenant sous échantillonner notre signal et voir ce qu'il se produit sur son spectre.

Exercice :

1. Soit $(s_k)$, $k=0,\dots,N-1$, un signal discrétisé avec $N$ pair. À partir de la définition de la transformée de Fourier discrète : $$\hat s_k=\sum_{n=0}^{N-1}s_n\mathrm{e}^{-2\mathrm{i}\pi \frac nN},$$ montrer que la transformée de Fourier du signal $(s'_n)$, défini par $s'_k=s_{2k}$ pour tout $k=0,\dots,(N-2)/2$, on a $$\forall k=0,\dots, (N-2)/2,\quad \hat s'_k=\sum_{n=0}^{(N-2)/2}s_{2n}\mathrm{e}^{-2\mathrm{i}\pi \frac{2n}{N}}.$$
2. Construire des signaux échantillonnés à 16, 8 et 4Hz à partir du signal de départ (on pourra utiliser une boucle). On fera pour cela une fonction python qui échantillonne : d'après la question 1, tout en gardant la même taille de signal (128 ici), on peut extraire seulement certaines valeurs en complétant les autres par zéro (en prenant une valeur sur 2, 4 et 8 pour les différents cas).
3. Observer ce qu'il se produit sur le spectre. Que constate-t-on ?

Transformée de Fourier pour les images

On rappelle qu’on utilise alors la Transformée de Fourier Discrète (TFD) qui s'obtient par :

$$\hat{f}[k] = \frac{1}{N} \sum_{j=0}^{N-1} f[j] e^{-2i\pi j \frac{k}{N}}$$

pour $f$ un signal discret à $N$ valeurs.

Pour un signal bidimensionnel (comme une image) défini sur une grille $H\times W$, on peut généraliser cette définition et écrire la Transformée de Fourier Discrète 2D :

$$\hat{f}[k_x, k_y] = \frac{1}{H W} \sum_{p=0}^{H-1} \sum_{q=0}^{W-1} f[p, q] \exp\left(-2i\pi \left(\frac{p \cdot k_x}{H} + \frac{q \cdot k_y}{W}\right)\right)$$

Autrement dit, la transformée de Fourier d'une image $I$ est une image complexe $\hat{I}$ de mêmes dimensions.

On chargera les bibliothèques suivantes:

import matplotlib.pyplot as plt
import numpy as np
import skimage

plt.rcParams['figure.figsize'] = (10, 10)

Nous allons utiliser deux images d'exemple pour mieux interpréter l'effet de la transformée de Fourier :

• un rectangle blanc sur fond noir,
• une photo en niveaux de gris.

Pour cela, on utilisera le code suivant.

from skimage import data, color

img = color.rgb2gray(data.astronaut())
rect = np.zeros((100, 100))
rect[45:55, 40:60] = 1.
fig = plt.figure()
fig.add_subplot(121) and plt.imshow(rect, cmap="gray")
fig.add_subplot(122) and plt.imshow(img, cmap="gray") and plt.show()

Nous allons utilisez les fonctions fft2 et ifft2 correspondent ainsi à la transformée de Fourier rapide en 2D et à son inverse.

# Transformée de Fourier sur les images
from scipy.fft import fft2, ifft2
from scipy.fft import fftshift, ifftshift

fourier = fft2(rect)
print(f"La transformée de Fourier de l'image du carré est une matrice "
      f"de type {fourier.dtype} et de dimensions {fourier.shape}")

La transformée de Fourier de l'image du carré est une matrice de type complex128 et de dimensions (100, 100)

La transformée de Fourier produit un résultat en complexe. Ce qui nous intéresse le plus souvent est de savoir quelle est l'amplitude de chacun de ses coefficients (c'est-à-dire la norme du complexe associé) car celle-ci est liée à l'énergie de la fréquence considérée. La norme du complexe s'obtient avec np.abs.

fourier_amplitude = np.abs(fourier)
plt.imshow(fourier_amplitude, cmap="magma") and plt.show()

Cette visualisation n'est pas idéale car l'origine $(0,0)$ se trouve en haut à gauche (remarquez que les coins sont symétriques ! comme pour une transformée de Fourier 1D où l'on observe une symétrie par rapport à l'origine). La fonction fftshift permet de décaler la visualisation de sorte à mettre l'origine au centre de l'image.

On peut alors représenter la transformée avec la fonction fftshift https://numpy.org/doc/stable/reference/generated/numpy.fft.fftshift.html.

fourier_amplitude = np.abs(fftshift(fourier))
plt.imshow(fourier_amplitude, cmap="magma") and plt.show()

Il existe aussi la fonction ifftshift pour calculer la transformée de Fourier inverse à partir d’une transformée de Fourier shiftée par fftshift.

Question: De quoi pouvez-vous rapporcher cette image?

Maintenant que nous comprenons un peu mieux la généralisation de la transformée de Fourier aux signaux 2D, passons à notre image d’astronaute. Comme nous l'avons dit, la transformée de Fourier est une image complexe, on ne peut pas la visualiser directement. On peut s'inspirer des techniques de visualisation de spectre dans le cas 1D. On va donc calculer la norme de chaque coefficient du spectre 2D et en afficher le logarithme pour avoir une sorte d’équivalent des décibels.

Question : Représenter la transformée de Fourier de l’image d’astronaute.

Un (léger) problème de cette façon de faire est que la transformée de Fourier fait l'hypothèse que le signal observé correspond à une période d'un signal périodique. Ce n'est bien entendu pas le cas pour les images mais cela introduit une discontinuité (et donc des hautes fréquences) aux bords de l'image. La « croix » est provoquée par ces discontinuités. Une façon de s'en débarasser est d'utiliser un fenêtrage qui réduit l'impact des bords.

from skimage.filters import window
from skimage.util import img_as_float
image = img_as_float(img)

wimage = image * window('hann', image.shape)
plt.imshow(wimage, cmap="gray") and plt.show()

Question: Représenter la transformée de Fourier de cette nouvelle image? Est-ce que notre stratégie a fonctionné ? Pourquoi ?

Convolution avec un noyau

On rappelle que pour deux signaux réels, on note $*$ le produit de convolution défini par : $$(f*g)(x) = \int_{-\infty}^{+\infty} f(x-t)g(t) dt$$

et sa version discrète : $$(f*g)[k] = \sum_{p=-\infty}^{+\infty} f[x - p] g[p]$$ Une propriété intéressante de la transformée de Fourier est qu'elle transforme la convolution en produit, c'est-à-dire qu'en notant $\mathcal{F}$ la transformée de Fourier :

$$\mathcal{F}(f * g) = \mathcal{F}(f) \cdot \mathcal{F}(g)$$ Dans le cas 2D (pour des images, donc), on peut écrire la convolution bidimensionnelle comme le produit entre une image $I$ de dimensions $(H, W)$ et un noyau $K$ de dimensions $(k_h, k_w)$. On définit alors le filtre $\mathcal{K}$ comme étant le résultat de la convolution de $I$ par $K$:

$$\mathcal{K}(I)[m,n] = K * I[m, n] = \sum_{i=-p}^{+p} \sum_{j=-q}^{+q} I[m - i, n -j] \cdot K[i,j]$$

où $p = \frac{k_h - 1}{2}$ et $q = \frac{k_w - 1}{2}$.

Filtrage d'images numériques

Le filtrage est un élément fondamental du traitement du signal. Comme nous l'avons vu précédemment, une image numérique est un signal 2D discrétisé sur lequel nous pouvons utiliser toutes les techniques classiques d'analyse.

Gradients

Si l'on considère une image en niveaux de gris $I$ comme une fonction $(x,y) \rightarrow p \in \mathbb{R}$, alors on peut considérer les gradients d'une image, c'est-à-dire ses dérivées partielles selon $x$ et selon $y$. On note conventionnellement $\nabla I$ le gradient de $I$:

$$\nabla I=\begin{bmatrix} g_{x} \\ g_{y} \end{bmatrix} = \begin{bmatrix} \frac{\partial f}{\partial x} \\ \frac{\partial f}{\partial y} \end{bmatrix} $$

Néanmoins, comme $I$ prend ses valeurs sur une grille discrète (et finie) de pixels, il n'est pas possible de calculer le gradient réel. En revanche, il est possible d'utiliser la méthode des différences finies pour estimer le gradient :

$$g_x[i,j] = I[i+1, j] - I[i, j]$$ $$g_y[i,j] = I[i, j + 1] - I[i, j]$$

Le plus souvent, on va s'intéresser à la norme du gradient, c'est-à-dire à : $$\mathbf{G} = \| \nabla I \| = \sqrt{g_x^2 + g_y^2}$$

Pour des questions de symétrie et rester centré en $(i,j)$, on préfèrera les formules suivantes: $$g_x[i,j] = I[i+1, j] - I[i-1, j]$$ $$g_y[i,j] = I[i, j + 1] - I[i, j-1]$$

Pour la dérivée partielle horizontale, ceci peut être formulé à l’aide d’un filtre $3\times1$ avec un 0 au centre : $$G:=\begin{bmatrix}-1\\0\\+1\end{bmatrix}$$

Question : Implémenter une fonction qui réalise le filtre par produit de convolution d’une image en niveau de gris.

def conv2d(I,G):
    n,m = I.shape
    (kp,kq) = G.shape
    p=int((kp-1)/2)
    q=int((kq-1)/2)
    print(p)
    print(q)    
    O = np.zeros_like(I)
    for i in range(n):
        for j in range(m):
            for u in range(-p, p+1):
                for v in range(-q, q+1):

                    O[i,j] += (
                        I[(i+u) % n, (j+v) % m]
                        * G[u+p, v+q]
                    )
    return O

G = np.array([[-1],[0],[1]])
gradImage = conv2d(image,G)
plt.imshow(gradImage, cmap="gray") and plt.show()

1
0

our calculer le gradient, on peut réaliser la convolution avec les noyaux ci-dessus mais on va généralement utiliser un noyau de la forme :

$$\mathbf{G}_{\text{horizontal}} = \begin{bmatrix} -1 & -1 & -1 \\ 0 & 0 & 0 \\ +1 & +1 & +1 \end{bmatrix} * \mathbf{I}$$

Ce noyau est appelé filtre de Prewitt (du nom de son inventrice).

Pour l'appliquer à notre image, scipy fournit des outils de base pour réaliser une convolution en 2D:

from scipy.signal import fftconvolve, convolve2d

Question : en utilisant convolve2d (https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.convolve2d.html), calculer la convolution de l’image de l’astronaute avec le noyau de Prewitt. Afficher l’image obtenue.

top_prewitt = np.array([[1, 1, 1], [0, 0, 0], [-1, -1, -1]])
filtered = convolve2d(img, top_prewitt)
plt.imshow(filtered, cmap="gray") and plt.show()

Question : représenter la transformée de Fourier de l’image obtenue.

fourierA = fft2(filtered)
fourierA_amplitude = np.abs(fftshift(fourierA))
plt.imshow(fourierA_amplitude, cmap="magma") and plt.show()

Question : représenter les valeurs de l’image ainsi filtrée qui dépasse un certain seuil (on pourra utiliser 0.3).

plt.imshow(filtered > 0.3, cmap="gray") and plt.show()

Filtre fréquentiel associé au gradient

Si on considère la fonction image en niveau de gris

$$ f : \mathbb R^2 \to \mathbb R $$

et sa transformée de Fourier :

$$ \hat f(\xi_x,\xi_y) = \int_{\mathbb R^2} f(x,y) e^{-i(x\xi_x+y\xi_y)} dx\,dy $$

Le gradient est :

$$ \nabla f = \left( \frac{\partial f}{\partial x}, \frac{\partial f}{\partial y} \right) $$

Maintenant, on peut regarder ce qu’il se passe dans le domaine fréquentielle.

On utilise alors la propriété fondamentale :

$$ \mathcal F\left( \frac{\partial f}{\partial x} \right) = i\xi_x \hat f $$

et :

$$ \mathcal F\left( \frac{\partial f}{\partial y} \right) = i\xi_y \hat f $$

Ce qui nous donne:

$$ \mathcal F(\nabla f) = \left( i\xi_x \hat f, i\xi_y \hat f \right) $$

Le gradient correspond donc au multiplicateur fréquentiel : $$ H(\xi_x,\xi_y) = i \begin{pmatrix} \xi_x \\ \xi_y \end{pmatrix} $$

Autrement dit :

$$ \nabla \quad \longleftrightarrow \quad i\xi $$

avec :

$$ \xi = \begin{pmatrix} \xi_x \\ \xi_y \end{pmatrix} $$

Mise en place pratique

Pour passer au domaine discret et implémenter le filtre fréquentiel correspondant à la dérivée partielle suivant $y$ (pour comparer avec l’exemple du filtre par convolution), il nous suffit de construire le filtre partiel: $$ \mathcal F\left( \frac{\partial f}{\partial y} \right) = i\xi_y \hat f $$

Question : Implémenter la fonction python dont la signature est:

{python}
def gy_fourier(I):

qui pourra suivre la structure suivante:

1. récupération des dimensions de l’image I
2. calcule de la transformée de Fourier de l’image (np.fft.fft2)
3. calcule des fréquences en $x$ et en $y$ avec la fonction np.fft.fftfreq() qui donne à partir de la dimension le tableau des fréquences discrètes correspondantes (d’abord les positives, puis les négatives)
4. fabrication de la matrice par répétition des fréquences (horizontalement et verticalement) avec np.meshgrid
5. fabrication du filtre Hy par multiplication de la transformée de Fourier de l’image par $i\xi_y$
6. obtention de la matrice de l’image correspondante par prise de la valeur réelle de la transformée de Fourier inverse (np.fft.ifft2)

def gy_fourier(I):
    n, m = I.shape

    # FFT de l'image
    F = np.fft.fft2(I)

    # Fréquences
    kx = 2 * np.pi * np.fft.fftfreq(m)
    ky = 2 * np.pi * np.fft.fftfreq(n)

    # Grille fréquentielle
    KX, KY = np.meshgrid(kx, ky)
    

    # Filtres de dérivation
    Hy = 1j * KY

    # Dérivées fréquentielles
    Fy = Hy * F

    # Retour espace
    gy = np.real(np.fft.ifft2(Fy))


    return gy

filtered2 = grad_squared_fourier(img)
plt.imshow(filtered2, cmap="gray") and plt.show()

Question: représenter la transformée de Fourier de l’image obtenue.

fourierA = fft2(filtered2)
fourierA_amplitude = np.abs(fftshift(fourierA))
plt.imshow(fourierA_amplitude, cmap="magma") and plt.show()

Le phénomène de Gibbs

On se propose de regarder numériquement l'approche par série de Fourier d'un signal rectangulaire périodique.

Exercice : Définissez un fonction python qui code un signal $s(\cdot)$ carré périodique de période $T=2$ et d'amplitude $A=3$, c'est-à-dire que pour $t\in[0,T/2[$, $s(t)=A$ et pour $t\in[T/2,T[$, $s(t)=-A$.

On pourra mettre ces valeurs en paramètres optionels de la fonction. Celle-ci pourra avoir la forme suivante :

def carre(t,A=3,T=2):
    # instruction 
    return valeur

On pourra utiliser le modulo pour gérer la périodicité de la fonction : pour deux réels $a$ et $b$, le modulo en python s'écrit a%b.

Exercice : Représentez cette fonction. On propose de compléter le code suivant :

N=100
t=np.arange(N)*4/N
liste = []
for i in t:
    # à compléter

en utilisant la fonction .append et la transformation de liste en np,ndarray.

Exercice : Montrez que la série de Fourier de $s$ est $$s(t)=\sum_{k=0}^{+\infty}\frac{-4A}{(2k+1)\pi}\sin\left((2k+1)t\frac{2\pi}{T}\right),$$ puis pour différentes valeurs de plus en plus élevées de $N$ tracer les fonctions $$s_N(t)=\sum_{k=0}^{N}\frac{-4A}{(2k+1)\pi}\sin\left((2k+1)t\frac{2\pi}{T}\right).$$ Qu'observe-t-on ?

Exercice :

1. Implémentez une fonction qui calcule numériquement la décomposition en série de Fourier (jusqu'à au rang $N$ fixé) d'une fonction réelle périodique quelconque. On utilisera les fonction d'intégration numérique de numpy.
2. Vérifiez sur le cas de la fonction créneau que votre fonction est correcte.