Analyse fréquentielle d'un signal par transformée de Fourier¶

La transformée de Fourier permet de représenter le spectre de fréquence d'un signal non périodique.

Note

Cette partie s'intéresse à un signal à une dimension.

Signal à une dimension¶

Un signal unidimensionnel est par exemple le signal sonore. Il peut être vu comme une fonction définie dans le domaine temporel :

$x:t\rightarrow x(t)$

Dans le cas du traitement numérique du signal, ce dernier n'est pas continu dans le temps, mais échantillonné.
Le signal échantillonné est obtenu en effectuant le produit du signal $x(t)$ par un peigne de Dirac de période $Te$ :

$x_e(t)=x(t)\sum\limits_{k=-\infty}^{+\infty}\delta(t-kT_e)$

Attention

La fréquence d'échantillonnage d'un signal doit respecter le théorème de Shannon-Nyquist qui indique que la fréquence $Fe$ d'échantillonnage doit être au moins le double de la fréquence maximale $f$ du signal à échantillonner :

$Fe>2\times f$

Transformée de Fourier Rapide (notée FFT)¶

La transformée de Fourier rapide est un algorithme qui permet de calculer les transformées de Fourier discrète d'un signal échantillonné.
On note pour la suite $X(f)$ la FFT du signal $x_e(t)$ .

Il existe plusieurs implantations dans Python de la FFT :

Ici nous allons utiliser numpy.fft pour calculer les transformées de Fourier.

FFT d'un sinus¶

Création du signal et échantillonnage¶

import numpy as np
import matplotlib.pyplot as plt


def x(t):
    # Calcul du signal x(t) = sin(2*pi*t)
    return np.sin(2*np.pi*t)


# Échantillonnage du signal
Durée = 1  # Durée du signal en secondes
Te = 0.1  # Période d'échantillonnage en seconde

N = int(Durée/Te) + 1  # Nombre de points du signal échantillonné

te = np.linspace(0, Durée, N)  # Temps des échantillons
t = np.linspace(0, Durée, 2000)  # Temps pour le signal non échantillonné

x_e = x(te)  # Calcul de l'échantillonnage

# Tracé du signal
plt.scatter(te, x_e, color='orange', label="Signal échantillonné")
plt.plot(t, x(t), '--', label="Signal réel")
plt.grid()
plt.xlabel(r"$t$ (s)")
plt.ylabel(r"$x(t)$")
plt.title(r"Échantillonnage d'un signal $x(t$)")
plt.legend()
plt.show()

png

Cas extrême où f=Fe¶

import numpy as np
import matplotlib.pyplot as plt


def x(t):
    # Calcul du signal x(t) = sin(2*pi*t)
    return np.sin(2*np.pi*t)


# Échantillonnage du signal
Durée = 1  # Durée du signal en secondes
Te = 1/2  # Période d'échantillonnage en seconde

N = int(Durée/Te) + 1  # Nombre de points du signal échantillonné

t_echantillons = np.linspace(0, Durée, N)  # Temps des échantillons
t = np.linspace(0, Durée, 2000)  # Temps pour le signal non échantillonné

# Tracé du signal
plt.scatter(t_echantillons, x(t_echantillons),
            color='orange', label="Signal échantillonné")
plt.plot(t, x(t), '--', label="Signal réel")
plt.grid()
plt.xlabel(r"$t$ (s)")
plt.ylabel(r"$x(t)$")
plt.title(r"Échantillonnage d'un signal $x(t$) à $Fe=2\times f$")
plt.legend()
plt.show()

png

Calcul de la transformée de Fourier¶

# Création du signal

import numpy as np
import matplotlib.pyplot as plt


def x(t):
    # Calcul du signal x(t) = sin(2*pi*t)
    f = 1  # Fréquence du signal
    A = 1  # Amplitude du signal
    return A*np.sin(2*np.pi*f*t)


# Échantillonnage du signal
Durée = 3  # Durée du signal en secondes
Te = 0.01  # Période d'échantillonnage en seconde
N = int(Durée/Te) + 1  # Nombre de points du signal échantillonné
te = np.linspace(0, Durée, N)  # Temps des échantillons
x_e = x(te)

# Tracé du signal
plt.scatter(te, x_e, label="Signal échantillonné")
plt.grid()
plt.xlabel(r"$t$ (s)")
plt.ylabel(r"$x(t)$")
plt.title(r"Signal échantillonné")
plt.show()

png

from numpy.fft import fft, fftfreq

# Calcul FFT
X = fft(x_e)  # Transformée de fourier
freq = fftfreq(x_e.size, d=Te)  # Fréquences de la transformée de Fourier

plt.subplot(2,1,1)
plt.plot(freq, X.real, label="Partie réel")
plt.plot(freq, X.imag, label="Partie imaginaire")
plt.grid()
plt.legend()
plt.xlabel(r"Fréquence (Hz)")
plt.ylabel(r"Amplitude $X(f)$")
plt.title("Transformée de Fourier")

plt.subplot(2,1,2)
plt.plot(freq, X.real, label="Partie réel")
plt.plot(freq, X.imag, label="Partie imaginaire")
plt.xlim(-2,2) # Limite autour de la fréquence du signal
plt.grid()
plt.legend()
plt.xlabel(r"Fréquence (Hz)")
plt.ylabel(r"Amplitude $X(f)$")
plt.title("Transformée de Fourier autour de la fréquence du signal")
plt.tight_layout()
plt.show()

png

