TD1 - Analyse en fréquence de signaux sonores.¶
1. Introduction à la FFT.¶
On va utiliser la fft sur un signal 1D pour en regarder le contenu en fréquentiel.
import numpy as np
import scipy as sp
import pylab as pl
# Choose a grid x and a signal y
n = 1000
m = int(n/2)
print(m)
x = np.linspace(-1,1,n)
y = np.zeros_like(x)
y[m:] = 1
y[0:m] = -1
pl.plot(x,y)
pl.title("Une fonction créneau")
500
Text(0.5, 1.0, 'Une fonction créneau')
import scipy.fftpack as sf
# On utilise la fonction fft pour calculer la transformée de Fourier.
spectre = sf.fft(y)
# on utilise la fonction ifft pour calculer l'inverse.
inverse_du_spectre = sf.ifft(spectre)
Exercice 1:¶
Calculer la norme \(L^2\) du signal et la norme \(L^2\) de sa transformée de Fourier.
On vérifiera que la fft et l’inverse de la ftt sont bien inverses numériques.
# illustration.
pl.figure(figsize = (12,6))
pl.subplot(1,2,1)
pl.plot(np.abs(spectre))
pl.title(" Le spectre du signal ")
pl.subplot(1,2,2)
pl.plot(np.abs(sf.fftshift(spectre)))
pl.title("La fréquence nulle est au milieu.")
Text(0.5, 1.0, 'La fréquence nulle est au milieu.')
normL2 = \( \sqrt{\sum_i{|y_i^2|}}\)
def normL2(y):
yn = y/y.size # normalisation
return np.sqrt(np.sum(np.abs(yn*yn)))
print('Norme du signal : ', normL2(y))
print('Norme du spectre : ', normL2(spectre))
print('Norme de (inverse_du_spectre-y) : ', normL2(inverse_du_spectre - y), '<< 1')
print("fft et l'inverse de la ftt sont bien inverses numériques")
Norme du signal : 0.0316227766016838
Norme du spectre : 1.0000000000000002
Norme de (inverse_du_spectre-y) : 9.738639465863477e-18 << 1
fft et l'inverse de la ftt sont bien inverses numériques
Exercice 2:¶
lisez le code ci-dessous et expliquez ce qui est obtenu.
reduced = np.copy(spectre)
h = 40
reduced[h:n-h] = 0
reconstruct=sf.ifft(reduced)
pl.plot(np.real(reconstruct))
[<matplotlib.lines.Line2D at 0x7fc3ba8c2d10>]
On obtient une réprensentation du signal filtré à l’aide d’un filtre passe-bas. On observe le phénomène de Gibbs.
Exercice 3:¶
Créer une gaussienne en utilisant la formule \(e^{-x^2/\sigma^2}\) et tracer la transformée de Fourier de cette fonction.
def gaussienne(x, sigma):
return np.exp(-np.power(x, 2)/sigma**2)
t = np.linspace(-1,1,1000)
signal = gaussienne(t, 0.05)
tf_gaussienne = sf.fft(signal)
pl.figure(figsize = (12,6))
pl.subplot(1,2,1)
pl.plot(signal)
pl.title("Signal de la gaussienne")
pl.subplot(1,2,2)
pl.plot(np.abs(sf.fftshift(tf_gaussienne)))
pl.title("Spectre de la gaussienne")
Text(0.5, 1.0, 'Spectre de la gaussienne')
On va créer un signal de type sinusoidal numérique, qu’on va ensuite écouter.
fs = 44100 # fq d'échantillonage
T = 2 # duree du signal en secondes
fq = 440
def generateT(fs, T):
return np.arange(0, T * fs)/fs
def sinus(t, amp, fq):
return amp*np.sin(2*np.pi*fq * t)
t = generateT(fs, T)
signal = sinus(t, 20, fq)
tf_signal = sf.fft(signal)
pl.plot(t, signal)
[<matplotlib.lines.Line2D at 0x7fc3ba715750>]
from IPython.display import Audio
Audio(10*y, rate = 44100)
Exercice 4:¶
Représenter le signal sur 2 dixièmes de seconde et compter le nombre de maximums pour en vérifier approximativement la fréquence du signal.
En utilisant la fft, représenter le contenu fréquentiel du signal. On s’attend à avoir un pic localisé en la fréquence utilisée pour la définition du signal.
signal_2ds = signal[:int(2*(signal.size/T/10))] # signal sur 2 dixiemes de seconde
tf_signal_2ds = sf.fft(signal_2ds)
pl.figure(figsize = (20,6))
pl.subplot(1,3,1)
pl.plot(signal)
pl.title("Signal")
pl.subplot(1,3,2)
pl.plot(signal_2ds)
pl.title("Zoom")
pl.subplot(1,3,3)
pl.plot(np.abs(sf.fftshift(tf_signal_2ds)))
pl.title("La TF du signal")
Text(0.5, 1.0, 'La TF du signal')
def countSpike(x, epsilon):
max = np.max(x)
return np.argwhere(x> max - epsilon).size
e= 10**(-2)
estimation_fq = 1/ (2/10 / countSpike(signal_2ds, e))
print("La fréquence pour laquelle la transformée de Fourier est maximale est {} Hz".format(estimation_fq))
print("C'est une bonne approximation de la fréquence réel égale à 440 Hz")
La fréquence pour laquelle la transformée de Fourier est maximale est 439.99999999999994 Hz
C'est une bonne approximation de la fréquence réel égale à 440 Hz
Exercice 5:¶
En utilisant la fonction np.argmax, détecter la position du maximum de la fft du signal.
Traduire cette position en fréquence.
fq = np.argmax(tf_signal)/tf_signal.size * fs
print("La fréquence du maximum de la fft du signal est {} Hz".format(int(fq)))
La fréquence du maximum de la fft du signal est 43660 Hz
2. Étude de signaux réels:¶
On retrouve bien, dans ce cas, la fréquence du son qu’on a créé. On propose de faire la même chose pour un fichier son, qui est l’enregistrement d’une note de musique. On importe un fichier son, on lit la fréquence d’échantillonnage. Le signal étant stéréo, il y a, dans le même vecteur les deux sons concaténés:
from scipy.io import wavfile
fs, data = wavfile.read('piano.wav')
print("la frequence est ",fs," donc 1 seconde du signal correspond à ",fs, " points")
print("la longueur du signal numerique mono est ",len(data)//2)
## On remet le signal à la dimension du "vrai" son.
data_stereo = np.reshape(data,(np.size(data)//2,2))
signal_mono = 0.5 * np.sum(data_stereo,axis = 1)
print("la longueur du signal mono en secondes est ", len(signal_mono)//fs)
la frequence est 44100 donc 1 seconde du signal correspond à 44100 points
la longueur du signal numerique mono est 460800
la longueur du signal mono en secondes est 20
Exercice 6:¶
Écouter le signal.
Le représenter en fonction du temps.
Représenter sa transformée de Fourier.
Que constatez-vous sur le spectre du signal ?
pl.figure(figsize = (10,10))
pl.subplot(1,2,1)
pl.plot(signal_mono)
pl.subplot(1,2,2)
tf_signal_mono = sf.fft(signal_mono)
pl.plot(sf.fftshift(tf_signal_mono))
/opt/hostedtoolcache/Python/3.7.9/x64/lib/python3.7/site-packages/numpy/core/_asarray.py:83: ComplexWarning: Casting complex values to real discards the imaginary part
return array(a, dtype, copy=False, order=order)
[<matplotlib.lines.Line2D at 0x7fc3b9caafd0>]
Le signal est composé de plusieurs fréquences : la fréquence fondamentale et les harmoniques.
FFT sur une sous partie du signal.
Exercice 7:¶
Estimer l’instant où le son apparait.
Effectuer le même travail que précédemment pour une sous partie du signal: on prendra le signal fenêtré sur une longueur de \(T\) secondes. Par exemple, \(1.5\) secondes.
Pour éviter de polluer le spectre avec les discontinuités aux bords avec les conditions périodiques de la FFT, on multiplie le signal fenêtré avec une fonction qui décroit vers \(0\) aux bords de la fenêtre. Voir ci-dessous.
## On introduit une fenetre sur laquelle on va calculer la fft.
import scipy.signal as sig
start = 44100 * 2
length = int(44100 * 0.5)
print(length)### on considère un nombre d'échantillons d'une durée équivalente de 1.5 seconde
slices = slice(start,start+length,1)
window = sig.hamming(length)
print(np.size(window))
pl.plot(window)
pl.title("Hamming window")
pl.ylabel("Amplitude")
pl.xlabel("Sample")
pl.show()
22050
22050
Exercice 8:¶
Comparer la différence entre la fft du signal fenêtré et la fft du signal fenêtré et multiplié par la fonction de Hamming introduite ci-dessus.
Exercice 9: On considère le signal ci-dessous qui est la prononciation d’une conçonne, un “s”.¶
Que pouvez-vous dire du contenu fréquentiel du signal ?
Est-ce que la localisation d’une fréquence dominante a du sens dans ce cas ?
from scipy.io import wavfile
fs, data = wavfile.read('sss.wav')
print("la frequence est ",fs," donc 1 seconde du signal correspond à ",fs, " points")
print("la longueur du signal numerique mono est ",len(data)//2)
print("la forme du signal est "+str(np.shape(data)))
## On remet le signal à la dimension du "vrai" son.
data_stereo = np.reshape(data,(np.size(data)//2,2))
signal_mono = 0.5 * np.sum(data_stereo,axis = 1)
print("la longueur du signal mono en secondes est ", len(signal_mono)//fs)
la frequence est 44100 donc 1 seconde du signal correspond à 44100 points
la longueur du signal numerique mono est 25863
la forme du signal est (51727, 2)
la longueur du signal mono en secondes est 1
from IPython.display import Audio
Audio(signal_mono, rate = 44100)
tf_signal = sf.fft(signal_mono)
pl.plot(sf.fftshift(tf_signal))
/opt/hostedtoolcache/Python/3.7.9/x64/lib/python3.7/site-packages/numpy/core/_asarray.py:83: ComplexWarning: Casting complex values to real discards the imaginary part
return array(a, dtype, copy=False, order=order)
[<matplotlib.lines.Line2D at 0x7fc3b6a26910>]
3. Le théorème de Shannon¶
## On définit un signal
def mySignal(t):
fq = 0.9
fq2 = 0.3
return np.sin(2*np.pi * fq*t) + 0.2*np.cos(2*np.pi * fq2*t)
x = np.linspace(0,2*np.pi,1000)
pl.plot(x,mySignal(x))
pl.show()
Exercice 10:¶
Choisissez un temps d’interpolation \(T_e\) et montrer la reconstruction du signal en utilisant une interpolation linéaire, c’est-à-dire des segments entre les points de la courbe.
Implémenter la formule d’interpolation de Shannon et comparer avec l’interpolation linéaire; changer de forme de signal.
N = 15
Te = np.pi / N
sampled_signal = []
x0 = []
for i in range(0, N-1):
n = i * Te
sampled_signal.append(mySignal(n))
x0.append(n)
pl.plot(x0, sampled_signal)
pl.show()
def sinc(t):
x_temp = (np.pi*t) / Te
return 1 if t == 0 else np.sin(x_temp) / x_temp
def shannon_reconstruct(t):
signal_sum = 0
for i in range(0, N-1):
nTe = i * Te
signal_sum += sampled_signal[i] * sinc(t-nTe)
return signal_sum
shannon = np.vectorize(shannon_reconstruct)
x = np.linspace(0, 2 * np.pi, N-1)
x_analog = np.linspace(0, 2*np.pi, 1000)
pl.plot(x0, sampled_signal, 'r')
pl.plot(x_analog, shannon(x_analog), 'y')
pl.figure()
pl.plot(x_analog, mySignal(x_analog) - shannon(x_analog))
pl.show()
