2017-12-19 10 views
1

私はpyculibを使用して、Anaconda 3.5のマトリックスで3D FFTを実行します。私はちょうどウェブサイトに掲載されたthe example codeをたどった。しかし、私は何か興味深いものを見つけ出し、理由を理解していません。pyculibを使用した3D FFTの結果が間違っています

pyculibを使用して行列に3D FFTを実行するのは、numpy.arangeを使用して行列を作成した場合のみです。ここで

コードです:

from pyculib.fft.binding import Plan, CUFFT_C2C 
import numpy as np 
from numba import cuda 
data = np.random.rand(26, 256, 256).astype(np.complex64) 
orig = data.copy() 
d_data = cuda.to_device(data) 
fftplan = Plan.three(CUFFT_C2C, *data.shape) 
fftplan.forward(d_data, d_data) 
fftplan.inverse(d_data, d_data) 
d_data.copy_to_host(data) 
result = data/n 
np.allclose(orig, result.real) 

最後に、それは、であることが判明しました。 origおよび結果の違いはごくわずかではありません。 他のデータセット(乱数ではない)を試してみましたが、間違った結果が出ました。

また、私は逆FFTなしテスト:

from pyculib.fft.binding import Plan, CUFFT_C2C 
import numpy as np 
from numba import cuda 
from scipy.fftpack import fftn,ifftn 
data = np.random.rand(26,256,256).astype(np.complex64) 
orig = data.copy() 
orig_fft = fftn(orig) 
d_data = cuda.to_device(data) 
fftplan = Plan.three(CUFFT_C2C, *data.shape) 
fftplan.forward(d_data, d_data) 
d_data.copy_to_host(data) 
np.allclose(orig_fft, data) 

結果も間違っています。

ウェブサイトのテストコードは、numpy.arangeを使用してマトリックスを作成します。私は試しました:

n = 26*256*256 
data = np.arange(n, dtype=np.complex64).reshape(26,256,256) 

この行列のFFT結果は正しいです。

誰でも原因を指摘できますか?

答えて

0

私はCUDAを使用していませんが、あなたの問題は本質的に数字だと思います。違いは、使用している2つのデータセットにあります。 random.randのダイナミックレンジは0〜1、範囲は0〜26 * 256 * 256です。 FFTは、値の範囲/ポイント数のオーダーで空間周波数成分を分解しようと試みる。ある範囲については、これは1になり、FFTは数値的に正確です。 randの場合、これは1/26 * 256 * 256〜5.8e-7です。

CUDAを使用せずにnumpy配列でFFT/IFFTを実行するだけで、同様の違いが表示されます。

+0

私は同意する、それはいくつかの数値的な問題かもしれない。興味深いのは、今日、同じことをやっているFFTW3を試したことです。そして、pyfftwとnumpy.fftはnumpy.rand入力と同じであることが分かります。私はカフに何か問題があったのだろうかと思っています。 – billinair

関連する問題