Quelles sont les différences entre la DFT et la FFT qui rendent la FFT si rapide?

16

J'essaie de comprendre les FFT, voici ce que j'ai jusqu'à présent:

Afin de trouver l'amplitude des fréquences dans une forme d'onde, il faut les rechercher en multipliant l'onde par la fréquence qu'ils recherchent, en deux phases différentes (sin et cos) et en faisant la moyenne de chacune. La phase se trouve par sa relation aux deux, et le code pour cela est quelque chose comme ceci:

//simple pseudocode
var wave = [...];                //an array of floats representing amplitude of wave
var numSamples = wave.length;
var spectrum = [1,2,3,4,5,6...]  //all frequencies being tested for.  

function getMagnitudesOfSpectrum() {
   var magnitudesOut = [];
   var phasesOut = [];

   for(freq in spectrum) {
       var magnitudeSin = 0;
       var magnitudeCos = 0;

       for(sample in numSamples) {
          magnitudeSin += amplitudeSinAt(sample, freq) * wave[sample];
          magnitudeCos += amplitudeCosAt(sample, freq) * wave[sample];
       }

       magnitudesOut[freq] = (magnitudeSin + magnitudeCos)/numSamples;
       phasesOut[freq] = //based off magnitudeSin and magnitudeCos
   }

   return magnitudesOut and phasesOut;
}

Pour ce faire pour de très nombreuses fréquences très rapidement, les FFT utilisent de nombreuses astuces.

Quelles sont certaines des astuces utilisées pour rendre les FFT beaucoup plus rapides que les DFT?

PS J'ai essayé de regarder les algorithmes FFT terminés sur le Web, mais toutes les astuces ont tendance à être condensées en un seul beau morceau de code sans trop d'explications. Ce dont j'ai besoin en premier, avant de pouvoir comprendre le tout, c'est une introduction à chacun de ces changements efficaces en tant que concepts.

Je vous remercie.

Seph Reed
la source
7
"DFT" ne fait pas référence à un algorithme: il fait référence à une opération mathématique. "FFT" se réfère à une classe de méthodes pour calculer cette opération.
1
Je voulais juste souligner que l'utilisation de sudodans votre exemple de code pourrait prêter à confusion, car il s'agit d'une commande bien connue dans le monde informatique. Vous vouliez probablement dire psuedocode.
rwfeather
1
@nwfeather Il voulait probablement dire «pseudocode».
user207421

Réponses:

20

L'implémentation naïve d'une TFD à points est fondamentalement une multiplication par une matrice N × N. Il en résulte une complexité de O ( N 2 ) .NN×NO(N2)

L'un des algorithmes de FFT (Fast Fourier Transform) les plus courants est l'algorithme FFT de décimation dans le temps de Radix-2 Cooley-Tukey. Il s'agit d'une approche de base diviser pour mieux régner.

Définissez d'abord le "facteur de twiddle" comme: j