Mise en forme des résultats¶

La mise en forme des résultats consiste à ne garder que les fréquences positives et à calculer la valeur absolue de l'amplitude pour obtenir l'amplitude du spectre pour des fréquences positives.

L'amplitude est ensuite normalisée par rapport à la définition de la fonction fft.

# On prend la valeur absolue de l'amplitude uniquement pour les fréquences positives
X_abs = np.abs(X[:N//2])
# Normalisation de l'amplitude
X_norm = X_abs*2.0/N
# On garde uniquement les fréquences positives
freq_pos = freq[:N//2]

plt.plot(freq_pos, X_norm, label="Amplitude absolue")
plt.xlim(0, 10)  # On réduit la plage des fréquences à la zone utile
plt.grid()
plt.xlabel(r"Fréquence (Hz)")
plt.ylabel(r"Amplitude $|X(f)|$")
plt.title("Transformée de Fourier")
plt.show()

png

Cas d'un fichier audio¶

On va prendre le fichier audio suivant Cri Wilhelm au format wav et on va réaliser la FFT de ce signal.

import scipy.io.wavfile as wavfile
import matplotlib.pyplot as plt
import numpy as np

# Lecture du fichier
rate, data = wavfile.read('0477.wav')

x = data[:, 0]  # Sélection du canal 1
# Création de instants d'échantillons
t = np.linspace(0, data.shape[0]/rate, data.shape[0])

# Tracé du signal
plt.plot(t, x, label="Signal échantillonné")
plt.grid()
plt.xlabel(r"$t$ (s)")
plt.ylabel(r"Amplitude")
plt.title(r"Signal sonore")
plt.show()

png

from numpy.fft import fft, fftfreq

# Calcul FFT
X = fft(x)  # Transformée de fourier
freq = fftfreq(x.size, d=1/rate)  # Fréquences de la transformée de Fourier

# Calcul du nombre d'échantillon
N = x.size

# On prend la valeur absolue de l'amplitude uniquement pour les fréquences positives et normalisation
X_abs = np.abs(X[:N//2])*2.0/N
# On garde uniquement les fréquences positives
freq_pos = freq[:N//2]

plt.plot(freq_pos, X_abs, label="Amplitude absolue")
plt.xlim(0, 6000)  # On réduit la plage des fréquences à la zone utile
plt.grid()
plt.xlabel(r"Fréquence (Hz)")
plt.ylabel(r"Amplitude $|X(f)|$")
plt.title("Transformée de Fourier du Cri Whilhelm")
plt.show()

png

Spectrogramme d'un fichier audio¶

On repart du même fichier audio que précédemment.
Le spectrogramme permet de visualiser l'évolution des fréquences du signal au cours du temps.

import scipy.signal as signal
import scipy.io.wavfile as wavfile
import matplotlib.pyplot as plt
import numpy as np

# Lecture du fichier
rate, data = wavfile.read('0477.wav')

x = data[:, 0]  # Sélection du canal 1
# Création de instants d'échantillons
#t = np.linspace(0, data.shape[0]/rate, data.shape[0])

# Calcul du spectrogramme
f, t, Sxx = signal.spectrogram(x, rate)
# On limite aux fréquences présentent
Sxx_red = Sxx[np.where(f<6000)]
f_red = f[np.where(f<6000)]

# Affichage du spectrogramme
plt.pcolormesh(t, f_red, Sxx_red, shading='gouraud')
plt.ylabel('Fréquence (Hz)')
plt.xlabel('Temps (s)')
plt.title('Spectrogramme du Cri Whilhem')
plt.show()

png

Spectrogramme d'une mesure¶

On réalise une mesure d'accélération à l'aide d'un téléphone, qui peut mesurer par exemple les vibrations dues à un séisme. Et on va visualiser le spectrogramme de cette mesure.

Le fichier de mesure est le suivant Raw_Data.csv.

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

# Lecture du fichier
# Lecture des en-têtes des données avec comme délimiteur le point-virgule
head = np.loadtxt('Raw_Data.csv', delimiter=',', max_rows=1, dtype=np.str)
# Lecture des données au format float
data = np.loadtxt('Raw_Data.csv', delimiter=',', skiprows=1)

# print(head)

# Sélection de la colonne à traiter
x = data[:, 3]
te = data[:, 0]
Te = np.mean(np.diff(te))

# Calcul du spectrogramme
f, t, Sxx = signal.spectrogram(x, 1/Te, window=signal.get_window('hann',32))

# On limite aux fréquences présentent
freq_lim = 11
Sxx_red = Sxx[np.where(f < freq_lim)]
f_red = f[np.where(f < freq_lim)]

# Affichage

# Signal d'origine
plt.plot(te, x)
plt.ylabel('accélération (m/s²)')
plt.xlabel('Temps (s)')
plt.title('Signal')
plt.grid()
plt.show()

# Affichage du spectrogramme
plt.plot(te, [0]*len(x))
plt.pcolormesh(t, f_red, Sxx_red, shading='gouraud')
plt.ylabel('Fréquence (Hz)')
plt.xlabel('Temps (s)')
plt.title('Spectrogramme')
plt.show()

png

Attention

Ici vous remarquerez le paramètre window=signal.get_window('hann',32) qui a été rajouté lors du calcul du spectrogramme. Il permet de définir la fenêtre d'observation du signal, le chiffre 32 désigne ici la largeur (en nombre d'échantillons) d'observation pour le calcul de chaque segment du spectrogramme.