TFD et ses applications en audio

by Joseph Razik, on 2019-10-18
3-TFD_applications_py3

Quelques applications de la Transformée de Fourier Discrète

Pour ces exercices vous pouvez récupérer les exemples de sons "1.wav", "mystere.wav" et "etrange.wav" à l'adresse suivante http://razik.univ-tln.fr/Enseignements/doku.php?id=d32:d_32

In [1]:
%pylab inline

Populating the interactive namespace from numpy and matplotlib
In [2]:
# importation de module pour la manipulation de fichiers audios
from scipy.io.wavfile import read as wread
In [3]:
# on défini les deux fonctions vues dans l'exemple de manipulation précédent sur Fourier
def w_k(k, n, N=64):
    """ fonction élémentaire de la décomposition de Fourier """
    return exp(2j * pi * k * n / (N * 1.0))

def dft(s, N):
    """ Discrete Fourier Transform pour un signal s, sur N dimensions """
    return array([sum([ s[n] * w_k(-k, n, N) for n in range(0, N)]) for k in range(0, N) ])

Première application avec un fichier identifié

In [4]:
# lecture du fichier correspondant au son 1 DTMF du téléphone
fs_1, s_1 = wread('1.wav')
In [5]:
# fréquence d'échantillonage du signal
print("fréquence d'échantillonage : ", fs_1)
# nombre d'échantillons dans le signal
print("nombre d'échantillons : ", len(s_1))

fréquence d'échantillonage :  44100
nombre d'échantillons :  22003
In [6]:
# normalement on applique la DFT sur tout le signal donc avec N=len(s)
# cependant avec cette méthode de calcul non optimisée, c'est trop lent
# il faut sous-echantilloner, par exemple d'un facteur 8
ss_ratio = 8
sss_1 = s_1[::ss_ratio]
print("Nombre d'échantillons du signal sous-échantillone d'un facteur ", ss_ratio, " : ", len(sss_1))
# taille de la transformée de Fourier
N = len(sss_1)//2*2  # on prend une valeur paire pour N
print( "Taille de la DFT : ", N)
# Nouvelle fréquence d'échantillonage
NFS = fs_1/(ss_ratio*1.0)
print("Nouvelle Fréquence d'échantillonage : ", NFS)
# nouvelle resolution fréquentielle
# = une des valeurs de la DFT représente (englobe) combien de Hz
RF = NFS/(N*1.0)
print("Résolution fréquentielle: ", RF)

