Retour

Transformée de Fourier Discrète et Python


La transformée de Fourier

Nous avons déjà abordé dans une autre page l'analyse de Fourier, c'est à dire la décomposition d'une fonction continue périodique en sommes de sinus et de cosinus.
Nous avons vu que, munis des théorèmes de Fourier et de Parseval, nous pouvions calculer les coefficients de Fourier et calculer le spectre d'une fonction, continue et périodique.
Tout cela est fort bien, mais dans la pratique, il est plutôt rare de disposer d'un signal continu et périodique ! Prenons le cas par exemple d'un enregistrement sonore réalisé à l'aide d'un smartphone et dont vous voudriez établir le spectre. Ce signal n'est pas continu: c'est une suite de mesures de tension, dans le nombre dépend de la fréquence d'échantillonnage de votre smartphone. Il n'est pas périodique non plus. Il possède un début et une fin, le début et la fin de l'enregistrement. Il ne possède donc aucune des caractéristiques prérequises pour utiliser l'analyse de Fourier.
Et pourtant, vous savez que l'on peut facilement tracer le spectre d'un signal sonore enregistré par un smartphone. C'est possible grâce à un outil mathématique merveilleux, la transformée de Fourier et à sa fille la Transformée de Fourier Discrète (la TFD).

Il n'est pas question ici de faire une introduction mathématique à la transformée de Fourier. Sachez simplement qu'elle permet d'étendre aux fonctions non périodiques l'analyse de Fourier, ce qui permet donc de calculer le spectre d'une fonction non périodique.
Il reste que dans sa version originale, une transformée de Fourier ne peut s'appliquer qu'à une fonction intégrable, ou au mieux de carré sommable grâce au théorème de Plancherel. Heureusement, des mathématiciens, en particulier Laurent Schwartz, ont travaillé pour qu'elle soit utilisable sur des suites discrètes périodiques ou rendues périodiques.
Bien que cela ne soit pas le sujet, il faut savoir que la transformée de Fourier sert à bien d'autres choses que l'analyse harmonique et que c'est toujours un domaine de recherche vivace en mathématique et physique mathématique.
Une précision de langage : lorsque les mathématiciens parlent de fonction, les physiciens parlent de signal. Vous aurez noté que j'emploie les deux indifféremment. De même, ce que les physiciens appellent spectre correspond au résultat de la transformée de Fourier d'une fonction pour les mathématiciens.

De ces études est né l'outil qui nous intéresse aujourd'hui: la TFD. Elle permet de faire de l'analyse harmonique sur des signaux discrets, c'est à dire non continus. Et il ne vous aura pas échappé que les ordinateurs et autres instruments de mesures numériques ne savent produire et manipuler que des signaux discrets!

La Transformée de Fourier Discrète

En physique numérique, nous disposons presque toujours de signaux issus d'une acquisition électronique ou de résultats de mesures discrètes. Ces signaux et résultats ne sont pas infinis, et généralement sont non périodiques. Pas vraiment le domaine d'application de la transformée de Fourier !
La TFD a été conçue pour traiter ce genre de données. Elle réalise une décomposition d'un signal tronqué (fini, de durée T) et échantillonné (discret) en une série de Fourier, en le "périodisant", avec une période égale à T.

La TFD nous permet donc de calculer de façon approchée les coefficients de Fourier des différentes harmoniques d'un signal presque quelconque et ainsi d'obtenir son spectre. C'est dire son utilité en physique et en ingénierie.

Pour l'approche mathématique de la TFD (DFT en anglais), je vous renvoie à l'article Wiki correspondant et aussi à celui du manuel de fft de scipy, qui en présentent les principaux aspects.

L'algorithme de FFT

Ne pas confondre FFT et TFD ! C'est souvent le cas par abus de langage. La FFT (Fast Fourier Transform) est un algorithme alors que la TFD est un outil mathématique.
J'ai déjà abordé la FFT sur une autre page, je ne vais donc pas y revenir, en vous invitant à la consulter.

J'aimerais cependant revenir sur quelques précautions à prendre lors de l'usage de l'algorithme de FFT:

Nous aurons l'occasion de revenir sur ce dernier point, le fenêtrage, dans une autre page. C'est un des problèmes cruciaux de la TFD, avec des racines théoriques intéressantes. Nous avons vu que la TFD périodisait un signal non périodique, en créant une période fictive T égale à la fenêtre d'acquisition. Imaginons maintenant le signal soit périodique de période T0. Si T <> nT0, nous allons observer la création d'harmoniques parasites du plus mauvais effet ! Pour l'éviter, surtout si l'on ne connait pas T0, on choisit n assez grand, mais pas trop...
Pour l'instant, nous allons retenir que si la fenêtre d'acquisition est trop étroite, elle ne contient qu'une ou quelques périodes du signal, alors le résultat de la FFT sera mauvais parce que le nombre de points de calcul sera insuffisant. Par contre, si cette fenêtre est trop large, contenant beaucoup de périodes du signal, alors le résultat sera faux parce que le signal temporel sera imprécis (pas assez de points de calcul pour trop de temps). Il existe une limite à la précision du calcul. Cela ne vous rappelle-t-il rien?

