🦐 😫 👩🏻‍🔧 Transformation de Fourier. Le rapide et le furieux 👨‍👩‍👧‍👦 ✋🏿 🖕🏻

Souvent, lors du développement d'algorithmes, nous nous heurtons à la limite de la complexité de calcul, qui, semble-t-il, est impossible à surmonter. La transformée de Fourier a de la complexité

O (n^{2})

$O (n ^ 2)$ , et une version rapide, proposée vers 1805 par House ¹ (et réinventée en 1965 par James Cooley et John Tukey)

O (n l o g (n))

$O (nlog (n))$ . Dans cet article, je veux vous montrer que vous pouvez obtenir les résultats de la conversion en temps linéaire

O (n)

$O (n)$ ou même atteindre une difficulté constante

O (1)

$O (1)$ sous certaines conditions qui se trouvent dans de vrais problèmes.
Analyse de Fourier

Lorsque j'ai été confronté à la tâche d'écrire un programme pour analyser les fonctions de transfert des systèmes de son en temps réel, je me suis d'abord tourné vers la conversion rapide comme tout le monde. Tout allait bien, mais avec la grande taille de la fenêtre de temps, la charge du processeur est devenue indécemment importante et quelque chose devait être fait. Il a été décidé de faire une pause et d'étudier à nouveau la transformation, tout en cherchant des moyens de résoudre le problème. Revenons à la transformation originale de Joseph Fourier ² :

f (x) = s u m l i m i t s_{- i n f t y}^{+ i n f t y} c_{k} e^{2 p i i k x / T} c_{k} = f r a c 1 T i n t l i m i t s_{0}^{T} f (x) e^{- 2 p i i k x / T} d x

$f (x) = \ sum \ limits _ {- \ infty} ^ {+ \ infty} c_ke ^ {2 \ pi ikx / T} \\ c_k = \ frac {1} {T} \ int \ limits_0 ^ Tf ( x) e ^ {- 2 \ pi ikx / T} dx$

Nous examinerons attentivement ce qui se passe ici. Chaque valeur de sortie dans le domaine fréquentiel

c_{k}

$c_k$ est la somme de toutes les valeurs d'entrée du signal

f (x)

$f (x)$ multiplié par

e^{- 2 p i i k x / T}

$e ^ {- 2 \ pi ikx / T}$ . Pour effectuer des calculs, nous devons parcourir toutes les données d'entrée pour chaque valeur de sortie, c'est-à-dire pour remplir ces

n^{2}

$n ^ 2$ opérations.

Débarrassez-vous de n

Permettez-moi de vous rappeler qu'au départ, la tâche consistait à analyser les données sonores en temps réel. Pour ce faire, la fenêtre temporelle sélectionnée (essentiellement un tampon) de taille N est remplie de données avec une fréquence f _d correspondant à la fréquence d'échantillonnage. Avec la période T, les données d'entrée sont converties de la fenêtre temporelle en fenêtre fréquentielle. Si vous regardez les nombres réels, alors N varie de 2 ¹⁴ (16 384) à 2 ¹⁶ (65 536) échantillons (les valeurs sont héritées de la FFT, où la taille de la fenêtre doit être une puissance de deux). Temps T = 80 ms (12,5 fois par seconde), ce qui vous permet de voir les changements très facilement et de ne pas surcharger le CPU et le GPU. La fréquence d'échantillonnage f _{d est} standard et est de 48 kHz. Calculons la quantité de données dans la fenêtre de temps change entre les dimensions. Pendant le temps T, il entre dans le tampon

$inline$ échantillons. Ainsi, seulement 5% à 23% des données sont mises à jour dans la fenêtre. Dans le pire des cas, 95% (et au mieux 73%, ce qui est aussi beaucoup!) Des échantillons traités tomberont encore et encore dans la conversion, malgré le fait qu'ils ont déjà été traités dans les itérations précédentes.Le lecteur attentif à ce moment lèvera la main et dira: "attendez, mais qu'en est-il du coefficient

e^{- 2 p i i k x / T}

$e ^ {- 2 \ pi ikx / T}$ ? Après tout, à chaque nouvelle transformation, les mêmes données seront localisées aux nouvelles positions de la série et, par conséquent, auront des coefficients différents? » Pour chaque cinq pour leurs soins, rappelons un détail important de la transformation qui est souvent oublié. Dans l'étude des valeurs de fonction

f (t)

$f (t)$ sur l'intervalle de 0 à t, la fonction est considérée comme périodique, ce qui vous permet de décaler sans douleur la fonction vers la gauche ou la droite dans le temps. Par conséquent, nous avons le droit de ne pas insérer une nouvelle valeur à la fin et de supprimer l'ancienne valeur depuis le début, mais de remplacer cycliquement les données dans le tampon.

