Jouons à implémenter une transformée de Fourier rapide !

Un algorithme que vous utilisez probablement au quotidien.

a marqué ce sujet comme résolu.

Tout le monde se secoue ! :D

J’ai commencé (dimanche 23 mai 2021 à 14h53) la rédaction d’un tutoriel au doux nom de « Jouons à implémenter une transformée de Fourier rapide ! » et j’ai pour objectif de proposer en validation un texte aux petits oignons. Je fais donc appel à votre bonté sans limites pour dénicher le moindre pépin, que ce soit à propos du fond ou de la forme. Vous pourrez consulter la bêta à votre guise à l’adresse suivante :

Merci !

Hello, pour un projet sur mon temps libre (que je présenterais peut-être plus tard) j’ai eu besoin d’implémenter une FFT efficace. Après quelques semaines à potasser Numerical Recipes, je me suis dit que ce serait sympa d’écrire un petit contenu qui résume ce que j’ai pu (ré)apprendre. J’espère que le format de mini-tuto est adapté, mais je peux changer au besoin.

La totalité du fond est présente, mais j’ai des doutes sur la forme : je ne suis pas très expérimenté en terme de vulgarisation, et en plus je me suis lancé le défi d’écrire tout ça ce week-end. Il y a donc probablement beaucoup à redire ! Mais j’ai très envie d’apprendre à vulgariser correctement, et j’attends donc avec impatience vos retours. :)

+1 -0

Je trouve ça assez chouette comme genre de tutoriel, qui rentre pas mal dans les détails, mais avec plusieurs niveaux de “difficulté” qui arrivent progressivement.

J’ai fait une première lecture rapide, trop rapide pour pouvoir formuler encore des retours précis sur les détails, mais j’ai quand même quelques points qui me sont venus à l’esprit.

Je pense qu’il manque peut-être une motivation un peu générale ; tu l’as dit au tout début, la vocation de ce tutoriel n’est pas de “présenter la transformée de Fourier”. Je comprends que ce n’est pas forcément le lieu pour faire une introduction de l’objet mathématique, et que ce n’est pas l’objet du tutoriel. Je pense néanmoins que quelques points “historiques’ qui ont motivé la recherche d’un algo efficace pour calculer des transformées de Fourier pourraient quand même valoir le coup d’être mentionnés et rester pertinent pour ton public. Une manière de voir ça est que tout au long du tutoriel, l’étalon que tu utilises pour voir où en est ton implémentation est celui de la comparaison avec une implémentation de référence. Je trouve que c’est intéressant, et démystifie un peu les boîtes noires que peuvent être ce genre d’implémentation. Mais je pense aussi que ça pourrait gagner en pédagogie en expliquant typiquement pourquoi la première implémentation naïve n’est pas suffisante : en voyant les chiffres on pourrait se dire que quelques dizaines de millisecondes ça reste quand même assez rapide. Peut-être qu’une analyse de complexité avec des ordres de grandeur pour montrer que ce n’est pas suffisant si on veut appliquer des filtres sur une vidéo en temps réel ou je ne sais quelle autre application pourrait-être assez efficace.

Autre point, même en connaissant déjà l’idée sous-jacente, j’avoue avoir eu un peu de mal à lire les schémas qui motivent l’utilisation de la permutation inverse de bits. Je sais que c’est toujours plus facile à critiquer qu’à faire les schémas efficaces, mais ça vaudrait peut-être le coup aussi de soit l’alléger / le décomposer, soit détailler un peu plus la légende (peut-être en donnant un exemple de lecture des différents éléments ?).

Petit détails plus formels :

  • il manque un chapeau au moment du calcul de f^[n+N/2]\hat{f}[n+N/2]​ pour trouver la formule de la "première FFT"
  • je trouve (mais c’est très personnel et absolument discutable) la notation des parties réelles et imaginaires avec des caractères gothiques un peu lourde visuellement, surtout dans des séquences un peu longues de calcul où ils interviennent à pas mal de reprises (Section "Calcul en place" pour le cas d’un signal réel notamment), j’ai l’impression que Re(xi[k])\text{Re}(x_i[k]) se lit plus facilement que R{xi[k]}\mathfrak{R}\left\{x_i[k]\right\}

Voilà, c’était en vrac quelques idées comme ça. Au delà de ça, c’était intéressant, même si je connaissais déjà une partie du contenu, j’ai quand même appris des choses, notamment sur les astuces d’implémentation et d’évaluation des fonctions trigonométriques sur des suites arithmétiques. Je suis même assez impressionné que l’implémentation en Julia arrive au niveau d’une implémentation de référence de FFTW. Tu as une idée des raisons qui font que ton implémentation arrive à être meilleure sur certaines entrées d’ailleurs ? Je suis assez curieux… ^^