Autre chose : il ne faut pas confondre le nombre d'échantillons dans le signal physique, qui dépend de la durée du signal et de sa fréquence d'échantillonnage, et le nombre de points de calcul de l'algorithme de FFT. Ce dernier est toujours une puissance de 2 (512, 1024, 2048 points ou plus). C'est une contrainte imposée par la nature de l'algorithme. Si vous ne la respectez pas, le calcul sera très nettement plus long, ou bien le logiciel que vous utilisez vous jettera.

Vous trouverez plus d'information sur la fonction fft de Python sur la page consacrée du site scipy.

Expériences avec Python

Expérience 1 - Un signal composite simple

Le but de l'expérience

Dans cette première expérience, nous établirons le spectre d'un signal composé par la somme de deux sinusoïdes d'amplitude et de fréquence différentes. Puisque nous le construisons, nous connaitrons le signal. On pourra ainsi vérifier que le spectre obtenu correspond bien au signal analysé.
Le signal choisi est:

A0*sin(f0*K*t) + A1*sin(f1*K*t)

Les coefficients A0 et A1 sont les amplitudes respectives des deux sinusoïdes. Les coefficients f0 et f1 en sont les fréquences. Le facteur est égale à \( 2\pi \). Leurs valeurs sont déclarées en tête de programme.

K = 2*pi # facteur conversion période/fréquence

A0 = 4 # amplitude fréquence fondamentale

A1 = 8 # amplitude première harmonique

f0 = 2 # fréquence fondamentale (Hz)

f1 = 8 # fréquence première harmonique (Hz)

Lançons le script FFTExperience1-0 dans l'IDE Python et observons les résultats.

Les résultats

Le programme affiche deux courbes. La première retrace l'évolution du signal dans le temps. Comme vous le constatez, il n'est pas évident de séparer à l'oeil la composante de 2 Hz et celle de 8 Hz !

La fenêtre suivante affiche le spectre d'amplitude du signal, c'est à dire la répartition d'énergie en fonction de la fréquence. Et là, nous voyons parfaitement les deux composantes fréquentielles, l'un à 2 Hz et l'autre à 8 Hz. Comme la valeur des amplitudes le laisse supposer, la composante à 8 Hz renferme le plus d'énergie, son pic est plus important.
Comme le signal est pur, il n'existe aucune énergie au delà des deux pics à 2 et 8 Hz. Le spectre est plat en dehors de ces deux pics.

FFTExperience1
Le programme

Le code source du script python FFTexperience1-0 est disponible ici. Pensez à le décompresser avant de l'utiliser.

Voyons les principales caractéristiques de ce script qui n'a vraiment rien de compliqué.

Il débute comme tous les scripts python par l'importation des modules et fonctions utilisés. Notons que les fonctions de fft existent dans numpy et scipy, avec des petites différences. J'ai utilisé la version de scipy.
Puis viennent les déclarations des constantes du signal et des paramètres d'échantillonnage.
Le vecteur temps est défni par:

t = linspace(t0, t1, N)

Cette instruction créée un vecteur d'une dimension nommé t, qui est composé de N points répartis entre t0 et t1.

Les instructions suivantes lancent le calcul de la FFT :

# définition des données de FFT

FenAcq = signal.size # taille de la fenetre temporelle

# calcul de la TFD par l'algo de FFT

signal_FFT = abs(fft(signal)) # on ne récupère que les composantes réelles

La fonction fft() renvoit ses résultats sous forme d'un vecteur de complexes. Seules les composantes réelles de l'amplitude m'intéressent pour tracer le spectre d'amplitude, c'est l'objet de la fonction abs() que d'en extraire le module.
Pour tracer le spectre, j'ai besoin de connaitre les fréquences réparties sur l'axe des abscisses du spectre, celles pour lesquelles les coefficients de Fourier ont été calculés. C'est l'objet de la fonction fftfreq dans l'instruction :

signal_freq = fftfreq(FenAcq,PerEch)

J'ai alors deux vecteurs signal_freq et signal_FFT, qui contiennent respectivement les fréquences des différentes harmoniques du signal et l'amplitude de chaque harmonique. Comme vous l'avez peut-être lu dans la doc, le domaine fréquentiel (la valeur des fréquences récupérées par fftfreq) est symétrique autour de 0 et réparti en fréquences positives et fréquences négatives. Les fréquences négatives ne m'intéressent pas, je vais donc les éliminer en retaillant le vecteur signal_freq en ne gardant que sa première moitié. Idem pour le vecteur signal_FFT sinon python ne serait pas content (les vecteurs doivent être de même taille pour tracer les courbes). C'est l'objet des instructions :

