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

Tout le monde se secoue ! 

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.

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 $\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 $\text{Re}(x_i[k])$ se lit plus facilement que $\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…

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(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 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

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

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 ?

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

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.

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)

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 $x$ (dépendante) et $\nu$ (indépendante), mais tu parles également de "valeur" pour désigner $\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.

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.

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 !

Bonjour les agrumes !

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

Merci d’avance pour vos commentaires.

Salut,

Pour moi le tutoriel est clair, et je n’ai pas remarqué d’erreur de fond.

J’ai tout de même une question : dans la section optimisation des fonctions trigonométriques, tu donnes une méthode de calcul par récurrence en précisant que "Cette relation présente également des intérêts en terme de stabilité numérique." En quoi est-ce avantageux, par rapport à la relation "naïve", que tu utilises pour montrer la validité de cette relation ? Est-ce parce qu’on vient uniquement ajouter ou soustraire des petits nombres, au lieu de commencer par une multiplication ?

Par ailleurs, j’appuie la remarque de @zeqL concernant le public visé : ce serait bien de mettre une petite phrase dans l’introduction pour préciser ce que tu nous a dit ici.

Sinon, pour avoir suivi rapidement l’historique, je trouve que tu as nettement amélioré la structure et la lisibilité, ce qui est une bonne chose !

Salut, merci pour le retour !

je viens de me rendre compte que j’ai oublié d’intégrer la remarque sur le public visé dans l’intro…

Je plaide coupable d’avoir un peu pris Numerical Recipes au mot. Texto, voici ce qu’ils indiquent :

Donc de ce que j’en comprend, si on utilise la relation bien connue de $\cos(a+b)$ il y a des chances que l’on trouve des résultats aberrants pour les petites valeurs de $\delta$. Du coup je me suis demandé à quel point c’était vraiment critique, et pour des incréments de $2\pi/1024$ il faut vraiment le chercher pour faire dérailler la méthode naive, mais on finit quand meme par y arriver en utilisant des flottants demi précision (Float16) autour des points où l’une des deux fonctions s’annulent :

Je ne sais pas si ça vaut le coup d’insister dessus dans le tutoriel par contre ?

Le code pour générer la figure pour celles et ceux qui veulent jouer :

using Plots, LaTeXStrings,

background_color = RGB(250/255, 250/255, 250/255)
blue = RGB(11/255, 90/255, 127/255)
orange = RGB(248/255, 171/255, 48/255)

n = 2^10
T = Float16

cosines_nr, sines_nr = let
    res_cos = zeros(T, n)
    res_sin = zeros(T, n)
    δ = T(2π/n)
    cj, sj = T(1),T(0)
    α = 2sin(δ/2)^2
    β = sin(δ)
    for j ∈ 1:n
        res_cos[j] = cj
        res_sin[j] = sj
        cj, sj = cj - (α*cj + β*sj), sj - (α*sj-β*cj)
    end
    res_cos, res_sin
end

cosines_simple, sines_simple = let
    res_cos = zeros(T, n)
    res_sin = zeros(T, n)
    δ = T(2π/n)
    cj, sj = T(1),T(0)
    α = cos(δ)
    β = sin(δ)
    for j ∈ 1:n
        res_cos[j] = cj
        res_sin[j] = sj
        cj, sj = cj*α - sj*β, cj*β + sj*α
    end
    res_cos, res_sin
end

cosines_true, sines_true = begin
    δ = 2π/n
    θ = range(0, length=n, step=δ)
    cos.(θ), sin.(θ)
end

plot_sines = begin
  plot(
    θ, 
    [cosines_true cosines_simple cosines_nr sines_true sines_simple sines_nr],
    label=["True cosine" "Naive cosine" "NR cosine" "True sine" "Naive sine" "NR sine" ],
    linestyle=[:solid :dash :dot :solid :dash :dot],
    color=[blue blue blue orange orange orange],
    background_color=background_color, 
    border=:box,
    xlabel=L"\theta",
  )
  lens!([π-0.5, π+0.5], [-1.1, -0.9], inset=(1, bbox(0.5, 0.1, 0.4, 0.4)), border=:box)
end

savefig(plot_sines, "plot_trigo.png")

Hello,

J’ai été étonné par le nombre de changements apportés en si peu de temps, et pour le mieux !

Plus grand chose à dire, alors je mets juste les fautes que j’ai repérées.

il peut être intéressant d’implémenter et d'optimiser sa propre FFT. Une application qui a peu été discutée dans ce tutoriel est le calcul de produits de convolution. Une méthode efficace dans le cas où on convolue des signaux de longueurs comparable est