+1 -0

Hello,

Je pense que le point évoqué par Næ est assez central. Considérer la notion de transformée de Fourier comme acquise est tout à fait acceptable vu que ce n’est pas le sujet de ce tutoriel. Mais pourquoi alors épiloguer sur la définition de celle discrète, si ce n’est pour introduire, soit des applications pratiques ou un contexte historique. Rien n’empêcherait de résoudre le problème de manière purement symbolique et de calculer le résultat de la transformée à la demande en suivant d’autres algorithmes d’évaluation d’intégrales par exemple.

Il faudrait peut-être parler à la fin des différences entre les algorithmes proposés, que ce soit l’officiel FFTW! ou le dernier FFT4. Peut-être que la perte de temps vient du "padding" rajouter par FFTW ou du changement de format dans les données (voir la ligne: real.(b[1:end÷2]) ≈ c[1:2:end] && imag.(b[1:end÷2]) ≈ c[2:2:end]). Quid de leur stabilité numérique, de leur erreur intrinsèque, de complexité asymptotique (ne pas répondre O(NlogN)O(N \log N) =) )

Sinon, voici les remarques purement textuelles:

Nous utiliserons la suivante : pour une fonction f, sa transformée de Fourier f chapeau est

J’aurais rajouté un: "défini par/comme suit:"

Comme on l’a vu, la transformée de Fourier est définie de manière continue.

Non, on ne l’a pas vu ensemble =)

Or, notre ordinateur ne dispose que d’une taille finie de mémoire, donc pour représenter un signal quelconque, nous ne pouvons utiliser qu’un nombre fini de valeurs.

L’implication du "donc" est hors-sujet. Il faudrait reformuler du style "Or, pour représenter un signal quelconque, notre ordinateur ne dispose que d’une taille finie de mémoire"

On échantillonne (ou discrétise) le signal à analyser, c’est à dire que l’on stocke la valeur du signal pour une suite de valeurs de sa variable.

J’essayerais d’éviter la confusion entre la valeur de la variable, et l’"évaluation" (la valeur ou un autre terme) de la transformée.

Avec \deltaδ la distribution de Dirac.

Éventuellement rajouter la définition de la fonction de Dirac.

g = ш_T X f

Personnellement, je ne suis pas fan de symbole "X", j’aurais préféré un "*".

Le choix du fenêtrage n’est absolument pas anodin, et peut mener à des problèmes innatendus si on l’ignore.

Inattendus

alors sa transformée de Fourier inverse se périodique de période 1/νs

Se périodique ? Ca m’a fait penser à "rotationer" =)

On peut donc implémenter relativement ce calcul !

Mais oui c’est clair !

Dans le code julia un symbole @.

Je ne sais pas lire Julia mais c’est quoi "@." alors qu’il y a une précision sur le "+/- égale"

Notre implémentation est donc vraiment lente et possède une empreinte mémoire très élevée par comparaison avec l’implémentation de référence ! Pour améliorer cela, nous allons implémenter la transformée de Fourier rapide.

Peut-être préciser que le facteur est 5000 ?

Cela qu’en calculant deux transformées de Fourier de longueur N/2,

Mais oui c’est clair !

cela dessine naturelement une implémentation récursive de la FFT.

naturellement

La seule subtilité est qu’il ne faut réaliser l’inversion qu’une fois par case

Cela m’a fait buggé le mot "case" parce qu’il s’agit du premier emploi du terme sans spécifier ce qu’il désigne.

(cela signifie que le compilateur place les quelques variables intermédiaires dans des registres du compilateur)

Quelques doutes sur cette affirmation. Il y a des chances que le stack - une partie de la mémoire - soit employé, non ?

Afin d’économiser de l’espace de stockage, on peut penser à utiliser cette moitié de tableau pour stocker noss nombres complexes.

Outre la type "nos" nombres complexes. La phrase juste avant précise qu’on s’intéresse (le plus souvent) aux signaux réels alors pourquoi conserver des nombres complexes ? La suite des phrases n’est pas cohérente.

Numerical Recipes on peut lire (section 5.4 "Recurrence Relations and Clenshaw’s Recurrence Formula", page 219 de la troisième édition)

Petite référence du livre =) Que ce soit dans l’introduction ou dans le texte.

my_fft_4(x)

Inversion des benchmarks par rapport à tous ceux présentés précédemment.

Im et Re