Nombre d'échantillons du signal sous-échantillone d'un facteur  8  :  2751
Taille de la DFT :  2750
Nouvelle Fréquence d'échantillonage :  5512.5
Résolution fréquentielle:  2.0045454545454544
In [7]:
# on calcul la DFT du signal sous-echantillone
X_1 = dft(sss_1, N)
In [8]:
# on affiche le spectre d'amplitude, que sur la moitié à cause de la symétrie fréquentielle
# on passe les amplitudes en dB
figure(figsize=(12,12))
plot([RF*i for i in range(0, N//2)], 10*log10(abs(X_1[0:N//2])))
xlim(0, NFS//2)
Out[8]:
(0, 2756.0)
In [9]:
# quelles sont les fréquences de forte réponse ?
m = max(10*log10(abs(X_1[0:N//2])))
pics_1 = nonzero(10*log10(abs(X_1[0:N//2])) > m*0.9)[0]
freq_1 = RF*pics_1
print(pics_1, freq_1)

[347 348 349 603] [  695.57727273   697.58181818   699.58636364  1208.74090909]
In [10]:
# Les vraies valeurs à trouver sont 697 Hz et 1 209 Hz

Manipulation avec un fichier dont le contenu est inconnu

In [11]:
# même chose avec le fichier mystère
fs_mystere , s_mystere = wread('mystere.wav')
# fréquence d'échantillonage
print("fréquence d'échantillonage : ", fs_mystere)
# nombre d'echantillon dans le signal
print("nombre d'échantillons : ", len(s_mystere))

fréquence d'échantillonage :  44100
nombre d'échantillons :  22051
In [12]:
# normalement on applique la DFT sur tout le signal donc avec N=len(s)
# cependant avec cette méthode de calcul non optimisée, c'est trop lent
# il faut sous-echantilloner, par exemple d'un facteur 8
ss_ratio_mystere = 8
sss_mystere = s_mystere[::ss_ratio_mystere]
print("Nombre d'échantillons du signal sous-échantilloné d'un facteur ", ss_ratio_mystere, " : ", len(sss_mystere))
# taille de la transformée de Fourier
N = len(sss_mystere)//2*2  # on prend une valeur paire pour N
print("Taille de la DFT : ", N)
# Nouvelle fréquence d'échantillonage
NFS = fs_mystere/(ss_ratio_mystere*1.0)
print("Nouvelle Fréquence d'échantillonage : ", NFS)
# nouvelle résolution fréquentielle
# = une des valeurs de la DFT représente (englobe) combien de Hz
RF = NFS/(N*1.0)
print("Résolution fréquentielle: ", RF)

Nombre d'échantillons du signal sous-échantilloné d'un facteur  8  :  2757
Taille de la DFT :  2756
Nouvelle Fréquence d'échantillonage :  5512.5
Résolution fréquentielle:  2.0001814223512335
In [13]:
# on calcul la DFT du signal sous-echantillone
X = dft(sss_mystere, N)
In [14]:
# on affiche le spectre d'amplitude, que sur la moitie de la symetrie frequentielle
# on passe les amplitudes en dB
figure(figsize=(12,12))
plot([RF*i for i in range(0, N//2)], 10*log10(abs(X[0:N//2])))
xlim(0, NFS//2)
Out[14]:
(0, 2756.0)
In [15]:
# quelles sont les fréquences de forte réponse ?
m = max(10*log10(abs(X[0:N//2])))
pics = nonzero(10*log10(abs(X[0:N//2])) > m*0.9)[0]
freq = RF*pics
print(pics, freq)

[385 668] [  770.06984761  1336.12119013]
In [16]:
# On en déduit que le signal correspond à la touche 5 car on trouve exactement les mêmes valeurs 770 1336

Manipulation avec un autre fichier dont le contenu est inconnu

In [17]:
# même chose avec le fichier "etrange"
fs_etrange , s_etrange = wread('etrange.wav')
# fréquence d'échantillonage
print("fréquence d'échantillonage : ", fs_etrange)
# nombre d'echantillon dans le signal
print("nombre d'échantillons : ", len(s_etrange))

fréquence d'échantillonage :  44100
nombre d'échantillons :  44055
In [18]:
# normalement on applique la DFT sur tout le signal donc avec N=len(s)
# cependant avec cette méthode de calcul non optimisée, c'est trop lent
# il faut sous-echantilloner, par exemple d'un facteur 16
ss_ratio_etrange = 16
sss_etrange = s_etrange[::ss_ratio_etrange]
print("Nombre d'échantillons du signal sous-échantilloné d'un facteur ", ss_ratio_etrange, " : ", len(sss_etrange))
# taille de la transformée de Fourier
N = len(sss_etrange)//2*2  # on prend une valeur paire pour N
print("Taille de la DFT : ", N)
# Nouvelle fréquence d'échantillonage
NFS = fs_etrange/(ss_ratio_etrange*1.0)
print("Nouvelle Fréquence d'échantillonage : ", NFS)
# nouvelle resolution fréquentielle
# = une des valeurs de la DFT représente (englobe) combien de Hz
RF = NFS/(N*1.0)
print("Résolution fréquentielle: ", RF)

Nombre d'échantillons du signal sous-échantilloné d'un facteur  16  :  2754
Taille de la DFT :  2754
Nouvelle Fréquence d'échantillonage :  2756.25
Résolution fréquentielle:  1.0008169934640523
In [19]:
# on calcul la DFT du signal sous-échantilloné
X = dft(sss_etrange, N)
In [20]:
# on affiche le spectre d'amplitude, que sur la moitié de la symétrie fréquentielle
# on passe les amplitudes en dB
figure(figsize=(12,12))
plot([RF*i for i in range(0, N//2)], 10*log10(abs(X[0:N//2])))
xlim(0, NFS//2)
Out[20]:
(0, 1378.0)
In [21]:
# quelles sont les fréquences de forte réponse ?
m = max(10*log10(abs(X[0:N//2])))
pics = nonzero(10*log10(abs(X[0:N//2])) > m*0.9)[0]
freq = RF*pics
print(pics, floor(freq))

[ 695  696  697  698  768  769  770  771  772 1205 1207 1208 1209 1211 1277
 1278 1279 1281] [  695.   696.   697.   698.   768.   769.   770.   771.   772.  1205.
  1207.  1208.  1209.  1211.  1278.  1279.  1280.  1282.]

Conclusion: ce n'est pas un signal DTMF d'une touche mais semble être la composée de deux touches

Analyse du problème de sous-échantionnage

In [22]:
# problème de sous-échantionnage trop important
# fréquence d'échantillonage
print("fréquence d'échantillonage : ", fs_etrange)
# nombre d'échantillon dans le signal
print("nombre d'échantillons : ", len(s_etrange))

fréquence d'échantillonage :  44100
nombre d'échantillons :  44055
In [23]:
# normalement on applique la DFT sur tout le signal donc avec N=len(s)
# cependant avec cette méthode de calcul non optimisée, c'est trop lent
# il faut sous-echantilloner, par exemple d'un facteur 32
ss_ratio = 32
sss = s_etrange[::ss_ratio]
print("Nombre d'échantillons du signal sous-échantillone d'un facteur ", ss_ratio, " : ", len(sss))
      # taille de la transformee de Fourier
N = len(sss)//2*2  # on prend une valeur paire pour N
print("Taille de la DFT : ", N)
# Nouvelle fréquence d'échantillonage
NFS = fs_etrange/(ss_ratio*1.0)
print("Nouvelle Fréquence d'échantillonage : ", NFS)
# nouvelle resolution fréquentielle
# = une des valeurs de la DFT représente (englobe) combien de Hz
RF = NFS/(N*1.0)
print("Résolution fréquentielle: ", RF)

Nombre d'échantillons du signal sous-échantillone d'un facteur  32  :  1377
Taille de la DFT :  1376
Nouvelle Fréquence d'échantillonage :  1378.125
Résolution fréquentielle:  1.001544331395349
In [24]:
# on calcul la DFT du signal sous-echantillone
X = dft(sss, N)
In [25]:
# on affiche le spectre d'amplitude, que sur la moitie de la symetrie frequentielle
# on passe les amplitudes en dB
figure(figsize=(12,12))
plot([RF*i for i in range(0, N)], 10*log10(abs(X[0:N])))
xlim(0, NFS//2)
Out[25]:
(0, 689.0)
In [26]:
# quelles sont les fréquences de forte réponse ?
m = max(10*log10(abs(X[0:N//2])))
pics = nonzero(10*log10(abs(X[0:N//2])) > m*0.9)[0]
freq = floor(RF*pics)
print(pics, freq)

[ 96  98  99 100 168 169 170 606 607 608 610 679 680 681 683] [  96.   98.   99.  100.  168.  169.  170.  606.  607.  608.  610.  680.
  681.  682.  684.]

Conclusion: Des pics de fréquences apparaissent à des fréquences en dehors des plages DTMF, c'est le phénomène de repliment

Utilisation des outils de calcul de la DFT optimisés: la FFT

In [27]:
# On peut utiliser la fonction optimisée fft du module numpy.fft pour accélérer les calculs
from numpy.fft import fft
In [28]:
# problème de sous-échantionnage trop important
# fréquence d'échantillonage
print("fréquence d'échantillonage : ", fs_etrange)
# nombre d'échantillon dans le signal
print("nombre d'échantillons : ", len(s_etrange))
# résolution fréquentielle
RF = fs_etrange/len(s_etrange)
print("résolution fréquentielle : ", RF)
spect = fft(s_etrange)
N = len(spect)

fréquence d'échantillonage :  44100
nombre d'échantillons :  44055
résolution fréquentielle :  1.0010214504596526
In [29]:
figure(figsize=(12,12))
plot([RF*i for i in range(0, N//2)], 10*log10(abs(spect[0:N//2])))
xlim(0, fs_etrange//2)
Out[29]:
(0, 22050)
In [30]:
# en zoomant un peu sur la partie qui nous intéresse
plot([RF*i for i in range(0, N//2)], 10*log10(abs(spect[0:N//2])))
xlim(0, 2205)
ylim(60,90)
Out[30]:
(60, 90)
In [31]:
# quelles sont les fréquences de forte réponse ?
m = max(10*log10(abs(spect[0:N//2])))
pics = nonzero(10*log10(abs(spect[0:N//2])) > m*0.9)[0]
freq = floor(RF*pics)
print(pics, freq)

[ 693  695  696  697  698  699  766  768  769  770  772 1205 1207 1208 1209
 1211 1473 1474 1475 1476 1477 1478] [  693.   695.   696.   697.   698.   699.   766.   768.   769.   770.
   772.  1206.  1208.  1209.  1210.  1212.  1474.  1475.  1476.  1477.
  1478.  1479.]

Vers l'analyse spectro-temporelle: le spectrogramme

Ceci ne correspond pas à une touche DTMF. Peut être est-ce la succession temporelle de 2 touches. Pour en le déterminer nous allons analyser par fenêtre temporelle glissant afin d'obtenir une meilleure résolution temporelle. Ceci s'appelle un spectrogramme.

In [32]:
# définition de la fonction calculant un spectrogramme sur un signal en le découpant en N fenêtres non recouvrantes.
def spectro(s, N):
    fft_len = shape(s)[0]//N//2*2
    sgram = zeros((fft_len, N))
    for k in range(N):
        lspec = 10*log10(abs(fft(s[k*fft_len:(k+1)*fft_len])))
        sgram[:, k] = lspec
        
    return sgram
In [33]:
print(type(s_etrange))
print(s_etrange.shape)
sp = spectro(s_etrange, 2)
spsh = sp.shape
print(spsh)
N = spsh[0]
# résolution fréquentielle
RF = fs_etrange/(N*1.0)
print(RF)
# durée
DUREE = s_etrange.shape[0] / (fs_etrange*1.0)
print(DUREE)

<class 'numpy.ndarray'>
(44055,)
(22026, 2)
2.002179242713157
0.9989795918367347
In [34]:
# on affiche le spectrogramme
#f, a = subplots()
#figsize(12,12)
#a.pcolormesh(array(xrange(0, spsh[1]+1)),array([RF*i for i in xrange(0, spsh[0])]), sp) 
#imshow(sp[:spsh[0]/2], aspect='auto', origin='lower')
#ylim((0,fs_etrange/2.0))
figure(figsize=(12,12))
imshow(sp, extent=[0, DUREE, 0, fs_etrange],
           interpolation='nearest', origin='lower', aspect='auto')
colorbar()
ylim((0,fs_etrange/2.0))
Out[34]:
(0, 22050.0)
In [35]:
# on zoom un peu sur la partie des fréquences qui nous intéressent
# on affiche le spectrogramme
#f, a = subplots()
#figsize(12,12)
#a.pcolormesh(array(xrange(0, spsh[1]+1)),array([RF*i for i in xrange(0, spsh[0])]), sp) 
#imshow(sp[:spsh[0]/2], aspect='auto', origin='lower')
#ylim((0,2000))
figure(figsize=(12,12))
imshow(sp, extent=[0, DUREE, 0, fs_etrange],
           interpolation='nearest', origin='lower', aspect='auto')
colorbar()
ylim((0,2000))
Out[35]:
(0, 2000)
In [36]:
sp = spectro(s_etrange, 3)
spsh = sp.shape
print(spsh)
N = spsh[0]
RF = fs_etrange/(N*1.0)
print(RF)
# on zoom un peu sur la partie des fréquences qui nous intéressent
# on affiche le spectrogramme
figure(figsize=(12,12))
imshow(sp, extent=[0, DUREE, 0, fs_etrange],
           interpolation='nearest', origin='lower', aspect='auto')
colorbar()
ylim((0,2000))

(14684, 3)
3.0032688640697356
Out[36]:
(0, 2000)
In [37]:
#figure(figsize=(12,12))
#imshow(sp[:22027/2], aspect='auto', origin='lower', interpolation='nearest', extent=[0,4, 0,22000])
#ylim((0,1000))

sp = spectro(s_etrange, 4)
spsh = sp.shape
print(spsh)
N = spsh[0]
RF = fs_etrange/(N*1.0)
print(RF)
# on zoom un peu sur la partie des fréquences qui nous intéressent
# on affiche le spectrogramme
figure(figsize=(12,12))
imshow(sp, extent=[0, DUREE, 0, fs_etrange],
           interpolation='nearest', origin='lower', aspect='auto')
colorbar()
ylim((0,2000))

(11012, 4)
4.004722121322194
Out[37]:
(0, 2000)
In [38]:
#figure(figsize=(12,12))
#imshow(sp[:22027/2], aspect='auto', origin='lower', extent=[0,100,0,500])
#ylim((0,100))

sp = spectro(s_etrange, 100)
spsh = sp.shape
print(spsh)
N = spsh[0]
RF = fs_etrange/(N*1.0)
print(RF)
# on zoom un peu sur la partie des fréquences qui nous intéressent
# on affiche le spectrogramme
figure(figsize=(12,12))
imshow(sp, extent=[0, DUREE, 0, fs_etrange],
           interpolation='nearest', origin='lower', aspect='auto')
colorbar()
ylim((0,2000))

(440, 100)
100.22727272727273
Out[38]:
(0, 2000)

On remarque que l'on est plus précis en temps (on a découpé en 100 l'échelle de temps) mais moins en fréquence (résolution fréquentielle de 100 Hz).