def test_against_fft_3d(self):
from pyculib.fft.binding import Plan, CUFFT_R2C
depth = 2
colsize = 2
rowsize = 64
N = depth * colsize * rowsize
x = np.arange(N, dtype=np.float32).reshape(depth, colsize, rowsize)
xf = np.fft.fftn(x)
d_x_gpu = cuda.to_device(x)
xf_gpu = np.zeros(shape=(depth, colsize, rowsize//2 + 1), dtype=np.complex64)
d_xf_gpu = cuda.to_device(xf_gpu)
plan = Plan.many(x.shape, CUFFT_R2C)
plan.forward(d_x_gpu, d_xf_gpu)
d_xf_gpu.copy_to_host(xf_gpu)
self.assertTrue(np.allclose(xf[:, :, 0:rowsize//2+1], xf_gpu, atol=1e-6))
python类fft()的实例源码
def test_fft_1d_single(self):
from pyculib.fft import fft
N = 32
x = np.asarray(np.arange(N), dtype=np.float32)
xf = np.fft.fft(x)
xf_gpu = np.empty(shape=N//2 + 1, dtype=np.complex64)
fft(x, xf_gpu)
self.assertTrue( np.allclose(xf[0:N//2+1], xf_gpu, atol=1e-6) )
def test_fft_1d_double(self):
from pyculib.fft import fft
N = 32
x = np.asarray(np.arange(N), dtype=np.float64)
xf = np.fft.fft(x)
xf_gpu = np.zeros(shape=N//2 + 1, dtype=np.complex128)
fft(x, xf_gpu)
self.assertTrue( np.allclose(xf[0:N//2+1], xf_gpu, atol=1e-6) )
def test_fft_2d_single(self):
from pyculib.fft import fft
N2 = 2
N1 = 32
N = N1 * N2
x = np.asarray(np.arange(N), dtype=np.float32).reshape(N2, N1)
xf = np.fft.fft2(x)
xf_gpu = np.empty(shape=(N2, N1//2 + 1), dtype=np.complex64)
fft(x, xf_gpu)
self.assertTrue( np.allclose(xf[:, 0:N1//2+1], xf_gpu, atol=1e-6) )
def test_fft_2d_single_col_major(self):
from pyculib.fft import fft
N2 = 2
N1 = 8
N = N1 * N2
x = np.asarray(np.arange(N), dtype=np.float32).reshape((N2, N1), order='F')
xf_ref = np.fft.rfft2(x)
xf = np.empty(shape=(N2, N1//2 + 1), dtype=np.complex64, order='F')
fft(x, xf)
self.assertTrue( np.allclose(xf_ref, xf, atol=1e-6) )
def test_fft_2d_double(self):
from pyculib.fft import fft
N2 = 2
N1 = 32
N = N1 * N2
x = np.asarray(np.arange(N), dtype=np.float64).reshape(N2, N1)
xf = np.fft.fft2(x)
xf_gpu = np.empty(shape=(N2, N1//2 + 1), dtype=np.complex128)
fft(x, xf_gpu)
self.assertTrue( np.allclose(xf[:, 0:N1//2+1], xf_gpu, atol=1e-6) )
def test_fft_3d_single(self):
from pyculib.fft import fft
N3 = 2
N2 = 2
N1 = 32
N = N1 * N2 * N3
x = np.asarray(np.arange(N), dtype=np.float32).reshape(N3, N2, N1)
xf = np.fft.fftn(x)
xf_gpu = np.empty(shape=(N3, N2, N1//2 + 1), dtype=np.complex64)
fft(x, xf_gpu)
self.assertTrue( np.allclose(xf[:, :, 0:N1//2+1], xf_gpu, atol=1e-6) )
def test_fft_3d_double(self):
from pyculib.fft import fft
N3 = 2
N2 = 2
N1 = 32
N = N1 * N2 * N3
x = np.asarray(np.arange(N), dtype=np.float64).reshape(N3, N2, N1)
xf = np.fft.fftn(x)
xf_gpu = np.empty(shape=(N3, N2, N1//2 + 1), dtype=np.complex128)
fft(x, xf_gpu)
self.assertTrue( np.allclose(xf[:, :, 0:N1//2+1], xf_gpu, atol=1e-6) )
def test_fft_1d_roundtrip_single(self):
from pyculib.fft import fft, ifft
N = 32
x = np.asarray(np.arange(N), dtype=np.float32)
x0 = x.copy()
xf_gpu = np.empty(shape=N//2 + 1, dtype=np.complex64)
fft(x, xf_gpu)
ifft(xf_gpu, x)
self.assertTrue( np.allclose(x / N, x0, atol=1e-6) )
def test_fft_1d_roundtrip_double(self):
from pyculib.fft import fft, ifft
N = 32
x = np.asarray(np.arange(N), dtype=np.float64)
x0 = x.copy()
xf_gpu = np.empty(shape=N//2 + 1, dtype=np.complex128)
fft(x, xf_gpu)
ifft(xf_gpu, x)
self.assertTrue( np.allclose(x / N, x0, atol=1e-6) )
def test_fft_3d_roundtrip_single(self):
from pyculib.fft import fft, ifft
N3 = 2
N2 = 2
N1 = 32
N = N3 * N2 * N1
x = np.asarray(np.arange(N), dtype=np.float32).reshape(N3, N2, N1)
x0 = x.copy()
xf_gpu = np.empty(shape=(N3, N2, N1//2 + 1), dtype=np.complex64)
fft(x, xf_gpu)
ifft(xf_gpu, x)
self.assertTrue( np.allclose(x / N, x0, atol=1e-6) )
def test_fft_inplace_1d_single(self):
from pyculib.fft import fft_inplace
N = 32
x = np.asarray(np.arange(N), dtype=np.complex64)
xf = np.fft.fft(x)
fft_inplace(x)
self.assertTrue( np.allclose(xf, x, atol=1e-6) )
def test_fft_inplace_1d_double(self):
from pyculib.fft import fft_inplace
N = 32
x = np.asarray(np.arange(N), dtype=np.complex128)
xf = np.fft.fft(x)
fft_inplace(x)
self.assertTrue( np.allclose(xf, x, atol=1e-6) )
def test_fft_inplace_2d_single(self):
from pyculib.fft import fft_inplace
N1 = 32
N2 = 2
N = N1 * N2
x = np.asarray(np.arange(N), dtype=np.complex64).reshape(N2, N1)
xf = np.fft.fft2(x)
fft_inplace(x)
self.assertTrue( np.allclose(xf, x, atol=1e-6) )
def test_fft_inplace_2d_double(self):
from pyculib.fft import fft_inplace
N1 = 32
N2 = 2
N = N1 * N2
x = np.asarray(np.arange(N), dtype=np.complex128).reshape(N2, N1)
xf = np.fft.fft2(x)
fft_inplace(x)
self.assertTrue( np.allclose(xf, x, atol=1e-6) )
def test_fft_1d_roundtrip_double_2(self):
from pyculib.fft import fft_inplace, ifft_inplace
N = 32
x = np.asarray(np.arange(N), dtype=np.complex128)
x0 = x.copy()
fft_inplace(x)
ifft_inplace(x)
self.assertTrue( np.allclose(x / N, x0, atol=1e-6) )
def test_fft_2d_roundtrip_single_2(self):
from pyculib.fft import fft_inplace, ifft_inplace
N2 = 2
N1 = 32
N = N1 * N2
x = np.asarray(np.arange(N), dtype=np.complex64).reshape(N2, N1)
x0 = x.copy()
fft_inplace(x)
ifft_inplace(x)
self.assertTrue( np.allclose(x / N, x0, atol=1e-6) )
def test_fft_3d_roundtrip_double(self):
from pyculib.fft import fft_inplace, ifft_inplace
N3 = 2
N2 = 2
N1 = 8
N = N3 * N2 * N1
x = np.asarray(np.arange(N), dtype=np.complex128).reshape(N3, N2, N1)
x0 = x.copy()
fft_inplace(x)
ifft_inplace(x)
self.assertTrue( np.allclose(x / N, x0, atol=1e-6) )
def _fft_module(da):
if da.chunks:
return dsar.fft
else:
return np.fft
def fast_autocorrelate(x):
"""
Compute the autocorrelation of the signal, based on the properties of the
power spectral density of the signal.
Note that the input's length may be reduced before the correlation is
performed due to a pathalogical case in numpy.fft:
http://stackoverflow.com/a/23531074/679081
> The FFT algorithm used in np.fft performs very well (meaning O(n log n))
> when the input length has many small prime factors, and very bad
> (meaning a naive DFT requiring O(n^2)) when the input size is a prime
> number.
"""
# This is one simple way to ensure that the input array
# has a length with many small prime factors, although it
# doesn't guarantee that (also hopefully we don't chop too much)
optimal_input_length = int(numpy.sqrt(len(x))) ** 2
x = x[:optimal_input_length]
xp = x - numpy.mean(x)
f = numpy.fft.fft(xp)
p = numpy.absolute(numpy.power(f, 2))
pi = numpy.fft.ifft(p)
result = numpy.real(pi)[:x.size / 2] / numpy.sum(numpy.power(xp, 2))
return result
