Präzise OpenCL-FFT

[Richard (Gast)]

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).

[HalfBit (Gast)]

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

#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.

[Andreas S. (andreas) (Admin)]

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?

[Max G. (l0wside)]

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