Dans cette première partie, je vous propose de découvrir comment construire la transformée de Fourier discrète puis comprendre pourquoi la transformée de FOurier rapide est utile.

Nous utiliserons la suivante : (étant donné ?) pour une fonction f, sa transformée de Fourier \hat{f} est définie par:

De la transformée de Fourier à la transformée de Fourier discrète

Préciser qu’elle est continue ?

Par exemple, si on regarde un signal temporel comme la valeur d’une tension lue sur un voltmètre, on pourrait enregistrer la valeur à chaque tic d’une montre.

On échantillonne la transformée de Fourier du signal pour obtenir la transformée de FOurier discrète.

Avec δ la distribution de Dirac. Voici le graphe que l’on peut obtenir si l’on représente f et g ensemble.

g est défini juste après.

les cloches de la transformée du signal échantillonné se supperposent.

Il est clair que l’on souhaite que l’échantillonnage soit le plus "fin" possible, pour ne manque aucun détail de la transformée de Fourier** **!

De manière générale, je ne suis pas convaincu que échantillonn* prennent bien 2 n partout dans le texte.

Or, la transformée de Fourier inverse (l’opération qui permet de retrouver le signal à partir de sa transformée de Fourier)

Il faudrait ici justifier que l’on peut inverser les symboles signes somme et intégrale

n effet la fonction exp et les opérateurs de division et multiplication sont définis pour des scalaires. Cette macro permet d’informer Julia qu’il doit appliquer les opérations scalaires termes à terme.

La transformée de Fourier a d’abord une foultitude d’applications théoriques, que ce soit pour résoudre des équations différentielles, en traitement du signal ou en physique quantique. Elle possède également des applications pratique en optique et en spectroscopie.

100 000 000 000 octets de mémoire préciser 100Go ? (un nombre complexe étant est stocké sur 2 flottants)

637espace 537 millisecondes, soit plus de 10 minutesespace!

De manière générale, j’ai l’impression que ton éditeur de texte est en anglais et colle le point d’exclamation au dernier mot de la phrase.

D’après notre benchmark, l’algorithme traite une entrée de 1024 points en 23.785µms.

Au plus, on a donc log_2(N) tableaux alloués avec des tailles divisées par deux à chaque fois.

Je sais qu’il y a débat sur l’usage de la virgule après des adverbes dits "charnières" ou des conjonctions de coordination. En fonction de s’ils sont inclus dans le groupement des mots qui suit.

Autrement dit, cela donne l’idée que l’on pourrait s’épargner toutes ces allocations de tableaux

Il y a deux petites subtilités dues à Julia : les tableaux commencent leur numérotation à 1, et on utilise la macro @inbounds pour accélérer un peu le code en désactivant les vérifications de débordement de tableaux.

Pourquoi ce n’est pas dans un encart warning comme tes précédentes remarques sur le plus ou moins égale ou @. ?

Le tableau h, qui est présenté précédemment, est à valeurs complexes

Merci pour la réponse !

Je ne pense pas non plus qu’il faille trop s’y étendre dans le tutoriel. Maintenant, tu as le recul et les données qui confirment l’affirmation du bouquin, et si des gens se posent des questions là-dessus il sera facile de les rediriger sur ton post.

J’allais relever quelques typos mais @Gawaboumga l’a fait mieux que moi pendant que je rédigeais.

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 pour vos retours, je pense que j’ai tout intégré, sauf la précision sur la continuité vu que c’est fait juste après, et les points d’exclamation parce que j’ai bien un espace avant dans le code Markdown, donc je suppose que c’est l’outil de rendu de Zeste de Savoir qui les retire.

J’ai quelques remarques à te faire sur la forme, qui sont détails avant de passer à la validation :

• la sous-partie Résumé des performances ressemble pas mal à une conclusion, et par ailleurs tu n’as pas de conclusion qui suit : je te suggère donc de passer tout ou partie de cette sous-partie dans la conclusion de la partie;
• le tutoriel n’a pas de miniature, or c’est obligatoire : tu devrais pouvoir reprendre une image de ton tuto en la recadrant un peu;
• je pense que d’autres sections (comme Programmation et algorithmique) seraient adaptées à ton tuto.

C’est à peu près tout pour moi.

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 voilà ! J’espère que je n’ai pas trop spammé l’équipe de validation en mettant à jour le contenu ^^'

Bonjour,

La bêta du contenu « Jouons à implémenter une transformée de Fourier rapide ! » a été désactivée (tuto publié).