Pour plus de clarté, vous pouvez écrire sous forme de tableau comment le tampon va changer:

t = 0	f (0)	f (1)	f (2)	f (3)	f (4)	f (5)	f (6)	f (7)	f (8)	f (9)
t = 1	f (10)	f (1)	f (2)	f (3)	f (4)	f (5)	f (6)	f (7)	f (8)	f (9)
t = 2	f (10)	f (11)	f (2)	f (3)	f (4)	f (5)	f (6)	f (7)	f (8)	f (9)
t = 3	f (10)	f (11)	f (12)	f (3)	f (4)	f (5)	f (6)	f (7)	f (8)	f (9)
t = 4	f (10)	f (11)	f (12)	f (13)	f (4)	f (5)	f (6)	f (7)	f (8)	f (9)

Vous pouvez écrire comment la transformation dans le temps passe de t ₁ à t ₂ :

F_{t} = F_{t - 1} + D e l t a F D e l t a F : D e l t a c_{k} = f r a c 1 T i n t l i m i t s_{t_{1}}^{t_{2}} (f_{t} (x) - f_{t - 1} (x)) e^{- 2 p i i k x / T} d x

$F_t = F_ {t-1} + \ Delta F \\ \ Delta F: \ Delta c_k = \ frac {1} {T} \ int \ limits_ {t_1} ^ {t_2} (f_t (x) - f_ { t-1} (x)) e ^ {- 2 \ pi ikx / T} dx$

Valeur

F_{t - 1} (x)

$F_ {t-1} (x)$ est le résultat de la conversion précédente et la complexité du calcul

D e l t a f (x)

$\ Delta f (x)$ ne dépend pas de la taille de la fenêtre temporelle et est donc constant. En conséquence, la complexité de la conversion sera

O (n)

$O (n)$ ^* car il ne nous reste plus qu'à parcourir une fois la fenêtre de fréquence et à appliquer les changements pour les échantillons T qui ont changé au cours du temps. Je voudrais également attirer votre attention sur le fait que les probabilités

e^{- 2 p i i k x / T}

$e ^ {- 2 \ pi ikx / T}$ peut être calculé à l'avance, ce qui donne un gain supplémentaire de productivité, et il ne reste que deux opérations dans le cycle: soustraire des nombres réels et multiplier un nombre réel par un complexe, dans la pratique, ces deux opérations sont simples et bon marché.

Pour compléter le tableau, il ne reste plus qu'à indiquer l'état initial, mais ici tout est simple:

F_{0} (x) = 0

$F_0 (x) = 0$

* - bien sûr, la complexité finale de toute la transformation restera ainsi

O (n^{2})

$O (n ^ 2)$ , mais il sera exécuté progressivement, sur n itérations, pendant la mise à jour du tampon.

O (n)

$O (n)$ - c'est la complexité de la mise à jour des données, mais c'est exactement ce dont nous avons besoin (lors de l'utilisation de la FFT, la complexité de chaque transformation

O (n l o g (n))

$O (nlog (n))$ )

Mais que faire si vous creusez plus profondément. Ou se débarrasser du deuxième n

Je veux tout de suite faire une réservation pour que les prochaines étapes ne soient applicables que si vous ne prévoyez pas d'effectuer la transformation inverse pour le résultat (afin de corriger le signal ou d'obtenir une réponse impulsionnelle). Pour commencer, je tiens à vous rappeler qu'à la suite de la conversion, nous obtenons un tableau de données symétrique, ce qui nous permet immédiatement de réduire de moitié le nombre de conversions.

Analysons maintenant l'ensemble de données résultant, compte tenu des conditions du problème. Nous avons un ensemble de nombres complexes, chacun décrivant l'amplitude et la phase des oscillations à une fréquence particulière. La fréquence peut être déterminée par la formule:

f [j] = j f r a c f d N

$f [j] = j \ frac {fd} {N}$ pour

j < f r a c N 2

$j <\ frac {N} {2}$ . Évaluons l'étape de la fenêtre de fréquence sur nos données:

D e l t a f = f r a c f d N

