Ich versuche FFTs zu verstehen, hier ist was ich bisher habe:
Um die Größe von Frequenzen in einer Wellenform zu finden, muss man nach ihnen suchen, indem man die Welle mit der Frequenz, nach der sie suchen, in zwei verschiedenen Phasen (sin und cos) multipliziert und jeweils einen Durchschnitt bildet. Die Phase wird durch ihre Beziehung zu den beiden gefunden, und der Code dafür ist ungefähr so:
//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;
}
Um dies für sehr viele Frequenzen sehr schnell zu tun, verwenden FFTs viele Tricks.
Mit welchen Tricks sind FFTs so viel schneller als DFT?
PS Ich habe versucht, fertige FFT-Algorithmen im Web zu betrachten, aber alle Tricks lassen sich ohne viel Erklärung zu einem schönen Stück Code zusammenfassen. Was ich zuerst brauche, um das Ganze zu verstehen, ist eine Einführung in jede dieser effizienten Änderungen als Konzepte.
Vielen Dank.
quelle
sudo
in Ihrem Codebeispiel verwirrend sein kann, da dies ein in der Computerwelt bekannter Befehl ist. Du hast wahrscheinlich Pseudocode gemeint.Antworten:
Die naive Implementierung einer Punkt-DFT ist im Grunde genommen eine Multiplikation mit einer N × N- Matrix. Dies führt zu einer Komplexität von O ( N 2 ) .N N×N O(N2)
Einer der gebräuchlichsten Fast Fourier Transform (FFT) -Algorithmen ist der Cooley-Tukey-Decimation-in-Time-FFT-Algorithmus ( Radix-2) . Dies ist ein grundlegender Ansatz zum Teilen und Erobern.
Definieren Sie zuerst den "Twiddle-Faktor" als: woj≜√
quelle
W
,j
,X()
,N
undk
noch nicht haben Definitionen für mich.http://nbviewer.jupyter.org/gist/leftaroundabout/83df89a7d3bdc24373ea470fb50be629
DFT, Größe 16
FFT, Größe 16
Der Unterschied in der Komplexität ist ziemlich offensichtlich, nicht wahr?
So verstehe ich FFT.
Die gemessenen Daten müssen jedoch nicht unbedingt einer physikalischen Grundgröße entsprechen. Wenn Sie beispielsweise eine Lichtintensität messen, messen Sie tatsächlich nur die Amplitude einer elektromagnetischen Welle, deren Frequenz selbst zu hoch ist, um mit einem ADC abgetastet zu werden. Es ist jedoch klar, dass Sie auch die DFT eines abgetasteten Lichtintensitätssignals berechnen können, und das zu einem günstigen Preis, trotz der irrsinnigen Frequenz der Lichtwelle.
Dies kann als der wichtigste Grund für die Kostengünstigkeit von FFT angesehen werden:
Versuchen Sie nicht, die einzelnen Schwingungszyklen von der höchsten Ebene aus zu sehen. Transformieren Sie stattdessen nur Informationen auf höherer Ebene, die bereits lokal vorverarbeitet wurden.
Das ist jedoch noch nicht alles. Das Tolle an FFT ist, dass Sie immer noch alle Informationen erhalten, die eine vollständige DFT liefern würde . Dh alle Informationen, die Sie auch erhalten, wenn Sie die exakte elektromagnetische Welle eines Lichtstrahls abtasten. Kann dies durch Transformieren eines Fotodiodensignals erreicht werden? - Können Sie daraus die genaue Lichtfrequenz messen?
Durch eine insgesamt längere Zeitspanne sollten wir auch in der Lage sein, die Frequenzunsicherheit einzugrenzen. Und dies ist tatsächlich möglich, wenn Sie nicht nur die Grobfrequenz, sondern auch die Phase der Welle lokal messen . Sie wissen, dass ein 1000-Hz-Signal genau dieselbe Phase hat, wenn Sie es eine Sekunde später betrachten. Während ein 1000,5-Hz-Signal auf der kurzen Skala nicht zu unterscheiden ist, hat es eine Sekunde später eine invertierte Phase.
Glücklicherweise kann diese Phaseninformation sehr gut in einer einzigen komplexen Zahl gespeichert werden. Und so funktioniert FFT! Es beginnt mit vielen kleinen, lokalen Veränderungen. Diese sind billig - zum einen offensichtlich, weil sie nur eine geringe Datenmenge verwenden, zum anderen, weil sie wissen, dass sie die Frequenz aufgrund der kurzen Zeitspanne ohnehin nicht sehr genau auflösen können - und daher auch für Sie noch erschwinglich sind mache eine ganze Menge solcher Transformationen.
Diese zeichnen jedoch auch die Phase auf , woraufhin Sie die Frequenzauflösung auf der obersten Ebene genauer einstellen können. Die erforderliche Transformation ist wieder günstig, da sie sich nicht um hochfrequente Schwingungen kümmert, sondern nur um die vorverarbeiteten niederfrequenten Daten.
† Ja, meine Argumentation ist an dieser Stelle etwas kreisförmig. Nennen wir es einfach rekursiv und es geht uns gut ...
‡ Diese Beziehung ist nicht quantenmechanisch, aber die Heisenbergsche Unsicherheit hat tatsächlich den gleichen fundamentalen Grund.
quelle
Beachten Sie den gezeigten Pfad und die Gleichung darunter zeigt das Ergebnis für das Frequenzfach X (1), wie es durch Roberts Gleichung gegeben ist.
Gestrichelte Linien unterscheiden sich nicht von durchgezogenen Linien, um zu verdeutlichen, wo sich die Summationsverknüpfungen befinden.
quelle
Bei der Berechnung der naiven DFT direkt aus der Summe gilt:
quelle
Ich bin eine visuelle Person. Ich stelle mir die FFT lieber als Matrixtrick als als Summationstrick vor.
Auf hohem Niveau erklären:
Eine naive DFT berechnet jede Ausgangsabtastung unabhängig und verwendet jede Eingangsabtastung in jeder Berechnung (klassischer N²-Algorithmus).
Eine übliche FFT verwendet Symmetrien und Muster in der DFT-Definition, um die Berechnung in "Schichten" (log N Schichten) durchzuführen, wobei jede Schicht mit konstantem Zeitbedarf pro Stichprobe einen N log N-Algorithmus erzeugt.
Weitere Einzelheiten:
Eine Möglichkeit, diese Symmetrien zu visualisieren, besteht darin, die DFT als 1 × N-Matrixeingang multipliziert mit einer N × N-Matrix zu betrachten aller Ihrer komplexen Exponentialfunktionen. Beginnen wir mit dem Fall "radix 2". Wir werden die geraden und ungeraden Zeilen der Matrix aufteilen (entsprechend den geraden und ungeraden Eingangsabtastwerten) und sie als zwei separate Matrixmultiplikationen betrachten, die sich addieren, um das gleiche Endergebnis zu erhalten.
Schauen Sie sich nun diese Matrizen an: In der ersten Hälfte ist die linke mit der rechten Hälfte identisch. In der anderen ist die rechte Hälfte die linke Hälfte x −1. Das heißt, wir müssen wirklich nur die linke Hälfte dieser Matrizen für die Multiplikation verwenden und die rechte Hälfte kostengünstig durch Multiplikation mit 1 oder –1 erstellen. Beachten Sie als Nächstes, dass sich die zweite Matrix von der ersten Matrix durch in jeder Spalte identische Faktoren unterscheidet, sodass wir diese herausrechnen und in die Eingabe multiplizieren können, sodass jetzt sowohl gerade als auch ungerade Samples dieselbe Matrix verwenden, jedoch einen Multiplikator erfordern zuerst. Und der letzte Schritt besteht darin, zu beobachten, dass diese resultierende N / 2 × N / 2-Matrix mit einer N / 2-DFT-Matrix identisch ist, und wir können dies immer wieder tun, bis wir eine 1 × 1-Matrix erreichen, bei der die DFT eine Identitätsfunktion ist.
Um über die Basis 2 hinaus zu verallgemeinern, können Sie das Teilen jeder dritten Zeile und das Betrachten von drei Spaltenblöcken oder jedem vierten usw. betrachten.
Für den Fall von Eingaben mit Primzahlgröße gibt es eine Methode zum korrekten Nullstellen, FFT und Abschneiden, die jedoch den Rahmen dieser Antwort sprengt.
Siehe: http://whoiskylefinn.com/MatrixFFT.html
quelle
Die DFT multipliziert eine Brute-Force-N ^ 2-Matrix.
FFTs machen clevere Tricks, indem sie die Eigenschaften der Matrix ausnutzen (degeneralisieren die Matrix multiplizieren), um die Rechenkosten zu senken.
Schauen wir uns zunächst eine kleine DFT an:
W = fft (Auge (4));
x = rand (4,1) + 1j * rand (4,1);
X_ref = fft (x);
X = W * x;
assert (max (abs (X-X_ref)) <1e-7)
Gut, dass wir MATLABs, die die FFTW-Bibliothek aufrufen, durch eine kleine (komplexe) 4x4-Matrixmultiplikation ersetzen können, indem wir eine Matrix aus der FFT-Funktion füllen. Wie sieht diese Matrix aus?
N = 4,
Wn = exp (-1j · 2 · pi / N),
f = ((0: N-1) '* (0: N-1))
f =
W = Wn. ^ F
W =
1 1 1 1
1 -i -1 i
1 -1 1 -1
1 i -1 -i
Jedes Element ist entweder +1, -1, + 1j oder -1j. Dies bedeutet natürlich, dass wir komplexe Multiplikationen vermeiden können. Ferner ist die erste Spalte identisch, was bedeutet, dass wir das erste Element von x immer wieder mit demselben Faktor multiplizieren.
Es stellt sich heraus, dass Kronecker-Tensorprodukte, "Twiddle-Faktoren" und eine Permutationsmatrix, in der der Index entsprechend der gespiegelten binären Repräsentation geändert wird, kompakt sind und eine alternative Perspektive für die Berechnung von FFTs als Satz von Operationen mit geringer Matrixdichte bieten.
Die folgenden Zeilen zeigen eine einfache DIF-Vorwärts-FFT (Decimation in Frequency) mit der Basis 2. Die Schritte mögen mühsam erscheinen, es ist jedoch zweckmäßig, sie für Forward / Inverse FFT, Radix4 / Split-Radix oder Decimation-in-Time wiederzuverwenden, während die Implementierung von In-Place-FFTs in der realen Welt angemessen dargestellt wird. Ich glaube.
N = 4;
x = Randn (N, 1) + 1j * Randn (N, 1);
T1 = exp (-1j · 2 · pi · ([Nullen (1, N / 2), 0: (N / 2-1)]). '/ N),
M0 = kron (Auge (2), fft (Auge (2))
M1 = Kronen (fft (Auge (2)), Auge (2)),
X = Bitreihenfolge (x. '* M1 * Diag (T1) * M0),
X_ref = fft (x)
assert (max (abs (X (:) - X_ref (:))) <1e-6)
CF Van Loan hat ein großartiges Buch zu diesem Thema.
quelle
Wenn du aus dem Feuer der Weisheit trinken willst, schlage ich vor:
"Schnelle Transformationen - Algorithmen, Analysen, Anwendungen" von Douglas F. Elliott, K. Ramamohan Rao
Es deckt FFT, Hartley, Winograd und Anwendungen ab.
Ein starker Punkt ist, dass es sich bei der FFT um eine Reihe von Faktorisierungen mit geringer Matrixdichte mit Bitumkehrreihenfolge handelt.
quelle