PrograZine #7 - La transformée de Fourier rapide

[ edito | sommaire | contribution | contacts ]

[>La transformée de Fourier rapide<]

Don Cross
dcross@intersrv.com
traduit par Jean Debord (JDebord@compuserve.com)

>Introduction à l'audio numérique

>Information fréquentielle d'une fonction du temps

cochlée

>La transformée de Fourier en tant que concept mathématique

₀

>L'algorithme de la transformée de Fourier rapide

Vos tableaux d'entrée et de sortie doivent avoir la même taille, que nous appellerons n.

La valeur de n doit être une puissance entière positive de 2. Par exemple, un tableau de taille 1024 est permis, mais pas un de taille 1000. La plus petite taille permise est 2. Il n'y a pas de limite supérieure à la taille de n à part les limitations liées à l'allocation de mémoire pour les tableaux (entrée{réel, imaginaire}, sortie{réel, imaginaire}) et au temps d'exécution de cet algorithme en O(n*log(n)).

Le tableau realOut contient les coefficients des cosinus dans la formule de Fourier.

Le tableau imagOut contient les coefficients des sinus dans la formule de Fourier.

L'ordre des fréquences dans les tableaux de sortie (realOut et imagOut) est un peu étrange, parce que ces fréquences peuvent être positives aussi bien que négatives. Les deux sont nécessaires pour des raisons mathématiques quand les entrées ont des valeurs complexes (c'est-à-dire quand l'une au moins des entrées a une composante imaginaire non nulle). La plupart du temps, la FFT s'utilise sur des données réelles, notamment dans le cas des signaux audio. La FFT, lorsqu'elle reçoit une entrée à valeurs réelles, fournit une sortie dont les fréquences positives et négatives sont redondantes, en ce sens que les valeurs correspondantes sont conjuguées complexes l'une de l'autre, ce qui signifie que leurs parties réelles sont égales et que leurs parties imaginaires sont opposées. Si vos entrées ont toutes des valeurs réelles, toute l'information dont vous avez besoin se trouve dans la première moitié des tableaux de sortie.

Le premier élément en sortie, realOut[0] et imagOut[0], contient la valeur moyenne de toutes les entrées. Pour les indices i = 1, 2, 3, ..., n/2, la valeur de la fréquence exprimée en Hz est f = samplingRate * i / n. Les contreparties négatives des fréquences positives en i = 1, 2, 3, ..., n/2 - 1 correspondent aux indices i' = n - i. Voici un exemple : supposons que la fréquence d'échantillonnage soit de 44100 Hz, et que nous utilisions des tampons de 1024 nombres complexes pour l'entrée et la sortie. La fréquence en i = 1 serait (44100 Hz) * 1 / 1024 = 43.07 Hz. La fréquence négative correspondante serait en i = 1024 - 1 = 1023 (soit le dernier élément du tableau), avec une valeur de -43.07 Hz. De même, i = 17 correspondrait à une fréquence de (43.07 Hz)*17 = 732.13 Hz, tandis que i = 1024 - 17 = 1007 correspondrait à -732.13 Hz, etc.

L'indice n/2 est un cas particulier : il correspond à la fréquence de Nyquist, qui est toujours égale à la moitié de la fréquence d'échantillonnage dans un enregistrement PCM. Par exemple, la fréquence de Nyquist d'un CD audio vaut (44100 Hz)/2 = 22050 Hz. La fréquence de Nyquist est importante car c'est la plus haute fréquence qu'un enregistrement PCM peut reproduire. Le PCM ne peut rien représenter au-dessus de cette fréquence. Si vous essayez d'échantillonner un signal de fréquence supérieure à celle de Nyquist, il s'ensuit une une sévère distorsion, connue sous le nom d'aliasing, et analogue aux illusions d'optique des vieux westerns qui font que les roues semblent tourner à l'envers. Les ingénieurs du son savent qu'ils doivent filtrer tous les signaux analogiques au-dessus de la fréquence de Nyquist avant de les numériser, afin d'éviter les problèmes d'aliasing. Notez que cette limitation est inhérente à l'enregistrement numérique lui-même, qu'il fasse ou non l'objet de la FFT.