WNej2πN
est l'unité imaginaire, alors la DFTX[k]dex[n]est donnée par X[k]= N - 1 n = 0 x[j1X[k]x[n] Si N est pair (et N
X[k]=n=0N1x[n]WNkn.
N est un entier), la somme peut alors être divisée en deux sommes comme suit X[k]= N / 2 - 1 n=0x[2n]W 2 k n N + N / 2 - 1 n=0x[2n+1]W k ( 2 n + 1 ) NN2
X[k]=n=0N/21x[2n]WN2kn+n=0N/21x[2n+1]WNk(2n+1)
où la première somme traite des échantillons pairs de et la seconde des échantillons impairs de x [ n x o [ n ] x [ 2 n + 1 ] et utiliser le fait quex[n] . Définissant x e [ n ] x [ 2 n ] etx[n]xe[n]x[2n]xo[n]x[2n+1]
  1. , etWNk(2n+1)=WN2knWNk
  2. WN2kn=WN/2kn

X[k]=n=0N/21xe[n]WN/2kn+WNkn=0N/21xo[n]WN/2kn=Xe[k]+WNkXo[k]
Xe[k]Xo[k]N2x[n]NN2
2(N2)2+N<N2
N>2

O(NlogN)O(N2)

anpar
la source
seriez-vous prêt à énumérer ce que représente chacune des variables? Je suis plutôt nouveau à cela, donc W, j, X(), Net kn'ont pas encore des définitions pour moi.
Seph Reed
Wkn
anpar
19

http://nbviewer.jupyter.org/gist/leftaroundabout/83df89a7d3bdc24373ea470fb50be629

DFT, taille 16

Diagramme des opérations dans un DFT naïf de taille 16

FFT, taille 16

Diagramme des opérations dans une FFT radix-2 de taille 16

La différence de complexité est assez évidente, n'est-ce pas?


Voici comment je comprends la FFT.

FT:L2(R)L2(R)

RC mais uniquement sur l'espace des fonctions intégrables (Lebesgue-, square-) . Maintenant, cette intégrabilité n'est pas une propriété très forte (beaucoup plus faible que la différentiabilité, etc.), mais elle exige que la fonction devienne «localement discribable avec des informations dénombrables». Une telle description est donnée par les coefficients d'un transformée de Fourier à court terme . Le cas le plus simple est que votre fonction est continue et que vous la divisez en si petites régions qu'elle est fondamentalement constante dans chacune d'elles. Ensuite, chacun des STFT a le plus fortement un terme nul. Si vous ignorez les autres coefficients (de toute façon en décomposition), chaque domaine n'est qu'un seul point de données. De tous ces coefficients de limite de temps court-LF, vous pouvez prendre une transformée de Fourier discrète. En fait, c'est exactement ce que vous faites lorsque vous effectuez un FT sur des données réelles mesurées!

Cependant, les données mesurées ne doivent pas nécessairement correspondre à une quantité physique fondamentale. Par exemple, lorsque vous mesurez une certaine intensité lumineuse , vous mesurez simplement amplitude d'une onde électromagnétique dont la fréquence est elle-même trop élevée pour être échantillonnée avec un ADC. Mais il est clair que vous pouvez également calculer la DFT d'un signal d'intensité lumineuse échantillonné, et à moindre coût, malgré la fréquence insensée de l'onde lumineuse.

Cela pourrait être compris comme la raison la plus importante pour laquelle la FFT est bon marché:

Ne vous embêtez pas à essayer de voir les cycles d'oscillation individuels du plus haut niveau. Au lieu de cela, ne transformez que des informations quelque peu élevées qui ont déjà été prétraitées localement.

Mais ce n'est pas tout. La grande chose au sujet de la FFT est qu'elle vous donne toujours toutes les informations qu'une DFT complète donnerait . C'est-à-dire toutes les informations que vous obtiendriez également lors de l'échantillonnage de l'onde électromagnétique exacte d'un faisceau lumineux. Cela peut-il être accompli en transformant un signal de photodiode? - pouvez-vous mesurer la fréquence lumineuse exacte à partir de cela?


Δν=1/Δt

En ayant globalement une durée plus longue, nous devrions également pouvoir réduire l'incertitude de fréquence. Et cela est en effet possible, si vous mesurez localement non seulement la fréquence approximative mais aussi la phase de l'onde. Vous savez qu'un signal de 1000 Hz aura exactement la même phase si vous le regardez une seconde plus tard. Alors qu'un signal de 1000,5 Hz, tout en étant indiscernable à petite échelle, aura une phase inversée une seconde plus tard.

Heureusement, ces informations de phase peuvent très bien être stockées dans un seul numéro complexe. Et c'est ainsi que fonctionne la FFT! Cela commence avec beaucoup de petites transformations locales. Celles-ci sont bon marché - pour une chose évidemment parce qu'elles n'utilisent qu'une petite quantité de données, mais deuxièmement parce qu'elles savent qu'en raison du court laps de temps, elles ne peuvent pas résoudre la fréquence de manière très précise de toute façon - donc c'est toujours abordable même si vous faire beaucoup de ces transformations.

Ceux-ci, cependant, enregistrent également la phase , et à partir de là, vous pouvez ensuite rendre la résolution de fréquence plus exacte au niveau supérieur. La transformation requise est à nouveau bon marché, car elle ne se soucie pas elle-même des oscillations haute fréquence, mais uniquement des données basse fréquence prétraitées.


Ouaip, mon argumentation est un peu circulaire à ce stade. Appelons ça récursif et tout va bien ...

Cette relation n'est pas mécanique quantique, mais l' incertitude de Heisenberg a en fait la même raison fondamentale.

à gauche
la source
2
belle représentation picturale du problème. :-)
robert bristow-johnson
2
N'aimez-
1
J'ai compris l'image après avoir lu la réponse d'Anpar.
JDługosz
15

