2D FFT

[simon...@gmail.com]

Hi

I am having a little trouble using the FFT function. I wonder if you can show me where I am bring stupid. I am trying to FFT a 1024 by 1024 numpy array both by numpy and reikna but they seem to give a different result (H and result). Is this to be expected?

Here is a snippet of the code.

import numpy as np
import reikna.cluda as cluda
from reikna.cluda import dtypes, functions
import reikna.fft as cudafft

api=cluda.cuda_api()
thr=api.Thread.create()

I1 = image
I1_dev=thr.to_device(I1)
outcomplex=thr.array((1024, 1024), dtype=np.complex128)
outreal = thr.array((1024, 1024), dtype=np.float64)

ffts = cudafft.FFTShift(I1_dev)
fftsc = ffts.compile(thr)

fft=cudafft.FFT(out)
fftc=fft.compile(thr)

for n in range(0,1): # Iterations to optimize the phase hologram

fftsc(outreal,I1_dev)
test=outreal.get()
test=test.astype(np.complex128)
test_dev=thr.to_device(test)
fftc(outcomplex,test_dev)
result = outcomplex

H = ((np.fft.fft2(np.fft.fftshift(I1))))

[Bogdan Opanchuk]

No, the result should be the same. Your code is not actually executable, so I cannot verify that you are checking the same thing as I do, but with minimal changes I get this:

import numpy as np
import reikna.cluda as cluda
from reikna.cluda import dtypes, functions
import reikna.fft as cudafft


api=cluda.cuda_api()
thr=api.Thread.create()

I1 = np.random.normal(size=(1024, 1024))
I1_dev = thr.to_device(I1)
outcomplex = thr.array((1024, 1024), dtype=np.complex128)

outreal = thr.array((1024, 1024), dtype=np.float64)

ffts = cudafft.FFTShift(I1_dev)
fftsc = ffts.compile(thr)

fft = cudafft.FFT(outcomplex)
fftc = fft.compile(thr)

fftsc(outreal,I1_dev)
test=outreal.get()
test=test.astype(np.complex128)
test_dev=thr.to_device(test)
fftc(outcomplex,test_dev)

result = outcomplex.get()


H = ((np.fft.fft2(np.fft.fftshift(I1))))

print(np.linalg.norm(H - result) / np.linalg.norm(H))

This gives the relative difference of the order of 10^(-15) for me. Could you run it and see what you get?

[simon...@gmail.com]

Hi thanks. Tried that code and it did work nicely. However when I compared the results using the np.angle() command I got very different results. Sorry for my previous code being a snippet, hopfully this should run for you. The difference between the j and g arrays should be visible on plotting. and the max(p) is significant. Is this to be expected?

import numpy as np
import reikna.cluda as cluda
from reikna.cluda import dtypes, functions
import reikna.fft as cudafft

api=cluda.cuda_api()
thr=api.Thread.create()

size=1024
I=np.zeros([size,size])
#set up a square for the sake of imaging.
x=np.arange(-size/2,size/2,1) #make a range that is X and Y size
y=np.arange(-size/2,size/2,1)
xv, yv = np.meshgrid(x, y, sparse=True)
xv=np.transpose(yv)
yv=np.transpose(xv)
test=(abs(xv)<50/2) & (abs(yv)<50/2)

I[test]=1
I1_dev = thr.to_device(I)

outcomplex = thr.array((1024, 1024), dtype=np.complex128)
outreal = thr.array((1024, 1024), dtype=np.float64)

ffts = cudafft.FFTShift(I1_dev)
fftsc = ffts.compile(thr)

fft2s = cudafft.FFTShift(outcomplex)
fft2sc = fft2s.compile(thr)

fft = cudafft.FFT(outcomplex)
fftc = fft.compile(thr)

fftsc(outreal,I1_dev)
test=outreal.get()
test=test.astype(np.complex128)
test_dev=thr.to_device(test)
fftc(outcomplex,test_dev)

fft2sc(outcomplex,outcomplex)
result = outcomplex.get()
H = (np.fft.fftshift(np.fft.fft2(np.fft.fftshift(I))))

print(np.linalg.norm(H - result) / np.linalg.norm(H))

j=np.angle(result)
g=np.angle(H)
p=j-g
print(np.max(p))

[Bogdan Opanchuk]

I think I know what is happening here. There are some gotchas in comparing angles. Your maximum difference ("np.max(p)") is almost 2pi, which is the same as 0. If you get rid of this circularity and compare, e.g.

print(np.allclose(np.exp(1j*j), np.exp(1j*g)))

you will get "True" (at least I do when I run your code). That is, everything works as expected.

[simon...@gmail.com]

Thanks sorted now.