Si le signal d'entrée de la FFT est réel, le résultat correspondant à la fréquence de Nyquist, d'indice n/2, aura toujours une valeur réelle (ce qui signifie que sa partie imaginaire sera nulle, ou proche de zéro du fait des erreurs d'arrondi). En effet, d'après la formule précédente, i' = n - n/2 = n/2, donc la fréquence de Nyquist est égale à sa contrepartie négative. Le résultat correspondant doit donc être égal à son conjugué complexe, il ne peut donc s'agir que d'un nombre réel.

Parfois vous ne vous intéresserez qu'au module (magnitude) ou à l'argument (angle) du résultat, pour chaque composant de fréquence. Voici comment les calculer :
#include double magnitude = sqrt ( realOut[i]*realOut[i] + imagOut[i]*imagOut[i] ); double angle = atan2 ( imagOut[i], realOut[i] );
l
Pour faire la conversion inverse, de la magnitude et l'angle vers les parties réelles et imaginaires, utilisez le code suivant :
#include double real = magnitude * cos(angle); double imag = magnitude * sin(angle);
l

Pour une meilleure compréhension mathématique de la transformée de Fourier discrète, l'équation suivante donne la relation exacte entre l'entrée et la sortie. Dans cette équation, x_k est le k^ième nombre complexe en entrée (domaine temporel), y_p est le p^ième nombre complexe en sortie (domaine fréquentiel), et n = 2^N est le nombre total d'échantillons. Notez que k et p sont dans l'intervalle 0 .. n-1.

Cette formule définit la transformée de Fourier discrète, mais pas l'algorithme FFT. En effet, cette formule requiert un nombre d'opérations de l'ordre de n², alors que la FFT est de l'ordre de n*log₂(n). En d'autre termes, si vous utilisiez la formule ci-dessus, le calcul serait beaucoup plus long qu'avec l'algorithme FFT. Cependant, si vous n'avez besoin que d'un petit sous-ensemble du spectre de fréquences (disons deux ou trois points), ou si le nombre de points n'est pas une puissance de 2, cette formule combinée à quelques optimisations trigonométriques peut être utile. Notez que c'est exactement ce que j'ai fait pour mon code source en Turbo Pascal dans la procédure appelée CalcFrequency.

>Applications de la FFT

exemple de programme en C++

filtres temporels

>Annexes

>Code source en Pascal

fourier.pas

testfft.pas

testfft.zip

>Code source en C/C++

fft.zip

Numerical Recipes in C

fftdom.cpp

₀

>Autres sites sur la FFT

FFTW - The "Fastest Fourier Transform in the West" - développé au MIT par Matteo Frigo et Steven G. Johnson. Regardez ici si vous êtes intéressé par un code source FFT portable pour l'analyse multi-dimensionnelle, ou si vous avez un tampon dont la taille n'est pas une puissance de 2.

The Fourier Page par Douglas Nunn. Beaucoup de liens intéressants... certains techniques, d'autres historiques.

Public Domain FFT Code - Codes source du domaine public en Fortran et C. Avec une belle discussion théorique de la transformée de Fourier continue, dont dérivent toutes les transformées discrètes (y compris la FFT).

FFT for 2D matrices in Fortran par Hon W Yau, Bhaven Avalani et Alok Choudhary.

Numerical Recipes in C - La version online de cet ouvrage incontournable. Le chapitre 12 présente la théorie de la FFT, et le chapitre 13 discute les applications. La consultation requiert un visualiseur Adobe Acrobat ou PostScript.

La page de Numerical Recipes.

>Sites apparentés

Pitch Shifting FAQ - par Prosoniq Products

Amara's Wavelet Page - Nombreux liens sur les ondelettes.

A Basis for Functions: The Fourier Series - Présente quelques concepts fondamentaux..

techniques de filtrage dans le domaine temporel

Numeric section de la Turbo Pascal Programmers Page [utilise des cadres].

FXT Page de Jörg's "Useful and Ugly Pages"

Ce document est la propriété de Don Cross. La rédaction ainsi que l'auteur déclinent toutes résponsabilité vis-a-vis des conséquences que pourrait entrainer le contenu de cet article.