WNnkej2πnkN

Notez le chemin indiqué et l'équation ci-dessous montre le résultat pour la fréquence bin X (1), tel que donné par l'équation de Robert.

Les lignes pointillées ne sont pas différentes des lignes pleines juste pour indiquer clairement où se trouvent les jointures de sommation.

Mise en œuvre de la FFT

Dan Boschen
la source
8

essentiellement, dans le calcul de la DFT naïve directement à partir de la somme:

X[k]=n=0N1x[n]ej2πnkN

Nej2πnkNNN1X[k]kX[k+1]

  1. de sorte que la FFT conserve certaines données intermédiaires.
  2. la FFT utilisera également un peu le facteur twiddle afin que le même facteur puisse être utilisé pour une combinaison intermédiaire de données.
robert bristow-johnson
la source
4

Je suis une personne visuelle. Je préfère imaginer la FFT comme une astuce matricielle plutôt que comme une astuce de sommation.

Pour expliquer à un niveau élevé:

Un DFT naïf calcule chaque échantillon de sortie indépendamment et utilise chaque échantillon d'entrée dans chaque calcul (algorithme N² classique).

Une FFT commune utilise des symétries et des motifs dans la définition de la DFT pour effectuer le calcul en "couches" (couches log N), chaque couche avec une exigence de temps constant par échantillon créant un algorithme N log N.

Plus de détails:

Une façon de visualiser ces symétries est de considérer la DFT comme une entrée matricielle 1 × N multipliée par une matrice NxN de toutes vos exponentielles complexes. Commençons par le cas "radix 2". Nous allons diviser les lignes paires et impaires de la matrice (correspondant aux échantillons d'entrée pairs et impairs) et les considérer comme deux multiplications de matrice distinctes qui s'ajoutent pour obtenir le même résultat final.

Regardez maintenant ces matrices: dans la première, la moitié gauche est identique à la moitié droite. Dans l'autre, la moitié droite est la moitié gauche x -1. Cela signifie que nous n'avons vraiment besoin d'utiliser la moitié gauche de ces matrices pour la multiplication et de créer la moitié droite à moindre coût en multipliant par 1 ou -1. Ensuite, observez que la deuxième matrice diffère de la première matrice par des facteurs qui sont les mêmes dans chaque colonne, donc nous pouvons factoriser cela et le multiplier dans l'entrée alors maintenant les échantillons pairs et impairs utilisent la même matrice, mais nécessitent un multiplicateur premier. Et la dernière étape consiste à observer que cette matrice N / 2 × N / 2 résultante est identique à une matrice N / 2 DFT et nous pouvons le faire encore et encore jusqu'à atteindre une matrice 1 × 1 où la DFT est une fonction d'identité.

Pour généraliser au-delà de Radix 2, vous pouvez envisager de fractionner chaque troisième ligne et de regarder trois morceaux de colonnes, ou tous les 4, etc.

Dans le cas d'entrées de taille principale, il existe une méthode pour correctement mettre à zéro, FFT et tronquer, mais cela dépasse le cadre de cette réponse.

Voir: http://whoiskylefinn.com/MatrixFFT.html

kylefinn
la source
FFT prime , diverses FFT . L'utilisation du zéro-pad n'est pas la seule option. Désolé, je trouve juste le remplissage nul surutilisé. Une petite question, je ne comprends pas ce que vous entendez par "chaque couche avec une exigence de temps constant par échantillon", si vous pouviez expliquer, ce serait génial.
Evil
1
Désolé, je ne voulais pas dire que le remplissage zéro était LA voie, je voulais juste indiquer une lecture plus approfondie. Et "couche" signifiant une récursivité, ou une traduction d'un N DFT à 2 N / 2 DFT, avec un temps constant par échantillon signifiant que cette étape est O (N).
kylefinn
Jusqu'à présent, de toutes les descriptions, celle-ci semble la plus proche de simplifier un problème complexe. La grande chose qui manque, cependant, est un exemple de ces matrices. En auriez-vous un?
Seph Reed
Téléchargé cela, devrait aider: whoiskylefinn.com/MatrixFFT.html
kylefinn
1

La DFT multiplie une matrice N ^ 2 de force brute.

Les FFT font des astuces intelligentes, exploitant les propriétés de la matrice (dégénéralisant la multiplication de la matrice) afin de réduire les coûts de calcul.

Voyons d'abord un petit DFT:

W = fft (œil (4));

x = rand (4,1) + 1j * rand (4,1);

X_ref = fft (x);

X = W * x;

affirmer (max (abs (X-X_ref)) <1e-7)

Très bien, nous sommes donc en mesure de remplacer l'appel MATLABs à la bibliothèque FFTW par une petite multiplication matricielle 4x4 (complexe) en remplissant une matrice à partir de la fonction FFT. Alors à quoi ressemble cette matrice?

N = 4,

Wn = exp (-1j * 2 * pi / N),

f = ((0: N-1) '* (0: N-1))

f =

 0     0     0     0
 0     1     2     3
 0     2     4     6
 0     3     6     9

W = Wn. ^ F

W =

1 1 1 1

1 -i -1 i

1 -1 1 -1

1 i -1 -i

Chaque élément est soit +1, -1, + 1j ou -1j. Évidemment, cela signifie que nous pouvons éviter les multiplications complexes et complètes. De plus, la première colonne est identique, ce qui signifie que nous multiplions le premier élément de x encore et encore par le même facteur.

Il s'avère que les produits tensoriels de Kronecker, les «facteurs de torsion» et une matrice de permutation où l'indice est modifié en fonction de la représentation binaire inversée sont à la fois compacts et donnent une perspective alternative sur la façon dont les FFT sont calculées comme un ensemble d'opérations matricielles clairsemées.

Les lignes ci-dessous sont une simple FFT vers l'avant de décimation en fréquence (DIF). Bien que les étapes puissent sembler lourdes, il est pratique de les réutiliser pour la FFT directe / inverse, radix4 / split-radix ou la décimation dans le temps, tout en étant une représentation fidèle de la façon dont les FFT en place ont tendance à être mises en œuvre dans le monde réel, Je crois.

N = 4;

x = randn (N, 1) + 1j * randn (N, 1);

T1 = exp (-1j * 2 * pi * ([zéros (1, N / 2), 0: (N / 2-1)]). '/ N),

M0 = kron (œil (2), fft (œil (2))),

M1 = kron (fft (œil (2)), œil (2)),

X = bitorder (x. '* M1 * diag (T1) * M0),

X_ref = fft (x)

assert (max (abs (X (:) - X_ref (:))) <1e-6)

CF Van Loan a un excellent livre sur ce sujet.

Knut Inge
la source