Cet article fait suite au travail sur la transformée de Fourier rapide que j’ai dû réaliser pour la fac avec mon collègue Ivan. Je précise tout de suite que cette implémentation s’applique exclusivement aux images carrées, en niveaux de gris et dont la dimension est de la forme 2ⁿ.

L’utilisation de java n’est peut être pas idéale dans la mesure où il s’agit d’un algorithme ayant une complexité assez élevée, d’autant plus que l’API java ne propose pas de classe permettant de coder les nombres complexes (en tout cas pas à ma connaissance), mais dans la mesure où il nous fallait une interface graphique et où je n’ai pas encore trouvé le temps d’adopter une bibliothèque graphique en C++… Le côté positif est que la lecture du code en sera facilitée.

Le tout… n’est plus disponible du fait qu’un trop grand nombre d’étudiants médiocres semble s’être contenté de copier le code sans même dire merci.Hop.

Rappel.
Rappelons tout d’abord la formule de base de la transformée de Fourier discrète:

Concrètement, cette reformulation revient à séparer les pixels qui composent l’image en deux groupes, d’un côté les pixels dont la position est paire (à gauche sur ce schémas), de l’autre ceux dont la position est impaire (à droite).

Or les deux sommes qui composent notre nouvelle formule sont exprimées sous la même forme que la formule d’origine. On perçoit donc le premier intérêt de cette transformation: elle pourra être exprimée au moyen d’une fonction récursive. Le deuxième intérêt est que la complexité de cet algorithme est inférieure à l’originale: sur une image à deux dimensions on obtient une complexité de N²ln²(N) (contre N⁴ pour l’algorithme d’origine).

Pour appliquer cet algorithme sur une image à deux dimensions, il faut d’abord l’appliquer sur toutes les lignes de l’image, sauvegarder le résultat sous forme d’une image intermédiaire, puis l’appliquer de nouveau sur toutes les colonnes de ce résultat. Mais suivre à la lettre cette démarche nous obligerait à écrire deux fois l’algorithme de la transformation, une fois pour les lignes et une autre pour les colonnes, à cause de la gestion des indices. Pour contourner cette difficulté, nous travaillerons exclusivement sur des images carrées et nous effectuerons sur l’image intermédiaire une transformation qui inverse lignes et colonnes, procédé que nous appliquerons de nouveau sur le résultat final.

Implémentation de l’algorithme.
La première méthode qui vient à l’esprit pour implémenter cet algorithme consiste à regrouper dans des structures de données les pixels pairs d’un côté, les pixels impairs de l’autre, et à subdiviser ainsi chaque structures jusqu’à ce qu’elles ne contiennent plus qu’un pixel chacune.
Mais cette méthode est inutilement lourde; il est en effet possible de retrouver les pixels qui nous intéressent en sauvegardant seulement les quelques informations que sont la position du pixel de départ, le pas qui sépare les pixels les uns des autres et le nombre de pixels qui compose la subdivision.

Si l’on s’intéresse à la position des pixels dans le tableau d’origine, on constate que le pas d’un niveau donné est égal au pas du niveau supérieur que multiplie 2. On en déduit donc que, pour un niveau n donné (0 étant le sommet), les pixels seront séparés par un pas de 2ⁿ.

Maintenant que l’on connaît le pas qui sépare deux pixels, il nous faut déterminer la position du premier pixel de chaque subdivision. En fait, pour une subdivision paire, cette valeur sera égale à la position du premier pixel du niveau supérieur et, pour une subdivision impaire, à la position du premier pixel augmenté du pas du niveau supérieur.

Le nombre de pixels qui compose une subdivision est quant à lui égal au nombre de pixels qui compose la subdivision du niveau supérieur, divisé par deux.

On peut désormais écrire l’algorithme récursif de la transformée de Fourier rapide:

private Complexe fourierRapide(int size, int k0, int pas) {
	Complexe exposant = new Complexe(0, -2 * Math.PI * pixel / (size + 1));
	if (size > 0)
		return fourierRapide((size - 1) / 2, k0, pas * 2).add(exposant.exp().mult(fourierRapide((size - 1) / 2, k0 + pas, pas * 2)));
	else
		return exposant.mult(k0).exp().mult(image1D.getPixelComplexe(0, k0));
}

L’algorithme de la transformée inverse est en tout point identique, à l’exception du signe de la variable exposant qui devient positif.

Pour accélérer encore un peu plus la transformation, nous avons implémenté une version multithread de l’algorithme, chaque thread se partageant les lignes qui composent l’image. Un système de barrière permet de synchroniser les principales étapes de la transformation (calcul des lignes / passage des lignes aux colonnes / calcul des lignes / passage des lignes aux colonnes).