signal_FFT = signal_FFT[0:len(signal_FFT)//2]

signal_freq = signal_freq[0:len(signal_freq)//2]

Viennent après les instructions de tracé des courbes dans deux sous-fenêtres, qui ne présentent aucune difficulté.

Expérience 2 - Un signal composite bruité

Le but de l'expérience

Pour cette expérience, nous utiliserons le signal de l'expérience 1, auquel nous ajouterons une composante aléatoire de bruit. Nous ferons deux expériences: l'une en bruitant l'amplitude d'un signal de fréquence plus haute que le signal composite (50 Hz), l'autre en bruitant la fréquence du signal de bruit.

Dans les deux cas, les paramètres du signal de bruitage sont :

Ab = 20 # amplitude du bruit

fb = 50 # fréquence du bruit

Le signal bruité en amplitude sera défini par :

signal = signal + Ab*random.random(len(signal))*sin(fb*K*t)

Le signal bruité en fréquence sera défini par :

signal = signal + Ab*sin(random.random(len(signal))*fb*K*t)

L'instruction random utilisée ici produit un vecteur de longueur égale à celle du vecteur signal, rempli de nombres réels aléatoires dont la valeur est comprise entre 0 et 1.

Les résultats
Bruitage sur l'amplitude

Le calcul du spectre d'amplitude du signal bruité en amplitude donne le résultat suivant:

FFTExperience2-0

Sur le tracé du signal, il est pratiquement impossible de distinguer ses différentes composantes. L'amplitude varie aléatoirement. Mais le spectre fait très clairement ressortir les trois composantes à 2, 8 et 50 Hz, avec leur amplitude correspondante. La TFD met donc en évidence le bruit noyé dans le signal.
Sur le tracé du spectre, vous noterez de très petites fluctuations pour les fréquences autres que les trois fréquences composites du signal. Leur amplitude semble aléatoire. Il s'agit d'artefacts dus à l'algorithme de FFT, qui produit lui aussi des parasites !
Dans la pratique, ce type de bruitage correspond à une réalité bien connue : le parasitage de nos signaux électriques faibles par le champ électromagnétique créé par le courant EDF, de fréquence 50 Hz et d'amplitude faible mais persistante. Il ne reste plus qu'à trouver le moyen de l'éliminer, ce qui est le rôle des filtres.

Bruitage sur la fréquence

Le calcul du spectre d'amplitude du signal bruité en fréquence donne le résultat suivant :

FFTExperience2-1

Là aussi, il est pratiquement impossible de distinguer les différentes composantes du signal, mais on note toutefois une différence d'aspect avec le tracé précédent. L'amplitude du signal apparait moins aléatoire. Compréhensible, car l'amplitude du bruit est constante, seule sa fréquence change.
Le tracé du spectre ne comporte plus de pic à 50 Hz, contrairement au précedent. Mais on note la présence de fluctuations significatives à toutes les fréquences, dues aux variations aléatoires de la fréquence du signal de bruitage. L'amplitude de ces fluctuations est aléatoire.
En pratique, ce type de bruit est plus difficile à traiter que le bruit précédent. Dans le cas précédant, il suffisait de mettre en place un filtre qui éliminait le 50 Hz, sans connaitre nécessairement les fréquences utiles (2 et 8 Hz). Dans notre cas, il faut éliminer toutes les fréquences sauf les fréquences utiles, que l'on est donc supposer connaitre. Vous connaissez une application universelle de ce principe : votre transistor ! Lorsque vous voulez écouter France Info, votre transistor sélectionne parmi tous les signaux du champ électromagnétique ambiant, celui dont la fréquence centrale est 105.5 MHz (à Paris).

Le programme

Les code sources des scripts python FFTExperience2-0 et FFTExperience2-1 sont disponibles ici et .
Ils sont identiques au code du script FFTExperience1-0, à l'exception de la définition du bruit ajouté au signal. Je ne reviens donc pas sur le code.
ATTENTION : les scripts python sont disponibles sous format RAR pour des raisons techniques. pensez à les décompresser avant de les utiliser.

Pour conclure

Voilà une approche très succincte de la TFD et de son calcul avec Python. Nous aurons l'occasion de revenir sur le sujet dans la page de TangenteX. J'ai envie de travailler avec des signaux réels et de me pencher sur les moyens d'acquisition, faire un peu d'électronique quoi ! Je me suis procuré un microcontroleur Arduino, pas cher et facile à programmer. On doit pouvoir faire des choses! Je vous tiens au courant...


Contenu et design par Dominique Lefebvre - www.tangenteX.com mai 2016 -- Vous pouvez me joindre par mail ou sur PhysiqueX

Licence Creative Commons

Cette œuvre est mise à disposition selon les termes de la Licence Creative Commons Attribution - Pas d’Utilisation Commerciale - Pas de Modification 3.0 France .