115 views

Skip to first unread message

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

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

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

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

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.

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.

Mar 16, 2017, 3:14:08 PM3/16/17

to reikna

Thanks sorted now.

Reply all

Reply to author

Forward

0 new messages

Search

Clear search

Close search

Google apps

Main menu