Je conforte ce que dit Næ, l’écriture gothique est illisible. J’ai du mal à interpréter la formule de récurrence des termes.

Bonjour les agrumes !

La bêta a été mise à jour et décante sa pulpe à l’adresse suivante :

Merci d’avance pour vos commentaires.


Merci à tous les deux pour vos retours, j’ai essayé de les incorporer au mieux !

Je suis même assez impressionné que l’implémentation en Julia arrive au niveau d’une implémentation de référence de FFTW. Tu as une idée des raisons qui font que ton implémentation arrive à être meilleure sur certaines entrées d’ailleurs ? Je suis assez curieux… ^^

Je ne connais pas en profondeur le code source de FFTW, mais comme ça j’ai trois hypothèses :

  1. (à mon avis la plus probable, en tout cas celles qui joue le plus) L’implémentation du reverse bit ordering ad-hoc dans mon code. En pratique cela constitue un très gros raccourci dont FFTW ne peut probablement pas savoir qu’il peut l’utiliser (condition sur la taille des tableaux etc);
  2. Je ne peut traiter que des tableaux dont la taille est une puissance de 2;
  3. FFTW est une bibliothèque externe en C que Julia doit appeler, il y a peut-être un très léger surcoût

Quelques doutes sur cette affirmation. Il y a des chances que le stack - une partie de la mémoire - soit employé, non ?

Gawaboumga

Oui tout à fait, j’avais encore la tête dans mon Arduino en écrivant ça je crois. ^^

J’essayerais d’éviter la confusion entre la valeur de la variable, et l’"évaluation" (la valeur ou un autre terme) de la transformée.

Gawaboumga

J’ai du mal à voir comment tourner ça. Je ne suis pas certain de comprendre si ce qui te dérange est que j’assimile la transformée de Fourier du signal à l’évaluation de cette transformée sur tous les points qui nous intéressent, ou si c’est plus simplement la tournure de ma phrase qui n’est pas claire ?

Personnellement, je ne suis pas fan de symbole "X", j’aurais préféré un "*".

Gawaboumga

Pour le coup je préfère ×\times, * est souvent associé au produit de convolution en physique et en traitement du signal.

Petite référence du livre =) Que ce soit dans l’introduction ou dans le texte.

Gawaboumga

Et… je viens de me rendre compte en rédigeant ce message que j’ai oublié de m’occuper de ça. Je le derais la modification demain.

Il faudrait peut-être parler à la fin des différences entre les algorithmes proposés, que ce soit l’officiel FFTW! ou le dernier FFT4. Peut-être que la perte de temps vient du "padding" rajouter par FFTW ou du changement de format dans les données (voir la ligne: real.(b[1:end÷2]) ≈ c[1:2:end] && imag.(b[1:end÷2]) ≈ c[2:2:end]). Quid de leur stabilité numérique, de leur erreur intrinsèque, de complexité asymptotique (ne pas répondre O(Nlog⁡N)O(N \log N)O(NlogN) =) )

Gawaboumga

Idem, j’ai oublié de faire la modification. Par contre je ne suis pas super à l’aise sur la deuxième partie de ta remarque, donc si tu as un peu de biblio vers laquelle m’orienter (ou même orienter le lecteur) je suis preneur.

Edit: Bon en fait il y a un papier qui parle du design de FFTW. Je vais essayer de parler de tout ça sans trop dire de bêtises (ça commence à s’éloigner de mon domaine d’expertise)

+0 -0

Hello,

Désolé, je n’ai fait que survoler les modifications que tu as apportées. Et, a priori, elles me conviennent !

J’ai du mal à voir comment tourner ça. Je ne suis pas certain de comprendre si ce qui te dérange est que j’assimile la transformée de Fourier du signal à l’évaluation de cette transformée sur tous les points qui nous intéressent, ou si c’est plus simplement la tournure de ma phrase qui n’est pas claire ?

Ce que je voulais exprimer, c’est que dans l’expression de la transformée de Fourier:

$$ \hat{f}(\nu) = \int_{-\infty}{+\infty} f(x) e{-2i\pi \nu x} dx $$

Tu as deux variables xx (dépendante) et ν\nu (indépendante), mais tu parles également de "valeur" pour désigner f^(ν)\hat{f}(\nu). Cf:

c’est à dire que l’on stocke la valeur du signal pour une suite de valeurs de sa variable. Dans le cas de la FFT, on échantillonne avec un pas constant. Par exemple si on regarde un signal temporel comme la valeur d’une tension lue sur un voltmètre, on pourait enregistrer la valeur à chaque tic d’une montre.

