1. Représentation de la décomposition en série de Fourier
La transformée de Fourier est basée sur la découverte que toute fonction périodique du temps x(t) peut être décomposée en une somme infinie de sinus et cosinus dont les fréquences commencent à zéro et augmentent par multiples entiers d’une fréquence de base f0 = 1/T, où T est la période de x(t). Ce développement se présente ainsi :
L’expression du membre de droite est appelée série de Fourier. Le rôle de la transformée de Fourier est de trouver toutes les valeurs ak et bk qui composent la série, connaissant la fréquence de base et la fonction x(t). Nous pouvons considérer le terme a0 à l’extérieur de la somme comme le coefficient du cosinus pour k=0. Il n’y a pas de coefficient correspondant b0 pour le sinus car le sinus de zéro vaut zéro, et donc ce coefficient n’aurait pas d’effet.
Bien sûr, aucun ordinateur réel ne peut calculer de sommes infinies ; il nous faut donc trouver un ensemble fini de sinus et cosinus. C’est facile à faire pour une entrée numérique échantillonnée, si nous stipulons qu’il y a autant de fréquences en sortie qu’il y a de valeurs temporelles en entrée. Nous profitons aussi du fait que tout enregistrement audio numérique a une longueur finie. Nous pouvons prétendre que la fonction x(t) est périodique, et que la période est égale à la longueur de l’enregistrement. En d’autres termes, imaginons que l’enregistrement se répète continuellement, et appelons x(t) cette fonction. La durée de la section répétée définit la fréquence de base f0 dans l’équation ci-dessus. En d’autres termes, f0 = sampling Rate / N, où N est le nombre d’échantillons dans l’enregistrement.
Par exemple, si nous utilisons un taux d’échantillonnage (sampling Rate)de 44100 échantillons / seconde, et que la longueur de l’enregistrement (N) est de 1024 échantillons, la durée représentée par l’enregistrement est 1024 / 44100 = 0.02322 seconde, de sorte que la fréquence de base f0 sera 1 / 0.02322 = 43.07 Hz. Si nous soumettons ces 1024 échantillons à la FFT, nous obtiendrons les coefficients ak et bk des sinus et cosinus pour les fréquences 43.07Hz, 2*43.07Hz, 3*43.07Hz, etc
L’applet suivante montre la décomposition en série de Fourier de signaux de forme simple comme le triangle, le carré et l’impulsion. Il est intéressant de visualiser l’effet du nombre de termes de la série sur la forme du signal.
2. Calcul de la transformée de Fourier
La transformée de Fourier discrète est un algorithme qui convertit une fonction du temps à valeurs complexes échantillonnées en une fonction à valeurs complexes de la fréquence, également échantillonnée. La plupart du temps, nous travaillons sur des fonctions à valeurs réelles, donc nous fixons toutes les parties imaginaires des entrées à zéro.
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, xk est le kième nombre complexe en entrée (domaine temporel), yp est le pième nombre complexe en sortie (domaine fréquentiel), et n = 2N est le nombre total d’échantillons. Notez que k et p sont dans l’intervalle 0 .. n-1.
Souvent on ne s’intéresse qu’au module (amplitude) ou à l’argument (phase) du résultat, pour chaque composant de fréquence.
amplitude = sqrt ( PartieReel*PartieReel + Partie Imaginaire*Partie Imaginaire );
Phase = arctan ( Partie Imaginaire, Partie Réel );
Pour faire la conversion inverse, de la magnitude et l’angle vers les parties réelles et imaginaires :
Partial Reel = amplitude * cos(phase);
Partie Imaginaire = amplitude * sin(phase);
3. FFT (Fast Fourier Transform)
La transformée de Fourier rapide est simplement un algorithme permettant de réduire le nombre d’opérations pour calculer la DFT (Discrete Fourier Transform).
| Rappelons la relation permettant de calculer la DFT : |  |
(1) |
| Par conséquent, les opérations à effectuer pour obtenir les N valeurs de la DFT sont : | ||
| (N.N) multiplications complexes | ||
| N(N-1) additions complexes | ||
| Bien entendu, les multiplications complexes ont une durée d’exécution beaucoup plus longue que les additions. | ||
4. Algorithme de Cooley-Tukey
L’algorithme de FFT le plus connu est celui de Cooley-Tukey, également appelé algorithme de réduction à base 2 dans le domaine temporel. Cet algorithme requiert un nombre d’échantillons qui soit une puissance de 2. Par exemple : N=128, 256, 4096, etc.
D’autres algorithmes FFT existent qui requièrent d’autres exigences, mais nous ne les considérerons pas ici.
Dans le cas d’une FFT selon l’algorithme de Cooley-Tukey, le nombre d’opérations est considérablement réduit :
 |
