(Le bruit aléatoire pure est bleu. Avec la superposition du signal : rouge. Vous voyez?...)
Pour le traitement de signaux (scientifiques, musique, parole, données sismographiques, etc.) et aussi des images (filtrage, compression, morphologie), la transformée de Fourier est l'outil numéro UN, et les techniques de calcul de la FT font partie de la panoplie de chaque scientifique et chaque ingénieur.
Je pense qu'une Ecole d'ingénieurs doit former les concepteurs et en général : des créateurs, y compris les créateurs des programmes, et algorithmes, voici donc une introduction à la technologie FFT (Fast-Fourier Transform), qui réduit la complexité de la transformée sur \(N\) points, de \(N^2\) à \(N \log N\). Ceci n'est pas une digression sans relation avec mon cours ! Le fix d'un satellite exploite la FFT tout le temps pour la synchro (le max. d'autocorrélation).
On aura besoin d'un peu d'algèbre, et un peu d'expérience avec la programmation, c'est tout. Dans la littérature il existe plusieurs conventions (surtout de normalisation), il faut être attentif. On n'abordera que la transformée discrète. Voici la formule en question.
La séquence \(x_k\), en général complexe peut avoir la longueur quelconque, mais pour faciliter la présentation, on choisira \(N = 2^K\), une puissance de 2. Actuellement on a la FFT généralisée pour un nombre quelconque de valeurs, mais la présentation devient considérablement plus compliquée. La FFT binaire est standard et facile.
En regardant la formule on constate que pour chaque élément de la suite - résultat, il faut parcourir le tableau source. Donc, on effectuera \(N\) fois quelques opérations algébriques, dont le coût temporel est considéré constant. Pour la transformée complète, \(N\) valeurs de sortie, on aura un nombre d'opérations de base proportionnel à \(N^2\).
Faisons quelques exemples concrets, pour voir la relation entre la source traitée comme une fonction discrète, et le résultat. Nous allons programmer en Python/numpy, en évitant les boucles Python.
Ceci est déjà implémenté dans quelques paquetages, notamment en numpy, car FFT est une technique cruciale. La fonction s'appelle fft.
Le rapport entre la formule, et la transformée continue, doit être mentionnée, car à cause des conventions, la forme du résultat n'est pas la même.
Mathématiquement on définit la transformée de Fourier par
$$
\hat{x}(\xi) = \int_{-\infty}^{\infty} x(t) e^{2\pi i t \xi} dt\,.
$$
Ceci est symétrique autour de zéro, tandis que les indices dans la forme discrète sont positifs. On interprète alors les indices grands, qui dépassent \(N/2\) comme négatifs. Observez que
\(
\exp(2\pi i k (N-m)/N) = \exp(-2\pi i k m/N)\,.
\)
Parfois on utilise pour la présentation une manipulation conventionnelle, qui échange les moitiés gauche et droite du tableau transformé ; les indices grands deviennent "négatifs".
Avec les fréquences entières, et les amplitudes \([(\zeta,a)]\) : [(0,1),(1,-0.8),(3,0.5),(7,0.2),(12,-0.3),(50,0.1)], la somme des signaux périodiques \(a \cdot \exp(2\pi i k \zeta /N)\), traçons le signal (real() et imag()) et sa transformée :
Voici l'algorithme. Séparons le tableau \(x\) en deux sous-séquences, avec les indices entiers pairs, et impairs. Définissons \(M=N/2\).
$$
f_m = \sum_{k=0}^{N-1} {x_k e^{-\frac{2\pi i}{N}m k}} =
\sum_{k=0}^{M-1} {x_{2k} e^{-\frac{2\pi i}{N}m (2k)}} + \sum_{k=0}^{M-1} {x_{2k+1} e^{-\frac{2\pi i}{N}m (2k+1)}} = \\
\sum_{k=0}^{M-1} {x_{2k} e^{-\frac{2\pi i}{M}m k}} +
e^{-\frac{2\pi i}{N}m} \cdot \sum_{k=0}^{M-1} {x_{2k+1} e^{-\frac{2\pi i}{M}m k}}\,.
$$
Les deux sommes sont des transformées de Fourier sur les moitiés de points ! Enfin, presque. L'indice \(m\) varie toujours entre 0 et N-1, non pas M-1, donc les espaces de l'original et de la transformée ne sont pas conformes. Mais l'exponentielle complexe présente quelques symétries utiles, qui simplifieront le calcul. Prenons la première moitié de la transformée, avec \(m = 0, 1, \ldots, M-1\), et l'autre moitié portera des indices \(m+M\).
Mais \( e^{-\frac{2\pi i}{M}(m+M) k} = e^{-\frac{2\pi i}{M}m k} \). Par contre, le facteur devant la somme :
$$
e^{-\frac{2\pi i}{N}(m+M)} = e^{-\frac{2\pi i}{N}m} \cdot e^{-\frac{2\pi i}{2M}M} = -e^{-\frac{2\pi i}{N}m}
$$
puisque \( e^{-\pi i} = -1 \), la racine primitive de 1. Donc
$$
f_{m+M} = \sum_{k=0}^{M-1} {x_{2k} e^{-\frac{2\pi i}{M}m k}} -
e^{-\frac{2\pi i}{N}m} \cdot \sum_{k=0}^{M-1} {x_{2k+1} e^{-\frac{2\pi i}{M}m k}}\,.
$$
... et la seconde moitié de la transformée peut réutiliser les deux sommes de la première moitié. Il reste d'effectuer la soustraction, ce qui ajoute un nombre d'opération proportionnel à \(M\).
Avant de passer au codage, ce qui est un exercice en programmation vectorielle (numpy), intéressant en soi, analysons la complexité calculatoire de la solution. Que le nombre d'opérations sur des tableaux de taille \(N\) soit noté comme \(C(N)\). La réduction récursive donne la relation suivante, avec une constante inconnue \(w\) : $$ C(2N) = 2 \cdot C(N) + 2w N\,. $$ On cherche le comportement de cette fonction pour \(N\) assez grands. Définissons \(D(N) = C(N)/N\). Alors $$ D(2N) = D(N) + w\,. $$ On peut facilement conclure, que le comportement de \(D\) est logarithmique, c'est la plus simple fonction avec la propriété que si on multiplie l'argument, on ajoute au résultat : \(\log(ax) = \log(a)+\log(x)\).
Donc, finalement, \(C(N) = c \cdot N \log(N)\). Les constantes multiplicatives, la base du logarithme, et autres détails accessoires restent inconnus, mais asymptotiquement le résultat est évident. Nous avons réduit la complexité de la transformation de manière très forte.
Tout d'abord, si \(N=2\), le facteur exponentiel se réduit à -1, et le résultat est la somme et la différence de deux éléments du tableau. On peut même descendre à \(N=1\), alors la transformée est triviale. La récursivité fait le reste. Le programme est vraiment très court, et il n'est pas optimisé.
def mfft(x):
N=len(x)
if N==1: return array([x[0]])
M=N//2; fc=exp(-1j*pi/M*arange(M))
fev=mfft(x[0::2]); fod=fc*mfft(x[1::2])
return concatenate((fev+fod,fev-fod))
On pourrait l'accélérer quelques dizaines de fois, sans changer son architecture. L'implémentation dans numpy est entre 50 et 200 fois plus rapide...
On peut aussi renoncer de la récursivité, en faveur d'itérations.
Observez la structure dynamique du programme : de manière indépendante à gauche et à droite, la première manipulation est la construction des deux segments, pair et impair, et ceci "descend" jusqu'à ce que la longueur des segments ne devient égal à 1.
Ensuite on combine les segments : sommes / différences, et on remonte la "pyramide récursive", en insérant les sommes et les différences dans des segments du tableau pré-alloués.
On peut d'abord déplacer les éléments pairs à gauche, impairs à droite, et répéter cela jusqu'à la longueur 2, et ensuite gérer les segments (slices) contigus du tableau, sans construire une multitude de tableaux secondaires.
Les optimisations très poussées dépassent les bornes de notre introduction. Cooley-Tukey n'est pas la seule implémentation possible. On peut procéder à l'envers, séparer les moitiés gauche et droite, et reconstruire les sgments pairs et impairs. L'usage du parallélisme est de rigueur. Pour N petit, on peut coder tout "in-line", sans procédure et sans boucles. Les FFT les filtrages, etc., c'est toute une technologie qui doit être connue, et pas seulement dans le cadre du traitement des signaux.
Quand même,deux mots sur le ré-arrangement des indices...
Nos appels récursifs contenaient : ... mfft(x[0::2]) ... mfft(x[1::2]) et cette manipulation qui déplace chaque élément plusieurs fois dans un nouveau tableau (ou sur place), est visiblement inefficace. On peut faire ce réarrangement en cascade une seule fois. On peut prouver mathématiquement que si l'indice de l'élément déplace, x[k] est k, le nouveau indice k' aura la représentation binaire obtenue de k par le renversement de ses bits. Par exemple, si N=16, k=11, signifie k='0b1011'. Alors k'=0b1101 = 13.
Voici la procédure en question, très courte. Il y a une boucle en Python, mais la longueur est le logarithme de N, pour un tableau de taille 1 million, c'est 20.
def rearr(ar): # res pré-alloué N=len(ar); c=zeros(N,dtype=int); L=1 while N>1: N//=2; c[L:2*L]=N+c[:L] L*=2 return ar[c] # [0,8,4,12,2,10,6,14,1,9,5,13,3,11,7,15]
Rappelons encore : l'algorithme FFT est hautement parallélisable, et il y a des systèmes multi-processeurs permettant d'effectuer les FFT de taille très considérable en temps réel.
Ceci a ses applications dans la gestion des données satellitaires.
Un classique exemple est le problème de synchronisation. Nous savons déjà que si on accumule plusieurs échantillons d'un signal bruité, le rapport données/bruit augmente. Mais il faut savoir où se trouve le vrai signal dans la séquence.
(Le bruit aléatoire pure est bleu. Avec la superposition du signal : rouge. Vous voyez?...)
On peut calculer la fonction d'auto-corrélation de ce signal. L'auto-corrélation de la séquence \(y\) est définie comme suit : $$ \Large R(k) = \sum_n y_n \cdot \bar{y}_{n-k}\,. $$ Avec la FFT, on calcule cette fonction par \(t=F(y)\,; S = t\cdot \bar{t}\,; R = F_{\text{inv}}(S)\,.\)
|
|