Präzise OpenCL-FFT

Richard Gast #2459772
Lesenswert?
• ▲
▼ 
Ich arbeite gerade an einer blockweisen Korrelationssuche im 
Frequenzbereich.
Das Verfahren funktioniert wunderbar mit einer "klassischen" CPU-FFT, 
der Cuda-FFT - aber versagt jämmerlich bei der OpenCL-FFT, wenn die 
FFT-Größe 1024 überschreitet (ich verwende die in JavaCL RC1 angebotene 
FFT).

Die Ursache liegt in der Ungenauigkeit der FFT-Berechnung. Genauer 
gesagt liefert
 iFFT(FFT(x)*conj(FFT(y)))
mit x und y als reelle Vektoren keinen reellen Ergebnisvektor.

Deshalb die Frage: gibt es irgendwo einen brauchbaren FFT-Kernel für 
Float-Felder, der möglichst auch die reel->komplex Transformation 
beherrscht - und trotzdem in der FFT-Länge variabel ist (Zweierpotenzen 
vorausgesetzt).
Antwort Zitieren Bearbeiten Thread verschieben Thread sperren Anmeldepflicht aktivieren Anpinnen Thread löschen Thread mit anderem zusammenführen Beitrag melden
HalfBit Gast #2461354
Lesenswert?
• ▲
▼ 
Das hängt mit der internen Genauigkeit der FFT-Implementierung zusammen.

Und da muss man leider länger ausholen: CUDA unterstützt den Datentyp 
double nativ (ab Fermi-GPU). Bei OpenCL wird bis heute der Typ double 
nur optional unterstützt, die Pragma ist
1
#pragma OPENCL EXTENSION cl_khr_fp64 : enable


Wie OpenCL damit umgeht, bleibt aber den Anbeiter überlassen. Sprich: 
Selbst wenn der JavaCL-Wrapper eine double-FFT anbietet, muss OpenCL 
keineswegs mit double arbeiten und konvertiert vor der FFT zu float.

Richard schrieb:
> Deshalb die Frage: gibt es irgendwo einen brauchbaren FFT-Kernel für
> Float-Felder, der möglichst auch die reel->komplex Transformation
> beherrscht - und trotzdem in der FFT-Länge variabel ist (Zweierpotenzen
> vorausgesetzt).

Ja, viele. Der Fachausdruck heißt "Split-Radix-FFT", die Länge des 
Vektors ist dann frei wählbar und ist noch nicht mal an 2er-Potenzen 
gebunden. Etwa in FFTW3, Parallel-FFTW, CuFFT, MathNet Numerics uvam. zu 
finden. Leider weiß ich nicht, ob es hierfür Java-Wrapper gibt.
Antwort Zitieren Bearbeiten Löschen Beitrag melden
Andreas S. Admin Persönliche Seite #2461689
Lesenswert?
• ▲
▼ 
Richard schrieb:
> Die Ursache liegt in der Ungenauigkeit der FFT-Berechnung. Genauer
> gesagt liefert
>  iFFT(FFT(x)*conj(FFT(y)))
> mit x und y als reelle Vektoren keinen reellen Ergebnisvektor.

Darauf darf man sich auch niemals verlassen, weil es bei beschränkter 
Rechengenauigkeit nicht garantiert werden kann.

Ist dein Problem wirklich die Genauigkeit, oder dass dein Verfahren (was 
auch immer es ist) nicht damit klarkommt wenn ein Wert != 0 im 
Imaginärteil auftaucht?
Antwort Zitieren Bearbeiten Löschen Beitrag melden
Max G. Persönliche Seite #2466622
Lesenswert?
• ▲
▼ 
Richard schrieb:

> Die Ursache liegt in der Ungenauigkeit der FFT-Berechnung. Genauer
> gesagt liefert
>  iFFT(FFT(x)*conj(FFT(y)))
> mit x und y als reelle Vektoren keinen reellen Ergebnisvektor.

Das ist aber nichts Ungewöhnliches, bei Matlab hat mich das Gleiche auch 
immer wieder geärgert. Bei näherer Betrachtung sind die imaginären 
Anteile dann irgendwo im Bereich von 10-15 Größenordnungen kleiner als 
die realen Anteile, so dass man sie getrost ignorieren kann.

Wenn du weißt, dass das Ergebnis reell sein muss, warum wirfst du dann 
den Imaginärteil nicht einfach weg (Korrektheit des Algorithmus 
natürlich vorausgesetzt)?

Max
Antwort Zitieren Bearbeiten Löschen Beitrag melden
Thread beobachten

Antwort schreiben

Bitte melde dich an, um einen Beitrag zu schreiben.

Bestehender Account

oder

Mit Google-Account einloggen

Noch kein Account?

Die Registrierung ist kostenlos und dauert nur eine Minute.

Jetzt registrieren