Dans plusieurs pages de TangenteX, j'ai abordé la TFD (Transformée de Fourier Discrète) et son algorithme de calcul le plus célèbre, la FFT. Dans le cas de l'analyse spectrale, en décrivant l'algorithme de FFT en FORTRAN ou encore en utilisant la librairie FFT de Python. J'ai par ailleurs abordé la décomposition d'un signal en séries de Fourier.
Ces pages semblent être lues car je reçois régulièrement des question à propos de la TFD, de son principe, et de la signification de ses paramètres. Alors, je me suis dit que je n'avais pas été assez clair sur le sujet et qu'il serait peut être souhaitable d'y revenir de manière plus détaillée. C'est l'objectif de cette page.
Je vais donc repartir de la notion de signal échantillonné puis étudier l'utilisation de l'algorithme de TFD, d'abord par un calcul direct sur un exemple simple puis sur un signal réél à l'aide de la librairie FFT de Python.
Pour aborder le processus d'échantillonnage, je vais partir d'un signal très simple, une sinusoïde \( s(t) = a\cos(\omega t)\), de pulsation \( \omega \), de phase nulle et d'amplitude maximum a. J'aurais pu indifférement choisir un sinus... La fonction s(t), notre signal dirons les physiciens, est une fonction continue, dérivable et intégrable sur \( \mathbb{R}\).
Mon objectif est de connaitre le spectre de ce signal, c'est à dire la somme des différents signaux élémentaires en sin et cos qui le constituent, puisque nous savons depuis Fourier et ses successeurs qu'un signal peut être, sous certaines conditions satisfaites ici, décomposé en une somme infinie de sin et cos. Ce qui m'intéresse, c'est de connaître la fréquence et l'amplitude de chacun de ces signaux élementaires, données qui forment le spectre du signal s(t).
Les mathématiques nous offrent un outil fantastique pour calculer le spectre d'une fonction intégrable comme notre signal s(t), c'est la transformée de Fourier. La transformée de Fourier \( \hat{f} \) d'une fonction f(t) est définie par \( \forall \nu \in \mathbb{R} \: \hat{f} = \int_{-\infty}^{+\infty} f(t) \exp(-i2 \pi \nu t) dt \).
Dans le cas de notre signal élémentaire, il est donc possible, par un simple calcul d'intégrale, de calculer la transformée de Fourier et donc le spectre de notre signal s(t). Vous noterez toutefois attentivement les hypothèses de calcul:
Ce qui nous pose au moins deux problèmes :
Dans la suite de cette page, nous allons voir comment traiter ces deux problèmes, et quelques autres..
Imagnons que nous voulions calculer le spectre d'un signal audio, le son produit par un instrument de musique par exemple. Il semble difficile d'en calculer la transformée de Fourier à la main et donc nous allons vouloir utiliser un ordinateur ou un analyseur de spectres, ce qui est pratiquement la même chose.
Il va d'abord falloir acquérir ce signal, c'est à dire le transformer en un signal traitable par un ordinateur, généralement un signal électrique. Ce sera le rôle du microphone, qui transformera le signal de variation de pression de l'air (le son) en un signal électrique, une variation de tension en fonction du temps. Nous appelerons ce signal électrique s(t), son unité étant le Volt donc. En passant, le microphone va ajouter à notre signal acoustique quelques bruits, acoustiques dus à l'environnement et électriques dus au transducteur acoustique/électrique.
Notre signal électrique s(t) est continu en première approche. Pour le rendre compréhensible par un outil numérique, il faut l'échantillonner, c'est à dire mesurer à intervalle constant la valeur de son amplitude. Nous allons transformer notre signal continue en une suite de valeurs numériques discrètes. L'intervalle de mesure est la période d'échantillonnage Te.
Imaginons que je fasse circuler mon signal électrique dans un fil muni d'un interrupteur. Lorsque l'interrupteur est fermé, la tension à ses bornes est nulle (il n'offre aucune résistance). J'ouvre maintenant mon interrupteur et je mesure immédiatement la tension à ses bornes. Je pourrai connaître ainsi l'amplitude de mon signal au moment même de l'ouverture de l'interrupteur.
Imaginons maintenant que j'invente un dispositif qui ouvre à intervalle constant Te, pendant un temps aussi bref que possible, l'interrupteur et que je puisse mesurer dans le même temps très bref de cette ouverture la tension aux bornes de l'interrupteur. J'obtiendrai une suite de mesures de l'amplitude de mon signal. Cette opération est effectuée par un dispositif électronique qui s'appelle un échantillonneur.
Mathématiquement, l'opération d'échantillonnage se traduit par un produit (une convolution plus exactement) de deux fonctions : notre signal s(t) et une fonction que l'on appelle chez les physiciens un "peigne de Dirac". Voyons de quoi il s'agit.
La suite numérique qui représente mon signal s(t), que je noterai S(t) est composée des valeurs de s(t) pour t = nTe, n étant borné par la durée de mon signal et de l'opération d'échantillonnage, d'où S(t) = s(nTe). Plus précisement, je peux écrire S(t) = s(t) * P(t), avec P(t) la fonction "peigne de Dirac".
Le peigne de Dirac est une fonction très particulière, dont les propriétés sont:
On l'a note classiquement \( P(t) = \sum_{n=-\infty}^\infty \delta (t - nTe)\).
Notre signal échantillonné s'écrit donc \( S(t) = \sum_{n=-\infty}^\infty s(nTe) \delta (t - nTe)\) . Nous n'avons plus un signal continu, mais un signal discret !
Une dernière chose concernant le peigne de Dirac, qui va nous servir dans la suite, il s'agit de sa décomposition en séries de Fourier. Je ne vais pas rentrer dans les détails du calcul, en vous la donnant: \( P(t) = \frac{t_0}{T_e} + 2\frac{t_0}{T_e} \cos (\omega_e t) + 2\frac{t_0}{T_e} \cos (2\omega_e t) + 2\frac{t_0}{T_e} \cos (3\omega_e t) + ...\).
Je suppose dans ce calcul que la durée d'ouverture de l'interrupteur to est très petite devant la période Te d'échantillonnage, ce qui est toujours vrai dans les applications pratiques.
Pour que ce signal soit exploitable, il reste une dernière chose à faire, c'est de quantifier pour numériser l'amplitude de chaque élément de notre signal échantillonné. C'est le travail d'un dispositif que l'on nomme CAN : Convertisseur Analogique Numérique.Il va transformer en nombres chaque valeur analogique de tension.
La taille de ces nombres, en bits, va déterminer la précision de notre CAN. Imaginons que nous voulions quantifier un signal dont l'amplitude varie entre -V et +V volts. Découpons l'intervalle [-V, +V] en p intervalles. L'amplitude v de chaque élement du signal s'exprimera donc sous la forme \( v = k*2V/2^p \).
En pratique, la valeur V est une caractéristique technique du CAN. Les CAN courants supportent des tensions d'entrée comprises entre -1V et +1V ou -5V et +5V. p désigne la taille des mots en bits générés par le CAN, sont c'est aussi une caractéristique. Les CAN les plus courants utilisent des mots de 8, 16, 24 ou 32 bits, signés ou non signés. Dans le cas de mots signés, un mot de p bits permet de coder (p-1) valeurs de tension. Pour fixer les idées, le CAN d'un iPhone utilise des mots de 16 bits signés.
Appliquons ce que nous venons de voir à notre signal sinusoïdal \( s(t) = a\cos(\omega t)\). Le signal échantillonné S(t) est obtenu par le produit de s(t) et de P(t), ou ce qui revient au même à la limite par le produit de leur décomposition de Fourier. Pour simplifier les calculs, nous nous limiterons à l'ordre 2, ce qui nous donne:
\( S(t) = (a\cos(\omega t))(\frac{t_0}{T_e} + 2\frac{t_0}{T_e} \cos (\omega_e t) + 2\frac{t_0}{T_e} \cos (2\omega_e t))\)
soit, en développant, en factorisant :
\( S(t) = a\frac{t_0}{T_e}\cos(\omega t) + 2a\frac{t_0}{T_e} \cos(\omega_e t)\cos(\omega t) + 2a\frac{t_0}{T_e} \cos(2\omega_e t)\cos(\omega t) \)
En appliquant la formule trigo de linéarisation \(\cos a \cos b = \frac{1}{2}(\cos(a + b) + \cos(a - b)) \), j'obtiens finalement :
\( S(t) = a\frac{t_0}{T_e}\cos(\omega t) + a\frac{t_0}{T_e}(\cos(\omega_e + \omega)t + \cos(\omega_e t - \omega)t) + a\frac{t_0}{T_e}(\cos(2\omega_e + \omega)t + \cos(2\omega_e t - \omega)t) \)
Pour que cette équation soit plus explicite, puisqu'il s'agit de calculer un spectre de fréquences, je vais employer la relation \( \omega = 2 \pi f \) pour l'exprimer avec des fréquences, ce qui nous donne :
\( S(t) = a\frac{t_0}{T_e}\cos(2 \pi f t) + a\frac{t_0}{T_e}(\cos(2 \pi (F_e + f))t + \cos(2 \pi (F_e - f)t) + a\frac{t_0}{T_e}(\cos(4 \pi (F_e + f))t + \cos(4 \pi (F_e - f))t) \)
En observant ce résultat, nous remarquons plusieurs propriétés importantes du spectre d'un signal échantillonné :
Ces observations sont généralisables à n'importe quel signal. Cela signifie que l'échantillonnage d'un signal provoque une reproduction de son spectre pour chaque multiple de la fréquence d'échantillonnage.
Imaginons que je procède à l'échantillonnage d'un signal composite qui contient une raie de fréquence Fmax. Son spectre va donc s'étaler entre nFe - Fmax et nFe + Fmax. Si Fe <= Fmax, je vais avoir un problème : les spectres vont se superposer et les résultats seront inexploitables ! Il existe donc une relation de contrainte entre Fe et Fmax. Il s'agit du célèbre théorème de Shannon, dont j'ai parlé plusieurs fois dans les pages précédentes, qui pose que Fe > 2*Fmax. Evidemment, cela implique que l'on connaisse Fmax, ce qui n'est toujours pas évident...
Le cas typique de problème lié à la connaissance de Fmax est la présence de bruit de large spectre. Prenons l'exemple d'un enregistrement vocale à l'aide d'un micro. La voix évolue dans un domaine de fréquences compris entre 10 Hz et 20 kHz. L'échantillonnage se fait généralement à 44 kHz, c'est à dire à plus de deux fois la fréquence maximum attendue. Parfait, sauf si votre micro produit du bruit à 25 ou 30 kHz, voire plus. Dans ce cas, votre choix de fréquence d'échantillonnage n'est plus bon. Pour éviter ce genre de mésaventure, on place à l'entrée de l'échantillonneur un filtre coupe-haut dit anti-repliement pour éliminer toute la partie du signal dont la fréquence est plus grande que Fmax attendue.
Lorsque dans un programme Python, je veux tracer une courbe, un cosinus en fonction du temps par exemple, je discrètise forcément le temps. Je le fais par exemple en utilisant l'instruction :
t = arange(t0,tf,pas)
avec laquelle je construis une suite entre t0 et tf, qui contient des valeurs discrètes de temps avec un intervalle de "pas". J'obtiens la courbe ci-dessous, qui vous semble continue, mais que ne l'est qu'artificiellement parce que le programme relie les points entre eux.
En réalité, il serait préférable de l'imaginer comme ça :
Par construction, chaque fois que nous utiliserons un vecteur temps ou autre dans un programme, il faudra penser qu'il s'agit d'un signal discret.
L'échantillonnage, qui revient à discrétiser le temps lors de la numérisation du signal, a des répercussions essentielles sur son spectre, dont il faut tenir compte pour déterminer la fréquence d'échantillonnage.
J'ai rappelé en début de page la définition de la transformée de Fourier : \( \hat{f} = \int_{-\infty}^{+\infty} f(t) \exp(-i2 \pi \nu t) dt \) et noté que celle-ci s'appliquait à des fonctions réelles continues. Or nous venons de voir que nous pouvions travailler avec un signal continu, puisque informatiquement ça n'existe pas.
Les physiciens ont donc été amenés à définir une nouvelle transformée de Fourier, qui s'appliquerait à des signaux discrets, la TFTD. Il s'agit de passer d'une intégrale entre plus et moins l'infini dans le domaine temps à une série comprise entre les mêmes bornes, et de préférence absolument convergente. La définition de la TFTD découle naturellement de cette observation:
\( S(f) = \sum_{-\infty}^{+\infty} f(t) \exp(-i2 \pi f t) \)
Appliquons cette définition à notre signal échantillonné \( s^*(t) = \sum_{n=-\infty}^\infty s(nTe) \delta (t - nTe)\). Après quelques calculs que je vous épargne, j'obtiens la définition suivante pour la TFTD d'un signal échantillonné:
\( S(f) = \sum_{n = -\infty}^{+\infty} S(nT_e) \exp(-i2 \pi n f T_e) \)
Par convention, dans les calculs, on normalise la période d'échantillonnage en posant Te = 1. Cela ne change rien au principe et simplifie les calculs. Après normalisation de Te, nous avons donc :
\( S(f) = \sum_{n = -\infty}^{+\infty} S(n) \exp(-i2 \pi n f) \)
Notons deux points :
Nous verrons ci-dessous comment résoudre ces problèmes. Mais d'abord, examinons quelques propriétés de la TFTD:
Tout ça est bien beau, mais pour calculer une TFTD il faut une suite infinie d'échantillons, ce qui est impossible. Et encore plus fort, nous utilisons pour le calcul une fonction de fréquence continue, alors que nous ne savons traiter que des nombres discrets. Voyons comment traiter ce problème.
Après avoir discrétiser le temps, il nous reste à discrétiser la fréquence et à traiter un signal de durée finie, c'est à dire avec un nombre fini d'échantillons. Alors imaginons que notre signal discret soit formé d'une suite de N échantillons {x(0), x(1),...,x(N-1)}. Au lieu de chercher à décomposer ce signal sur un intervalle continu de fréquences, comme dans le cas d'une TFTD, je vais choisir une suite discrète de fréquences de L éléments. Je pose donc \( f = k \Delta F \) avec \( \Delta F = \frac{1}{L}\) et k variant de 0 à L-1.
La définition de la TFD, la transformée de Fourier discrète devient alors:
\( S(\frac{k}{L}) = \sum_{n = 0}^{N-1} S(n) \exp(-i2 \pi \frac{k}{L} n) \)
Dans cette formule, n est l'indice de boucle temporelle, en fait l'indice sur le signal échantillonné, k est l'indice des fréquences. Dans la pratique, nous verrons plus loin que l'on pose la plupart du temps L = N.
Nous voilà en possession d'un outil, la TFD, qui nous permet d'analyser la décomposition fréquentielle d'un signal échantillonné. Voyons ce que cela peut donner en calcul numérique.
Pour illustrer nos petites expériences sur la TFD, je vais construire un signal échantillonné, qui sera formé par la somme de deux sinusoïdes d'amplitude et de fréquence différentes.
La fréquence maximum de ce signal sera de 20 Hz et je souhaite une résolution spectrale de 0.1 Hz. La fréquence d'échantillonnage Fe sera donc de 40 Hz minimum, pour respecter le critère de Nyquist. Dans le code du script TFDCalculBrut.py, je déclare ces informations par :
Fe = 40.0 # fréquence d'échantillonnage
dF = 0.1 Hz # résolution en fréquence de l'analyse spectrale
Le nombre de mesures de mon échantillon se déduit de ces deux informations :
N = int(Fe/dF) # nombre d'échantillons
Comme indiqué ci-dessus, je choisis le même nombre de mesures pour la discrétisation du temps et de la fréquence (hypothèse N = L)
Je peux donc calculer la durée de l'échantillonnage qui est donnée par Te = N/Fe d'où Te = 1/dF :
Te = 1.0/dF # durée du signal à échantillonner
Muni de ces informations, je peux maintenant construire mon signal discret :
t = arange(0.0,Te,Te/N) # vecteur temps discrétisé
A0 = 5.0; F0 = 10.0 # amplitude et fréquence de la sinusoïde 0
A1 = 2.0; F1 = 5.0 # amplitude et fréquence de la sinusoïde 1
s = A0*sin(2.0*pi*F0*t) + A1*sin(2.0*pi*F1*t)
Je définis également le vecteur de fréquences discrètisées, avec le même nombre de mesures que pour le signal :
f = dF*arange(N)
Le script TFDCalculBrut.py propose une technique de calcul brut, qui applique la formule de définition de la TFD. En fait, j'ai écris deux scripts, l'un qui fait un calcul avec une double boucle sur k et n, comme dans la formule (TFDCalculBrut.py) et l'autre qui utilise les capacités matricielles de Python (TFDCalculBrutVectorise.py). Ce dernier est basé sur la transformation matricielle de la définition de la TFD. Vous pourrez contrôler que ces deux scripts donnent le même résulat, mais dans des temps très différents ! Sur mon Mac, la double boucle calcule la TFD sur le signal s en 1 seconde environ, et la forme matricielle en 0,01 seconde environ ! La différence se passe de commentaire !
Le code de calcul du spectre avec une double boucle :
X = zeros(N,dtype= complex) # définition du vecteur spectral
ex = -2j*pi/N # part invariable de l'exposant
for k in range(0,N-1):
y = 0.0
for n in range(0,N-1):
y += s[n]*exp(ex*n*k)
X[k] = y
Et maintenant le code matriciel. Je commence par construire la matrice M des coefficients exponentiels :
n = arange(N)
k = n.reshape((N,1))
M = exp(-2j*pi*k*n/N)
Puis j'effectue le produit matriciel X = M.s :
X = dot(M,s)
Et je poursuis le calcul comme pour l'algorithme "double boucle". Avouez que c'est plus élégant et beaucoup plus rapide !
Je vais calculer le spectre du signal s en calculant l'amplitude de chaque composante de X, qui est, je vous le rappelle une matrice à coefficients complexes. La fonction abs() appliquée à une matrice complexe calcule le module de chaque élément complexe de la matrice.
FacteurNorm = 2.0/N
SP = FacteurNorm * abs(X)
Le code d'affichage du spectre SP n'appelle pas de commentaire.
Par les deux méthodes, ce que vous vérifierez, j'obtiens le spectre suivant :
Sur le spectre on aperçoit deux couples de raies, qui sont réparties symétriquement par rapport à 20 Hz, soit Fe/2. Ceci est conforme à la théorie que nous avons vu plus haut. Ainsi, si l'on veut afficher un spectre simple, il faut n'afficher que les composantes du spectre dans le domaine fréquentiel [0, Fe/2], ce qui est classiquement fait dans les programmes de calcul de TFD et que nous ferons plus loin.
Dans chaque couple, nous observons une raie à 10 Hz d'amplitude 5 unités et unes raie à 5 Hz d'amplitude 2 unités. Ce sont bien les valeurs attendues de construction de notre signal.
Vous noterez que les raies sont très étroites, parce que la résolution spectrale choisie est assez fine (dF = 0.1 Hz). Vous pouvez vous amuser à la modifier, pour voir les raies s'empâter ou bien au contraire s'affiner en fonction de la résolution en fréquence que vous choisirez.
A propos des unités de chaque axe :
Pour ceux qui voudraient afficher l'amplitude en dB, ils peuvent dans l'instruction plot() remplacer la variable SP par l'expression :
20*log10(SP/SP.max())
Voici le spectre obtenu :
En résumé, la méthode directe donne plutôt de bons résultats, surtout si l'on utilise la vectorisation permise par Python. L'utilisation d'une double boucle est à déconseiller formellement. Quoiqu'il en soit, cette méthode de calcul nécessite N^2 opérations pour calculer un TFD, ce qui est une complexité élevée.
Dans le calcul du spectre, vous avez sans doute remarqué que je multiplie le module de chaque composante par un facteur de normalisation égal à 2/N. La question m'est souvent posée de savoir pourquoi, d'autant que ce facteur peut varier selon les conventions de définition de la TFD. Les explications qui suivent valent donc pour la définition que je vous ai donné de la TFD.
Le facteur 2 - Il faut considérer deux choses : chaque élément du vecteur X[k] est un complexe ; et chaque élement du spectre affiché SP[k] est un réel, égal au module de X[k]. Si vous regardez à nouveau la formule de la TFD, vous noterez, en vous souvenant de l'expression d'un sinus sous forme complexe de la formule d'Euler, qu'il s'agit d'une demi-somme de sinus (seule la partie exp(-ix) est présente). Et vous vous souvenez aussi que dans la formule d'Euler, l'argument est divisé par 2. Ici, l'argument est l'amplitude de notre signal. Aussi, pour retrouver l'amplitude initiale, je multiplie par 2 l'amplitude calculée.
La division par N - Observez le code de calcul de X[k] : on l'obtient par la somme de 0 à N-1 des éléments s[n]. Or s[n] est de la forme A*sin(nTe), c'est à dire que X[k] est affecté d'un coefficient A*N, A étant l'amplitude du signal. Pour obtenir l'amplitude correcte lors du calcul du spectre, je divise l'amplitude totale par N.
Si vous comptez le nombre d'opérations nécessaires pour calculer le spectre d'un signal par la méthode brute, matricielle ou pas, vous trouverez N^2 multiplications et N(N-1) additions, soit un ordre de grandeur de N^2 opérations, et je vous rappelle qu'il s'agit d'opérations sur des complexes. On dit que l'algorithme brut est de complexité O(N^2). Dans la majorité des applications, la TFD est calculée sur des signaux échantillonnés de grande taille (plusieurs milliers de points) ce qui pose un réel problème de performances ! Surtout dans les années 60, avec les moyens de calcul du moment, c'est à dire des ordinateurs bien moins puissants que vos tablettes !
Aussi, en 1965, Cooley et Tuckey ont proposé un algorithme révolutionnaire qui calculait une TFD avec une complexité O(Nlog2(N)), soit un gain en nombre d'opérations de 3 ordres de grandeur (supérieur à 1000), et donc autant en performances. C'est le fameux algorithme dit "FFT" ou Transformée de Fourier Rapide.
Je ne vais pas le décrire ici. Ceux qui sont intéressés par plus de détails sur cet algorithme peuvent consulter le "Méthodes numériques appliquées pour le scientifique et l'ingénieur" de Jean-Philippe Grivet au chapitre 9 et bien les "Numerical Recipes" en C ou Fortran, deux ouvrages, entres autres, qui décrivent la technique de cet algorithme.
J'en profite pour aborder un point qui fait aussi l'objet de nombreuses questions. J'ai écris dans les pages précédentes qu'il convenait que le nombre d'échantillons du signal à traiter par FFT soit une puissance de 2. Ce n'est pas obligatoire, mais ceux qui ont pris le temps de comprendre le fonctionnement de l'algorithme auront compris que cet algo fonctionne de façon optimum lorsque cette règle est respectée. Le pire, c'est à dire un algo qui se traîne, est observé quand N est premier. Donc, dans la mesure du possible, respectez la règle N = 2^p . D'ailleurs, la documentation de Python sur la fonction fft() est assez claire sur le sujet. J'en donnerai un exemple plus loin.
La fonction fft() est disponible dans le module scipy.fftpack, que vous devrez importer pour pouvoir utiliser la fonction fft(). Elle comporte plusieurs paramètres :
En général, vous appellerez fft(x) ou parfois fft(x,n) si vous voulez optimiser votre calcul en déclarant n = 2^p. Attention, le remplissage par des zéros est une opération qui doit être bien réfléchie, au risque de donner des résultats surprenants.
Il est important de connaitre le format du vecteur de retour de fft(). Posons X = fft(x). La structure de X est particulière et explique certaines manipulations faites dans le code et qui posent parfois question.
Il faut savoir que fft() calcule le spectre du signal entre [-Fe/2, Fe/2]. X[0] contient le terme continu du signal. X[1:n/2] contient les termes de fréquences positives et X[n/2:n] contient les termes de fréquences négatives par ordre décroissant des fréquences.
Comme assez généralement c'est le spectre du signal entre 0 et Fe/2 qui nous intéresse, cela explique les lignes de code suivantes :
SP = FacteurNorm*abs(Xc[1:N/2])
Je récupère la partie comprise entre 1 et N/2, qui contient les fréquences positives par l'instruction Xc[1:N/2].
De même pour le vecteur fréquences dans cette instruction :
f = fftfreq(N,1.0/Fe)[1:N/2]
Le code est dans le script TFDCalculFFT.py. La définition du signal reste inchangée. Seule varie la fonction de calcul du spectre, qui est réduite à sa plus simple expression :
Xc = fft(s)
Puis la normalisation du spectre :
SP = FacteurNorm*abs(Xc[1:N/2])
avec le facteur de normalisation habituel. Notez que je calcule le spectre que pour la partie positive des fréquences.
Enfin, je construis le vecteur d'échantillons des fréquences pour afficher le spectre en utilisant la fonction fftfreq(), toujours en ne conservant que les fréquences positives comprises entre 0 et Fe/2.
f = fftfreq(N,1.0/Fe)[1:N/2]
Le code d'affichage du spectre est identique aux autres scripts.
Le spectre obtenu est identique à celui obtenu par le calcul brut, si ce n'est que je n'ai retenu que la partie du spectre comprise entre 0 et Fe/2.
![Spectre FFT d'un signal composite - [0,Fe/2]](images/TFDFFT1.png "Spectre FFT d'un signal composite - [0,Fe/2]")
Le temps de calcul est insignifiant, quelques millisecondes. Sur des processeurs spécialisés, les DSP (Digital Signal Processor), on calcule une TFD par l'algorithme FFT en quelques microsecondes, ce qui permet de calculer des spectres en temps réel sur des fenêtre assez larges.
Dans le script TFDCalculFFTShift.py, je vous propose une légère variante. En appliquant la fonction fftshift() à la fois au vecteur Xc résultant de fft() et au vecteur f des fréquences, je réorganise le vecteur Xc pour mettre les fréquences négatives et positives dans l'ordre croissant de [-Fe/2, Fe/2]. Cela nous permet de visualiser la symétrie du spectre, comme le montre la courbe ci-dessous.
![Spectre FFT d'un signal composite - [-Fe/2,Fe/2]](images/TFDFFT2.png "Spectre FFT d'un signal composite - [-Fe/2,Fe/2]")
Nous allons maitenant tester notre outil FFT pour visualiser le spectre d'un signal réel, l'enregistrement du bruit des vagues, enregistré avec un iPhone au Tréport.
Ce signal est enregistré dans un fichier au format wav, dont vous pourrez trouver le format sur wikipedia. Ce qu'il faut retenir, c'est que les paramètres d'enregistrement (fréquence d'échantillonnage, nombre de canaux) définissent les paramètres de la fonction de lecture des données dans le fichier. C'est ainsi que le script proposé est programmé pour lire un signal monophonique échantillonné à 8 kHz. Le signal électrique, la tension plus précisément, est numérisé sur des mots de 16 bits signés.
Du fichier original BruitWav.wav, j'ai extrait le bruit d'une vague, que j'ai enregistré dans le fichier BruitWav1.wav. La durée de l'enregistrement est de 5 secondes, très exactement 4,996 secondes. Attention, le fichier BruitWav1.wav doit être stocké dans le même répertoire que le script TFDSpectreSignalReel.py qui va nous permettre de calculer son spectre.
Une précision qui a son importance dans notre code : les fichiers wav ne contiennent pas les valeurs d'amplitude réelle mais les coefficients de quantification (le coefficient k que j'ai employé plus haut). Pour tracer le signal avec son amplitude réelle, et pas des amplitudes de 10 000 ou 15 000, il faut diviser ces coefficients par 2^(p-1). Dans le cas de l'iPhone, p = 16.
Je commence par ouvrir le fichier qui contient le signal. La fonction read() me retourne la fréquence d'échantillonnage Fe et un vecteur contenant le signal. ATTENTION : ce code correspond à la lecture d'un canal unique, un signal monophonique. Si vour utilisez un signal stéréo, il faudra adapter le code. Je récupère aussi le nombre de points de mon signal avec la méthode shape. Puis je calcule Te.
Fe,signal = read('BruitWav1.wav')
N = signal.shape[0] # nb de points dans l'échantillon
Te = 1.0/Fe # période d'échantillonnage
Puis je calcule l'amplitude réelle:
FacteurQuantification = 1.0/2**15
signal = FacteurQuantification*signal
Je construis mon vecteur temps t, le plus simplement du monde :
t = arange(0.0,N,1.0)*Te
puis je lance ma FFT. Le paramètre N est ici facultatif, mais je l'indique pour vérifier tout à l'heure l'influence de N sur les performances :
Xc = fft(signal,N)
je construis mon spectre normalisé comme ci-dessus :
FacteurNorm = 2.0/N
SP = FacteurNorm*abs(Xc[1:N/2])
je construis mon vecteur fréquence en utilisant fftfreq(), en lui allouant la même longueur que le vecteur SP (sinon gare aux erreurs dans le plot !) :
f = fftfreq(N,Te)[1:N/2]
et j'affiche le signal et son spectre comme d'habitude !
Voici la forme du signal lu dans BruitWav1.wav et son spectre tracé avec le code ci-dessus.
Je vous laisse découvrir ce spectre, qui est ramassé dans les basses fréquences avec quelques raies dominantes dont vous chercherez la signification (je ne la connait pas).
Sur mon Dell Precision M6400 (qui n'est plus un foudre de guerre !), le calcul de la FFT sur les 39965 points que comptent mon signal dure environ 0,29 secondes. Maintenant, fixons un nombre de points un peu inférieur, mais égal à une puissance de 2. Le plus proche est bien sur N = 32768 (2^15). Les instructions correspondantes sont :
N = 32768
Xc = fft(s,N)
Le spectre obtenu est :
Le signal a été tronqué à 32768 points. le spectre est presque quasi identique, mais le temps de calcul indiqué par Spyder est de 0.0 seconde : insignifiant ! Le fait d'utiliser un vecteur signal de taille égale à une puissance de 2 a une influence immédiate sur les performances ! Pensez-y...
J'ai laissé de côté un point abordé plus haut, un tout petit problème... Lorsque je construis un signal simulé, ou que j'acquière un signal avec un appreil numérique, non seulement j'échantillonne le signal, en créant une discrétisation temporelle, mais je le tronque : sa durée est finie. Je ne sais pas calculer un spectre d'un signal infini, parce qu'en physique expérimentale, un signal infini n'existe pas. Le problème tient dans le fait que la TFD n'est qu'une approximation de la transformée de Fourier et que son application sur un signal de durée finie a un impact plus ou moins négligeable sur le spectre. Les raies du spectre sont élargies et des raies parasites apparaissent.
Il existe un outil pour corriger plus ou moins complétement ce défaut : c'est le fenêtrage. Nous l'aborderons dans une prochaine page.
En analyse spectrale pratique, il est souvent nécessaire de distinguer des raies qui sont très proches, d'amplitude comparable ou bien très différente. Le physicien cherchera donc à améliorer la précision de sa TFD. Il s'agit d'améliorer la résolution fréquentielle pour pouvoir discriminer deux raies proches, et pour cela diminuer la largeur des raies (on dit aussi le lobe des raies). il s'agit aussi augmenter la résolution en amplitude afin de discriminer deux raies d'amplitude proche.
Plusieurs moyens sont possibles :
Et aussi adopter une autre fenêtre que la porte rectangulaire pour tronquer le signal. Nous verrons cela dans la page consacrée au fenêtrage.
Il existe des appareils numériques, les analyseurs de spectres, qui effectuent l'analyse spectrale d'un signal continu en "temps réel", c'est à dire à mesure de la production du signal. Nous les trouverons dans les laboratoires d'électronique pour mesurer la qualité d'un filtre, d'un émetteur ou d'un récepteur, dans les labo de mesures physique et dans presque tous les labo de physique.
Les analyseurs de spectre doivent leur existence à deux éléments : l'algorithme de FFT et les DSP, des processeurs spécialisés ultra-rapides. Ils procédent en continu sur le signal entrant à :
La fréquence du signal que peut traiter un analyseur de spectres, et la richesse et la précision du spectre produit dépendra bien sur de la puissance du DSP dont il dispose. Il existe aujourd'hui des analyseurs de spectres, qui par des moyens détournés, arrivent à tracer des spectres de signaux de plusieurs dizaines de Ghz. Voilà à quoi ressemble ce genre d'appareil (ici un analyseur Agilent - Ex HP, de 13 GHz) :
Il existe aussi des logiciels qui font ce genre de traitement, sur des signaux de fréquences beaucoup moins élevées. Ils procèdent de façon identique. Inutile de dire que vous avez intérêt à posséder un processeur puissant pour les faire tourner. Dans certains PC, on peut ajouter des co-processeurs ou des cartes spécialisées pour faire du traitement du signal.
Les scripts Python étudiés dans cette page sont disponibles dans le package TFDPython.zip :
J'espère que cette page vous aura éclairé un peu plus sur la TFD et l'algorithme de FFT, et répondu aux éventuelles questions que vous vous êtes posé en lisant mes pages précédentes.
Je reviendrai prochainement sur le fenêtrage temporel des signaux. N'hésitez pas à continuer à me poser vos questions par mail si vous ne trouvez pas votre réponse ici.
Contenu et design par Dominique Lefebvre - tangenteX.com avril 2016 -- Vous pouvez me joindre par ou sur PhysiqueX
Cette œuvre est mise à disposition selon les termes de la Licence Creative Commons Attribution - Pas d’Utilisation Commerciale - Pas de Modification 3.0 France.