2D FFT

101 views
Skip to first unread message

simon...@gmail.com

unread,
Mar 10, 2017, 11:43:30 AM3/10/17
to reikna
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

unread,
Mar 13, 2017, 7:35:35 PM3/13/17
to reikna
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?
Message has been deleted
Message has been deleted
Message has been deleted

simon...@gmail.com

unread,
Mar 14, 2017, 7:22:19 AM3/14/17
to reikna
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

unread,
Mar 16, 2017, 3:08:02 AM3/16/17
to reikna
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

unread,
Mar 16, 2017, 3:14:08 PM3/16/17
to reikna
Thanks sorted now.
Reply all
Reply to author
Forward
0 new messages