Donc, quand tu parles de la "valeur du signal pour des valeurs de sa variable", je ne trouve vraiment pas cela clair. Surtout que DFT permet de calculer la "représentation spectrale discrète du signal échantillonné", on ne s’intéresse pas au "signal" à proprement parler mais à son spectre. Quitte à remplacer "de sa variable" par ν\nu ou parler de l’évaluation du spectre du signal. Je ne sais pas si c’est plus clair.

PS: on pourrait.

Pour le coup je préfère ×\times, *∗ est souvent associé au produit de convolution en physique et en traitement du signal.

Tu as parfaitement raison, j’avais oublié la notation du produit de convolution !

Par contre je ne suis pas super à l’aise sur la deuxième partie de ta remarque, donc si tu as un peu de biblio vers laquelle m’orienter (ou même orienter le lecteur) je suis preneur.

Je dois admettre ne pas pouvoir aider davantage. Les algorithmes que tu as proposés et celui officiel sont tous différents. En ce sens, qu’au delà de la vitesse d’exécution ou de la consommation mémoire, tu évoquais à juste titre que:

Je ne peut traiter que des tableaux dont la taille est une puissance de 2;

Cela me semble très important comme différence à mentionner ! Je me demandais également si les résultats étaient "strictement égaux" pour une large variété de fonctions et pas seulement ton exemple. La fonction rand de Julia doit sans doute retourner des nombres entre 0 et 1 qui ont l’avantage de généralement très bien se comporter en arithmétique floatante. Chaque opération sur des floats peut introduire des erreurs et certaines peuvent en introduire plus que d’autres. C’était le sujet de ma réflexion mais si personne n’a le niveau pour répondre précisément à ce genre de questions, tant pis …


Ca m’a également fait tilt, mais éventuellement mentionner la matrice de Vandermonde avec la forme matricielle ?

Bonjour les agrumes !

La bêta a été mise à jour et décante sa pulpe à l’adresse suivante :

Merci d’avance pour vos commentaires.


Et voici les dernières remarques intégrées. :) S’il n’y a pas d’autres remarques d’ici à vendredi soir je l’enverrais en validation.

Ca m’a également fait tilt, mais éventuellement mentionner la matrice de Vandermonde avec la forme matricielle ?

Gawaboumga

Oui tout à fait ! D’ailleurs cette matrice de Vandermonde permet de visualiser la FFT comme un algorithme qui peut servir à multiplier deux polynômes. J’ai mis un lien vers cette vidéo que j’aime beaucoup dans le tutoriel.

+0 -0

Salut,

Je n’ai pas lu en détails car je n’ai pas le niveau mathématique suffisant, néanmoins quelques remarques :

  • A qui s’adresse cet article ? Quel est le public visé ?
  • Retravailler la forme pour que ce soit un peu plus aéré, plus facile à lire

Pour le premier point, ça serait pas mal de mettre une indication dans l’introduction. Là comme ça, quand je survole l’article, je vois plein d’équations à base de sigma, d’intégrales, pfiiuu ça fait un peu peur :-°
De même il y a quelques blocs de code avec de gros pavés textuels. On rejoint là la remarque sur la forme.

Je comprend que l’objectif n’est pas de faire une présentation détaillée de la FFT, sujet assez traité sur le net, mais ne mettre que des équations et aucun graphique alors qu’on parle d’échantillonnage, je trouve ça vraiment dommage.
C’est probablement un parti pris de rester dans la partie mathématique, mais les applications de la FFT sont concrètes, par exemple dans les oscilloscopes.

L’article est relativement long, mais c’est en partie dû aux nombreuses équations et blocs de code. Néanmoins mettre cet article en tutoriel avec deux parties distinctes : une première où tu introduit la transformée de Fourrier et la FFT au niveau mathématique et une seconde sur la partie code et optimisation; ça ferait quelque chose de plus clair et surtout enlèverait une sorte de poids sur la longueur de l’article.

Salut,

Merci pour ton retour, a priori le public visé a déjà utilisé une transformée de Fourier, sans pour autant l’avoir déjà implémenté. Je comprends les remarques sur la forme et j’ai quelques idées pour rendre ça plus visuel, surtout la première partie. Je peux potentiellement faire des GIF animés aussi, je verrais ce qui rend le mieux !

Connectez-vous pour pouvoir poster un message.
Connexion

Pas encore membre ?

Créez un compte en une minute pour profiter pleinement de toutes les fonctionnalités de Zeste de Savoir. Ici, tout est gratuit et sans publicité.
Créer un compte