Schnelle "Schnelle Fouriertransformation" gesucht
-
Hi.
Ich arbeite gerade an einem Programm, in dem oft eine 3D-Fouriertransformation gemacht wird. Gitterlänge des Gitters, auf dem die Fouriertransformation gemacht werden muss, kann ich zwar kontrollieren, die "natürliche Länge" ist aber oft ungünstig. Und zwar arbeite ich auf Gittern die mindestens die Größe (4n+1)(4n+1)(4n+1) haben müssen, wobei n eine kleine ganze Zahl zwischen 1 und 20 ist. Wie Ihr seht, kommt da als kleinstmögliches Gitter nicht gerade ein Gitter mit Seitenlänge 2^N heraus, wie es für schnelle Fouriertransformationen eigentlich typischerweise genutzt wird. Stattdessen kriege ich Seitenlängen von 37 (n=9), 41 (n=10), 53 (n=13) und so. Das heißt, ich habe oft Seitenlängen, die Primzahlen sind. Das ist natürlich ganz schlecht, denn AFAIK skalieren alle Abarten von FFTs für beliebige Seitenlängen bei Primzahlen nur mit O(n^2), bei n=2^N-Gittern hingegen mit O(n*log(n)). Für alle anderen Gitter skalieren sie mit irgendetwas dazwischen.
Die Frage ist, wie ich die Fouriertransformationen am Schnellsten hinkriege. Ich nutze momentan den "mixed-radix fast fourier transform algorithm" von R. Singleton. Meine Idee war es zuerst, das Gitter etwas zu vergrößern, damit der Algorithmus mit der Seitenlänge besser zurecht kommt. Zum Beispiel habe ich versucht, ein 41x41x41 Gitter durch ein 64x64x64 Gitter zu ersetzen. Mit dem Erfolg, dass es langsamer geworden ist. Von der Theorie her hätte es schneller werden müssen, aber in dem Fall kommen glaube ich auch noch so Effekte wie die Cache-Größe ins Spiel. Letztendlich hat die 41-Variante etwa 15 Sekunden gebraucht, die 64-Variante hingegen 70 Sekunden. Bei einem 42er-Gitter bin ich mit 10 Sekunden ausgekommen. So, ich muss aber mit 2-3 Sekunden auskommen oder wenigstens so schnell wie möglich sein. Dafür suche ich jetzt eine Lösung.
Die erste Frage, die mir diesbezüglich einfällt ist, wie konkurrenzfähig dieser FFT-Algorithmus von R. Singleton ist. Kennt sich jemand damit aus und kann sagen, wie der sich im Vergleich zu anderen FFT-Algorithmen verhält?
Und: Hat jemand eine Idee, wie man eine "Optimale Gittergröße" bestimmen kann? (Hierbei gehe ich mal davon aus, dass die Frage zu speziell ist.) Ich könnte mir vorstellen, dass man sich ein paar Gitter ansieht, die "etwas" größer als das Minimalgitter sind und dann guckt, wie die Primfaktoren der Gitterlängen aussehen. Ich gehe davon aus, dass die FFT schneller ist, wenn die Primfaktoren kleiner sind. Gibt es dafür eine Regel, mit der man das quantitativ abschätzen kann? Jenseits von Cache-Effekten natürlich.
Wär echt toll, wenn einer von Euch Ahnung auf dem Gebiet hat und mir ein paar wertvolle Hinweise geben kann. Wie sieht es mit der FFTW aus? Ist die wesentlich schneller als dieser Algorithmus von R. Singleton? Und kann die mit Gittern beliebiger Größe umgehen?
-
Gregor schrieb:
Das ist natürlich ganz schlecht, denn AFAIK skalieren alle Abarten von FFTs für beliebige Seitenlängen bei Primzahlen nur mit O(n^2),
Es gibt angeblich auch spezielle Algorithmen für Primzahllängen. Wie die funktionieren, weiss ich aber nicht.
Gregor schrieb:
Zum Beispiel habe ich versucht, ein 41x41x41 Gitter durch ein 64x64x64 Gitter zu ersetzen. Mit dem Erfolg, dass es langsamer geworden ist.
Also Zero-Padding ist erlaubt, ja?
Gregor schrieb:
Die erste Frage, die mir diesbezüglich einfällt ist, wie konkurrenzfähig dieser FFT-Algorithmus von R. Singleton ist. Kennt sich jemand damit aus und kann sagen, wie der sich im Vergleich zu anderen FFT-Algorithmen verhält?
Keine Ahnung. Probier's aus.
Gregor schrieb:
... "Optimale Gittergröße" ... Primfaktoren ... Wär echt toll, wenn einer von Euch Ahnung auf dem Gebiet hat und mir ein paar wertvolle Hinweise geben kann. Wie sieht es mit der FFTW aus? Ist die wesentlich schneller als dieser Algorithmus von R. Singleton? Und kann die mit Gittern beliebiger Größe umgehen?
Also... Algorithmen sind das eine, die Implementierungen das andere. Kann sein, dass Dein Algorithmus theoreitsch ganz toll ist, aber er schlecht implementiert wude. Die FFTW-Bibliothek unterstützt beliebige Größen, vorzugsweise mit vielen kleinen Primfaktoren wie 2,3 und 5. Bei der FFTW läuft das ganze so ab, dass man vorher einen "Plan" erstellen muss von dem, was man machen möchte. Beim Erzeugen des Plans werden typischerweise versch. Tests durchgeführt, mit denen dann bestimmt wird, welche Faktorisierung der DFT für die angegebene Dimension und gegebenen Speicherlayout am effektivsten funktioniert. Ist der Plan erstellt, kann man ihn relativ leicht (und oft) ausführen lassen. Ich würde auch mal stark davon ausgehen, dass die FFTW auch ordentlich SSE-optimiert ist. Zumindest soll man fftw_malloc und fftw_free statt malloc/free oder new[]/delete[] benutzen, da fftw_malloc garantiert einen Zeiger zurück gibt, der "SSE-aligned" ist.
Wenn Du Zero-Padding erlauben kannst, denke ich, dass es sich lohnt zB die 41 auf 48 zu "runden". 48 ist als Mittelwert zwischen zwei benachbarten Zweierpotenzen interessant, da sie auch aus nur kleinen Primfaktoren besteht, nämlich nur Zweien und einer Drei. Damit sollte die FFTW wunderbar klarkommen. Vielleicht verkraftet sie sogar Primzahlen. Guck doch mal in der Doku nach.
Ich würde Dir die FFTW empfehlen. Es macht nicht viel Sinn, das Rad neuzuerfinden, wenn es hier auf Geschwindigkeit ankommt.
-
krümelkacker schrieb:
Wenn Du Zero-Padding erlauben kannst, denke ich, dass es sich lohnt zB die 41 auf 48 zu "runden". 48 ist als Mittelwert zwischen zwei benachbarten Zweierpotenzen interessant, da sie auch aus nur kleinen Primfaktoren besteht, nämlich nur Zweien und eine Drei. Damit sollte die FFTW wunderbar klarkommen. Vielleicht verkraftet sie sogar Primzahlen. Guck doch mal in der Doku nach.
Das habe ich gemacht. Bei meinem jetzigen Algorithmus schlägt dann schon der Cache-Effekt ein. Mit einer 48er-Kantenlänge brauche ich so um die 20 Sekunden. Vermutlich werde ich die FFTW mal ausprobieren.
Zero-Padding ist insofern erlaubt, als ich Nullen in die Mitte einfügen darf. Ich habe in der Mitte des FFT-Gitters (im Frequenzraum) also ein Kreuz aus Nullen.
Und mich interessiert auch die Transformation von Frequenzraum in den Realraum.
-
Gregor schrieb:
Bei meinem jetzigen Algorithmus schlägt dann schon der Cache-Effekt ein. Mit einer 48er-Kantenlänge brauche ich so um die 20 Sekunden.
Klingt nach sehr viel -- zu viel. Gut, ob das mit kleinen Primfaktoren der 48 gegenüber 41 viel bringt (bei solch kleinen Zahlen), weiß ich nicht. Aber Du transformierst ja in mehreren Dimensionen. Beispielsweise könnte man bei einer 2D-Transformation, wo die Daten Zeilenweise im Speicher liegen (row major storage) mehrere Spalten gleichzeitig per SSE transformieren, um es (a) Cache-freundlicher zu machen und (b) mehrere Multiplikationen auf einen Schlag zu erledigen. Solche Tricks sind vermutlich in der FFTW eingebaut.
Gregor schrieb:
Vermutlich werde ich die FFTW mal ausprobieren.
Mach das.
-
krümelkacker schrieb:
Gregor schrieb:
Bei meinem jetzigen Algorithmus schlägt dann schon der Cache-Effekt ein. Mit einer 48er-Kantenlänge brauche ich so um die 20 Sekunden.
Klingt nach sehr viel -- zu viel.
Ich haette vielleicht dazu sagen sollen, dass es sich nicht um _eine_ FFT handelt, sondern um 2000.
Ausserdem scheint der Cache-Effekt mehr oder weniger meine Einbildung gewesen zu sein.krümelkacker schrieb:
Gregor schrieb:
Vermutlich werde ich die FFTW mal ausprobieren.
Mach das.
So, habe FFTW jetzt mal getestet. Sieht ganz gut aus, wenn auch nicht ganz so gut, wie zunaechst erhofft. Wen das Verhalten von FFTW bzw. von schnellen Fouriertransformationen im Allgemeinen in Abhaengigkeit von der Groesse des Gitters interessiert, der kann sich mal folgenden Graph angucken, der eben bei mir entstanden ist:
http://img688.imageshack.us/img688/5927/ffttimes.png
Sorry fuer die ganzen unnoetigen Informationen auf dem Graph
Bezueglich FFTW ist eigentlich nur die rote Kurve wichtg (die schwarze ist die bisher genutzte FFT). Die Ausreisser nach oben liegen bei den Primzahlen. Entsprechend habe ich 2 unterschiedliche Kurven bezueglich der Skalierung eingebaut. Da es eine 3D-Fouriertransformation ist, skalieren Gitter mit Laenge N=Primzahl gemaess O(N^4). Wenn man eine Laenge der Art N=2^n hat, dann skaliert die 3D-FFT wie O(N^3*ld(N)).