multiplications complexes | ||
Exemple :
Prenons 8 échantillons, qui ont les valeurs successives suivantes s0, s1, s2, s3, s4, s5, s6, s7.
La DFT se présente ainsi :
En séparant les échantillons pairs et impairs et en factorisant les nombres impairs, on peut mettre en évidence des similitudes entre les opérations pratiquées sur ces échantillons. Prenons comme exemple la ligne S1, en séparant les échantillons pairs/impairs, puis en factorisant les échantillons impairs.
Traduction en langage mathématique :
On voit que les termes W sont identiques pour les couples 0-1, 2-3, 4-5, 6-7 à un facteur exp(-jpi/4) près. On peut encore subdiviser les N/2 échantillons de nos parenthèses :
On voit qu’en factorisant, on arrive à former les couples s0-s4, s2-s6, s1-s5, s3-s7, puis à remonter dans les parenthèses pour arriver à S1.
Si on prend en compte les huit lignes, cela nous donne :
En langage plus mathématique :
En appliquant les règles (4) et (5), on simplifie aisément les facteurs WN nk. Le résultat est le suivant :
4.1 Opération “Butterfly”
L’opération “butterfly” ou papillon est la suivante :
 |
4.2 Application à la FFT
L’opération “Butterfly” est particulièrement adaptée à l’algorithme FFT, comme on peut le voir dans la Figure 2. Cette figure correspond directement à la résolution de l’énoncé (6).
Figure 2
4.3 L’incrément en “reverse carry”
Comme on peut le voir sur la Figure 2, les couples d’échantillons doivent être choisis au départ selon un ordre particulier : s0-s4, s2-s6, etc.
Cette incrémentation particulière est appelée “reverse carry” (retenue inverse). Elle consiste à additionner N/2 à l’indice, mais à reporter la retenue à droite plutôt qu’à gauche.
Exemple pour 8 échantillons (N/2 = 4) :
Les processeurs spécialisés type DSP proposent tous un mécanisme d’incrémentation de ce type pour leurs pointeurs.
4.4 Réduction dans le domaine fréquentiel
Les équations décrites dans (6) peuvent encore être arrangées autrement :
La structure de l’algorithme n’est pas différente, simplement ce sont les index qui changent :
les échantillons sont ordonnés linéairement, alors que les valeurs de fréquence sont ordonnés en “reverse carry” (situation inverse de la Figure 2).
Exemple d’algorithme FFT en “C”. Notons que l’on utilise deux Butterfly aux seins de la boucle principales.
void FFT( SHORT INT *fix_x ) {
/* Local declarations */
SHORT INT fix_xt, fix_t1, fix_t2;
SHORT INT fix_cos,fix_sin;
SHORT INT fix_temp,fix_temp1,fix_temp2;
SHORT INT *fix_src1,*fix_src2,*fix_src3,*fix_src4;
SHORT INT *fix_dst1;
…
/*** Compute digit reverse counter ***/
j = 0;
for ( i = 0; i < NB_SAMPLES – 1; i++ ) {
if ( i < j ) {
fix_xt = fix_x[j];
fix_x[j] = fix_x[i];
fix_x[i] = fix_xt;
}
k = NB_SAMPLES >> 1;
while ( k < (j+1) ) { j -= k; k >>= 1; }
j += k;
}
….
/*** butterflies ***/
n2 = 1;
n1 = 2;
n8 = NB_SAMPLES >> 1;
for ( k = 1; k < NB_STAGE; k++ ) {
n4 = n2;
n2 = n1;
n1 <<= 1;
n8 >>= 1;
for ( i = 0; i <= NB_SAMPLES-1; i += n1 ) {
fix_xt = fix_x[i] >> 1;
fix_temp = fix_x[i+n2] >> 1;
fix_x[i] = fix_xt + fix_temp;
fix_x[i+n2] = fix_xt – fix_temp;
fix_x[i+n4+n2] = – fix_x[i+n4+n2];
it = -1;
fix_src1 = &fix_x[i+1];
fix_src3 = &fix_x[i+n2+1];
fix_src2 = &fix_x[i+n2-1];
fix_src4 = &fix_x[i+n1-1];
for ( j = 1; j <= n4-1; j++ ) {
it += n8;
fix_cos = Table[NB_SAMPLES/4-2-it]; // Table COSINUS / SINUS
fix_sin = Table[it];
fix_temp1 = *fix_src3; // SOURCES
fix_temp2 = *fix_src4;
fix_t3 = ( fix_temp1 * fix_cos ); // (S1*COS) + (S2*SIN)
fix_t3 += ( fix_temp2 * fix_sin );
fix_t1 = fix_t3 >> 16;
fix_t4 = ( fix_temp1 * fix_sin ); //(S1*SIN)-(S2*COS)
fix_t4 -= ( fix_temp2 * fix_cos );
fix_t2 = fix_t4 >> 16;
fix_temp2 = *fix_src2 >> 1;
*fix_src4– = fix_temp2 – fix_t2; // Butterfly 1
*fix_src3++ = – fix_temp2 – fix_t2;
fix_temp2 = *fix_src1 >> 1; // Butterfly 2
*fix_src2– = fix_temp2 – fix_t1;
*fix_src1++ = fix_temp2 + fix_t1;
}
…
}
4.5 FFT sur base d’un nombre d’échantillons N quelconque
4.5.1 L’algorithme FFT pour N=k*j
A partir des règles (4) et (5), on peut créer des algorithmes de FTT pour N=k×j avec k et j entiers. L’efficacité est cependant moindre que lorsque N est une puissance de 2.
4.5.2 L’algorithme de base mixte
Selon §2.5.1, il est tout à fait possible de travailler avec une base (radix) différente de deux. Parfois, la base change au cours d’opération (split-radix) ! Par exemple, les premières “butterfly” travaillent en 2×2, puis les suivantes en 4×4.
4.5.3 L’algorithme de Rader
Lorsque N est premier, on peut imaginer que la DFT est irréductible. Ce n’est pas le cas, un algorithme particulier (l’algorithme de Rader) est au contraire très performant !
4.6 Nombre de points pour une FFT
La fréquence d’échantillonnage doit être au moins deux fois plus élevée que la bande passante du signal (théorème de nyquist). Le rapport de 2.56 est souvent choisi pour calculer des FFT. Imaginons que notre bande passante soit sur 100 points, le nombre d’échantillons sera 256, soit une puissance de 2 :
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 numérique. 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 numérique peut reproduire. L’enregistrement ne peut rien représenter au-dessus de cette fréquence. Si nous essayons d’interpréter 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. Notons que cette limitation est inhérente à l’enregistrement numérique lui-même, qu’il fasse ou non l’objet d’une FFT.
4.7 Inconvénient de la série de Fourier
Le principal inconvénient de la série de Fourier dans le traitement du son est que l’enregistrement numérique doit être divisé en blocs de n valeurs, or la série de Fourier associée se comporte comme étant une fonction périodique ce qui n’est pas le cas de nos blocs définis précédemment. La transformée de Fourier induit des phénomènes de bords qui sont indésirables et non représentatifs du signal d’entrée original. L’utilisation de fenêtres permet de réduire ces phénomènes de bords.