$\ Delta f = \ frac {fd} {N}$ Pour N = 2 ¹⁴ : 2,93 Hz (et pour 2 ¹⁶ : 0,73 Hz). Ainsi, dans la plage de 1 kHz à 2 kHz, nous obtenons 341 résultats. Essayez d'évaluer indépendamment la quantité de données dans la plage de 8 kHz à 16 kHz pour N = 65536. Beaucoup, non? Beaucoup! Avons-nous besoin de tant de données? Bien sûr, dans les problèmes d'affichage des caractéristiques de fréquence des systèmes sonores, la réponse est non. Et d'autre part, pour l'analyse dans la région des basses fréquences, un petit pas est très utile. N'oubliez pas qu'il y a encore un calendrier à prévoir pour ces volumes (

f r a c N 2

$\ frac {N} {2}$ ) convertir en une forme lisible par l'homme (moyennage, spline ou lissage) et les afficher à l'écran. Et à des fréquences élevées, même avec un écran 4K et l'affichage du graphique en mode plein écran avec l'axe de fréquence logarithmique, la taille du pas se révélera rapidement bien inférieure à 1 pixel.

Par expérience, vous pouvez découvrir qu'il suffit d'avoir seulement 48 points par octave, et pour avoir les données un peu plus lisses et plus moyennes, je suggère de s'arrêter à 96. Dans la plage de fréquences audio de 20 Hz à 20 kHz, il est facile de compter seulement 10 octaves: 20, 40, 80 , 160, 320, 640, 1280, 2560, 5120, 10240, 20480, chacune pouvant être divisée en un nombre donné de sous-bandes (n'oubliez pas que la partition doit être faite géométriquement et non arithmétiquement), par conséquent, il est plus que suffisant d'effectuer la conversion uniquement pour 960 fréquences pour obtenir Performan que dans 16 ... 65 fois plus petit que la version originale.

Ainsi, en combinant les deux approches, nous obtenons la complexité constante de l'algorithme de mise à jour des données

O (1)

$O (1)$ .

Miel carré et une cuillerée de goudron

Maintenant, nous pouvons dire que la complexité

O (n^{2})

$O (n ^ 2)$ nous sommes arrivés à la complexité

O (1)

$O (1)$ en utilisant deux astuces de vie simples:

Après avoir analysé le problème, nous avons remarqué que les données sont ajoutées progressivement et que la période de mise à jour complète de la fenêtre temporelle est beaucoup plus élevée que la période de transformations et nous avons ensuite calculé la différence de la transformée de Fourier.
passé de l'étape arithmétique dans la fenêtre de fréquence à limité uniquement par les valeurs spécifiées, ce qui peut réduire considérablement le nombre de conversions.

Mais, bien sûr, la vie serait vraiment un conte de fées, sinon un mais. L'application de ces deux approches nous a permis de vraiment décharger le CPU pour que deviner qu'il calcule la transformée de Fourier et affiche les résultats à l'écran même avec

N = 2^{20}

$N = 2 ^ {20}$ c'était difficile. Mais la punition ne s'est pas fait attendre lorsque vos signaux en réalité ne sont pas périodiques (et cela est nécessaire pour obtenir les résultats de conversion corrects) et qu'il n'est pas possible de sélectionner la taille de fenêtre appropriée, il devient nécessaire d'utiliser diverses fonctions de fenêtre, ce qui ne vous permet plus d'utiliser pleinement la première étape. La pratique a montré que l'utilisation des fonctions de fenêtre est critique dans l'étude des signaux avec une fréquence inférieure à

0, 1 f_{d}

$0,1f_d$ . Aux hautes fréquences, le nombre de périodes tombant dans la fenêtre temporelle atténue significativement les distorsions résultant de la présence d'un écart de premier ordre (entre f (0) et f (N-1)) dans la fonction d'origine.

Au final, j'ai également refusé la deuxième étape et je suis retourné à la FFT, car le gain dans cette tâche était déjà faible.

En conclusion

La première approche peut être appliquée si vos données sont de nature périodique prononcée et doivent être analysées dans le temps à l'aide d'une grande fenêtre de temps, qui, je le rappelle, n'a pas besoin d'être de degré 2, c'est-à-dire tout nombre naturel.
La deuxième approche est applicable (même en tenant compte des fonctions de fenêtre) si seul un petit ensemble de fréquences est analysé dans les données.

Hélas, pour moi dans ce problème, cela n'est resté qu'un petit divertissement mathématique, mais j'espère que cela vous inspirera à étudier d'autres algorithmes en vacances en termes de changements dans les données d'entrée au fil du temps :)

Littérature

Image tirée du manga de Michio Shibuya. “MATHÉMATIQUES EXCITANTES. Analyse de Fourier

Transformation de Fourier. Le rapide et le furieux

Débarrassez-vous de n

Mais que faire si vous creusez plus profondément. Ou se débarrasser du deuxième n

Miel carré et une cuillerée de goudron

En conclusion

Littérature

More articles: