test_pocketfft.py 9.6 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261
  1. from __future__ import division, absolute_import, print_function
  2. import numpy as np
  3. import pytest
  4. from numpy.random import random
  5. from numpy.testing import (
  6. assert_array_equal, assert_raises, assert_allclose
  7. )
  8. import threading
  9. import sys
  10. if sys.version_info[0] >= 3:
  11. import queue
  12. else:
  13. import Queue as queue
  14. def fft1(x):
  15. L = len(x)
  16. phase = -2j*np.pi*(np.arange(L)/float(L))
  17. phase = np.arange(L).reshape(-1, 1) * phase
  18. return np.sum(x*np.exp(phase), axis=1)
  19. class TestFFTShift(object):
  20. def test_fft_n(self):
  21. assert_raises(ValueError, np.fft.fft, [1, 2, 3], 0)
  22. class TestFFT1D(object):
  23. def test_identity(self):
  24. maxlen = 512
  25. x = random(maxlen) + 1j*random(maxlen)
  26. xr = random(maxlen)
  27. for i in range(1,maxlen):
  28. assert_allclose(np.fft.ifft(np.fft.fft(x[0:i])), x[0:i],
  29. atol=1e-12)
  30. assert_allclose(np.fft.irfft(np.fft.rfft(xr[0:i]),i),
  31. xr[0:i], atol=1e-12)
  32. def test_fft(self):
  33. x = random(30) + 1j*random(30)
  34. assert_allclose(fft1(x), np.fft.fft(x), atol=1e-6)
  35. assert_allclose(fft1(x) / np.sqrt(30),
  36. np.fft.fft(x, norm="ortho"), atol=1e-6)
  37. @pytest.mark.parametrize('norm', (None, 'ortho'))
  38. def test_ifft(self, norm):
  39. x = random(30) + 1j*random(30)
  40. assert_allclose(
  41. x, np.fft.ifft(np.fft.fft(x, norm=norm), norm=norm),
  42. atol=1e-6)
  43. # Ensure we get the correct error message
  44. with pytest.raises(ValueError,
  45. match='Invalid number of FFT data points'):
  46. np.fft.ifft([], norm=norm)
  47. def test_fft2(self):
  48. x = random((30, 20)) + 1j*random((30, 20))
  49. assert_allclose(np.fft.fft(np.fft.fft(x, axis=1), axis=0),
  50. np.fft.fft2(x), atol=1e-6)
  51. assert_allclose(np.fft.fft2(x) / np.sqrt(30 * 20),
  52. np.fft.fft2(x, norm="ortho"), atol=1e-6)
  53. def test_ifft2(self):
  54. x = random((30, 20)) + 1j*random((30, 20))
  55. assert_allclose(np.fft.ifft(np.fft.ifft(x, axis=1), axis=0),
  56. np.fft.ifft2(x), atol=1e-6)
  57. assert_allclose(np.fft.ifft2(x) * np.sqrt(30 * 20),
  58. np.fft.ifft2(x, norm="ortho"), atol=1e-6)
  59. def test_fftn(self):
  60. x = random((30, 20, 10)) + 1j*random((30, 20, 10))
  61. assert_allclose(
  62. np.fft.fft(np.fft.fft(np.fft.fft(x, axis=2), axis=1), axis=0),
  63. np.fft.fftn(x), atol=1e-6)
  64. assert_allclose(np.fft.fftn(x) / np.sqrt(30 * 20 * 10),
  65. np.fft.fftn(x, norm="ortho"), atol=1e-6)
  66. def test_ifftn(self):
  67. x = random((30, 20, 10)) + 1j*random((30, 20, 10))
  68. assert_allclose(
  69. np.fft.ifft(np.fft.ifft(np.fft.ifft(x, axis=2), axis=1), axis=0),
  70. np.fft.ifftn(x), atol=1e-6)
  71. assert_allclose(np.fft.ifftn(x) * np.sqrt(30 * 20 * 10),
  72. np.fft.ifftn(x, norm="ortho"), atol=1e-6)
  73. def test_rfft(self):
  74. x = random(30)
  75. for n in [x.size, 2*x.size]:
  76. for norm in [None, 'ortho']:
  77. assert_allclose(
  78. np.fft.fft(x, n=n, norm=norm)[:(n//2 + 1)],
  79. np.fft.rfft(x, n=n, norm=norm), atol=1e-6)
  80. assert_allclose(
  81. np.fft.rfft(x, n=n) / np.sqrt(n),
  82. np.fft.rfft(x, n=n, norm="ortho"), atol=1e-6)
  83. def test_irfft(self):
  84. x = random(30)
  85. assert_allclose(x, np.fft.irfft(np.fft.rfft(x)), atol=1e-6)
  86. assert_allclose(
  87. x, np.fft.irfft(np.fft.rfft(x, norm="ortho"), norm="ortho"), atol=1e-6)
  88. def test_rfft2(self):
  89. x = random((30, 20))
  90. assert_allclose(np.fft.fft2(x)[:, :11], np.fft.rfft2(x), atol=1e-6)
  91. assert_allclose(np.fft.rfft2(x) / np.sqrt(30 * 20),
  92. np.fft.rfft2(x, norm="ortho"), atol=1e-6)
  93. def test_irfft2(self):
  94. x = random((30, 20))
  95. assert_allclose(x, np.fft.irfft2(np.fft.rfft2(x)), atol=1e-6)
  96. assert_allclose(
  97. x, np.fft.irfft2(np.fft.rfft2(x, norm="ortho"), norm="ortho"), atol=1e-6)
  98. def test_rfftn(self):
  99. x = random((30, 20, 10))
  100. assert_allclose(np.fft.fftn(x)[:, :, :6], np.fft.rfftn(x), atol=1e-6)
  101. assert_allclose(np.fft.rfftn(x) / np.sqrt(30 * 20 * 10),
  102. np.fft.rfftn(x, norm="ortho"), atol=1e-6)
  103. def test_irfftn(self):
  104. x = random((30, 20, 10))
  105. assert_allclose(x, np.fft.irfftn(np.fft.rfftn(x)), atol=1e-6)
  106. assert_allclose(
  107. x, np.fft.irfftn(np.fft.rfftn(x, norm="ortho"), norm="ortho"), atol=1e-6)
  108. def test_hfft(self):
  109. x = random(14) + 1j*random(14)
  110. x_herm = np.concatenate((random(1), x, random(1)))
  111. x = np.concatenate((x_herm, x[::-1].conj()))
  112. assert_allclose(np.fft.fft(x), np.fft.hfft(x_herm), atol=1e-6)
  113. assert_allclose(np.fft.hfft(x_herm) / np.sqrt(30),
  114. np.fft.hfft(x_herm, norm="ortho"), atol=1e-6)
  115. def test_ihttf(self):
  116. x = random(14) + 1j*random(14)
  117. x_herm = np.concatenate((random(1), x, random(1)))
  118. x = np.concatenate((x_herm, x[::-1].conj()))
  119. assert_allclose(x_herm, np.fft.ihfft(np.fft.hfft(x_herm)), atol=1e-6)
  120. assert_allclose(
  121. x_herm, np.fft.ihfft(np.fft.hfft(x_herm, norm="ortho"),
  122. norm="ortho"), atol=1e-6)
  123. @pytest.mark.parametrize("op", [np.fft.fftn, np.fft.ifftn,
  124. np.fft.rfftn, np.fft.irfftn])
  125. def test_axes(self, op):
  126. x = random((30, 20, 10))
  127. axes = [(0, 1, 2), (0, 2, 1), (1, 0, 2), (1, 2, 0), (2, 0, 1), (2, 1, 0)]
  128. for a in axes:
  129. op_tr = op(np.transpose(x, a))
  130. tr_op = np.transpose(op(x, axes=a), a)
  131. assert_allclose(op_tr, tr_op, atol=1e-6)
  132. def test_all_1d_norm_preserving(self):
  133. # verify that round-trip transforms are norm-preserving
  134. x = random(30)
  135. x_norm = np.linalg.norm(x)
  136. n = x.size * 2
  137. func_pairs = [(np.fft.fft, np.fft.ifft),
  138. (np.fft.rfft, np.fft.irfft),
  139. # hfft: order so the first function takes x.size samples
  140. # (necessary for comparison to x_norm above)
  141. (np.fft.ihfft, np.fft.hfft),
  142. ]
  143. for forw, back in func_pairs:
  144. for n in [x.size, 2*x.size]:
  145. for norm in [None, 'ortho']:
  146. tmp = forw(x, n=n, norm=norm)
  147. tmp = back(tmp, n=n, norm=norm)
  148. assert_allclose(x_norm,
  149. np.linalg.norm(tmp), atol=1e-6)
  150. @pytest.mark.parametrize("dtype", [np.half, np.single, np.double,
  151. np.longdouble])
  152. def test_dtypes(self, dtype):
  153. # make sure that all input precisions are accepted and internally
  154. # converted to 64bit
  155. x = random(30).astype(dtype)
  156. assert_allclose(np.fft.ifft(np.fft.fft(x)), x, atol=1e-6)
  157. assert_allclose(np.fft.irfft(np.fft.rfft(x)), x, atol=1e-6)
  158. @pytest.mark.parametrize(
  159. "dtype",
  160. [np.float32, np.float64, np.complex64, np.complex128])
  161. @pytest.mark.parametrize("order", ["F", 'non-contiguous'])
  162. @pytest.mark.parametrize(
  163. "fft",
  164. [np.fft.fft, np.fft.fft2, np.fft.fftn,
  165. np.fft.ifft, np.fft.ifft2, np.fft.ifftn])
  166. def test_fft_with_order(dtype, order, fft):
  167. # Check that FFT/IFFT produces identical results for C, Fortran and
  168. # non contiguous arrays
  169. rng = np.random.RandomState(42)
  170. X = rng.rand(8, 7, 13).astype(dtype, copy=False)
  171. # See discussion in pull/14178
  172. _tol = 8.0 * np.sqrt(np.log2(X.size)) * np.finfo(X.dtype).eps
  173. if order == 'F':
  174. Y = np.asfortranarray(X)
  175. else:
  176. # Make a non contiguous array
  177. Y = X[::-1]
  178. X = np.ascontiguousarray(X[::-1])
  179. if fft.__name__.endswith('fft'):
  180. for axis in range(3):
  181. X_res = fft(X, axis=axis)
  182. Y_res = fft(Y, axis=axis)
  183. assert_allclose(X_res, Y_res, atol=_tol, rtol=_tol)
  184. elif fft.__name__.endswith(('fft2', 'fftn')):
  185. axes = [(0, 1), (1, 2), (0, 2)]
  186. if fft.__name__.endswith('fftn'):
  187. axes.extend([(0,), (1,), (2,), None])
  188. for ax in axes:
  189. X_res = fft(X, axes=ax)
  190. Y_res = fft(Y, axes=ax)
  191. assert_allclose(X_res, Y_res, atol=_tol, rtol=_tol)
  192. else:
  193. raise ValueError()
  194. class TestFFTThreadSafe(object):
  195. threads = 16
  196. input_shape = (800, 200)
  197. def _test_mtsame(self, func, *args):
  198. def worker(args, q):
  199. q.put(func(*args))
  200. q = queue.Queue()
  201. expected = func(*args)
  202. # Spin off a bunch of threads to call the same function simultaneously
  203. t = [threading.Thread(target=worker, args=(args, q))
  204. for i in range(self.threads)]
  205. [x.start() for x in t]
  206. [x.join() for x in t]
  207. # Make sure all threads returned the correct value
  208. for i in range(self.threads):
  209. assert_array_equal(q.get(timeout=5), expected,
  210. 'Function returned wrong value in multithreaded context')
  211. def test_fft(self):
  212. a = np.ones(self.input_shape) * 1+0j
  213. self._test_mtsame(np.fft.fft, a)
  214. def test_ifft(self):
  215. a = np.ones(self.input_shape) * 1+0j
  216. self._test_mtsame(np.fft.ifft, a)
  217. def test_rfft(self):
  218. a = np.ones(self.input_shape)
  219. self._test_mtsame(np.fft.rfft, a)
  220. def test_irfft(self):
  221. a = np.ones(self.input_shape) * 1+0j
  222. self._test_mtsame(np.fft.irfft, a)