test_basic.py 33 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913914915916917918919920921922923924925926927928929930931932933934935936937938939940941942943944945946947948949950951952953954955956957958959960961962963964965966967968969970971972
  1. # Created by Pearu Peterson, September 2002
  2. from __future__ import division, print_function, absolute_import
  3. __usage__ = """
  4. Build fftpack:
  5. python setup_fftpack.py build
  6. Run tests if scipy is installed:
  7. python -c 'import scipy;scipy.fftpack.test()'
  8. Run tests if fftpack is not installed:
  9. python tests/test_basic.py
  10. """
  11. from numpy.testing import (assert_, assert_equal, assert_array_almost_equal,
  12. assert_array_almost_equal_nulp, assert_array_less)
  13. import pytest
  14. from pytest import raises as assert_raises
  15. from scipy.fftpack import ifft, fft, fftn, ifftn, rfft, irfft, fft2
  16. from scipy.fftpack import _fftpack as fftpack
  17. from scipy.fftpack.basic import _is_safe_size
  18. from numpy import (arange, add, array, asarray, zeros, dot, exp, pi,
  19. swapaxes, double, cdouble)
  20. import numpy as np
  21. import numpy.fft
  22. from numpy.random import rand
  23. # "large" composite numbers supported by FFTPACK
  24. LARGE_COMPOSITE_SIZES = [
  25. 2**13,
  26. 2**5 * 3**5,
  27. 2**3 * 3**3 * 5**2,
  28. ]
  29. SMALL_COMPOSITE_SIZES = [
  30. 2,
  31. 2*3*5,
  32. 2*2*3*3,
  33. ]
  34. # prime
  35. LARGE_PRIME_SIZES = [
  36. 2011
  37. ]
  38. SMALL_PRIME_SIZES = [
  39. 29
  40. ]
  41. def _assert_close_in_norm(x, y, rtol, size, rdt):
  42. # helper function for testing
  43. err_msg = "size: %s rdt: %s" % (size, rdt)
  44. assert_array_less(np.linalg.norm(x - y), rtol*np.linalg.norm(x), err_msg)
  45. def random(size):
  46. return rand(*size)
  47. def get_mat(n):
  48. data = arange(n)
  49. data = add.outer(data, data)
  50. return data
  51. def direct_dft(x):
  52. x = asarray(x)
  53. n = len(x)
  54. y = zeros(n, dtype=cdouble)
  55. w = -arange(n)*(2j*pi/n)
  56. for i in range(n):
  57. y[i] = dot(exp(i*w), x)
  58. return y
  59. def direct_idft(x):
  60. x = asarray(x)
  61. n = len(x)
  62. y = zeros(n, dtype=cdouble)
  63. w = arange(n)*(2j*pi/n)
  64. for i in range(n):
  65. y[i] = dot(exp(i*w), x)/n
  66. return y
  67. def direct_dftn(x):
  68. x = asarray(x)
  69. for axis in range(len(x.shape)):
  70. x = fft(x, axis=axis)
  71. return x
  72. def direct_idftn(x):
  73. x = asarray(x)
  74. for axis in range(len(x.shape)):
  75. x = ifft(x, axis=axis)
  76. return x
  77. def direct_rdft(x):
  78. x = asarray(x)
  79. n = len(x)
  80. w = -arange(n)*(2j*pi/n)
  81. r = zeros(n, dtype=double)
  82. for i in range(n//2+1):
  83. y = dot(exp(i*w), x)
  84. if i:
  85. r[2*i-1] = y.real
  86. if 2*i < n:
  87. r[2*i] = y.imag
  88. else:
  89. r[0] = y.real
  90. return r
  91. def direct_irdft(x):
  92. x = asarray(x)
  93. n = len(x)
  94. x1 = zeros(n, dtype=cdouble)
  95. for i in range(n//2+1):
  96. if i:
  97. if 2*i < n:
  98. x1[i] = x[2*i-1] + 1j*x[2*i]
  99. x1[n-i] = x[2*i-1] - 1j*x[2*i]
  100. else:
  101. x1[i] = x[2*i-1]
  102. else:
  103. x1[0] = x[0]
  104. return direct_idft(x1).real
  105. class _TestFFTBase(object):
  106. def setup_method(self):
  107. self.cdt = None
  108. self.rdt = None
  109. np.random.seed(1234)
  110. def test_definition(self):
  111. x = np.array([1,2,3,4+1j,1,2,3,4+2j], dtype=self.cdt)
  112. y = fft(x)
  113. assert_equal(y.dtype, self.cdt)
  114. y1 = direct_dft(x)
  115. assert_array_almost_equal(y,y1)
  116. x = np.array([1,2,3,4+0j,5], dtype=self.cdt)
  117. assert_array_almost_equal(fft(x),direct_dft(x))
  118. def test_n_argument_real(self):
  119. x1 = np.array([1,2,3,4], dtype=self.rdt)
  120. x2 = np.array([1,2,3,4], dtype=self.rdt)
  121. y = fft([x1,x2],n=4)
  122. assert_equal(y.dtype, self.cdt)
  123. assert_equal(y.shape,(2,4))
  124. assert_array_almost_equal(y[0],direct_dft(x1))
  125. assert_array_almost_equal(y[1],direct_dft(x2))
  126. def _test_n_argument_complex(self):
  127. x1 = np.array([1,2,3,4+1j], dtype=self.cdt)
  128. x2 = np.array([1,2,3,4+1j], dtype=self.cdt)
  129. y = fft([x1,x2],n=4)
  130. assert_equal(y.dtype, self.cdt)
  131. assert_equal(y.shape,(2,4))
  132. assert_array_almost_equal(y[0],direct_dft(x1))
  133. assert_array_almost_equal(y[1],direct_dft(x2))
  134. def test_djbfft(self):
  135. for i in range(2,14):
  136. n = 2**i
  137. x = list(range(n))
  138. y = fftpack.zfft(x)
  139. y2 = numpy.fft.fft(x)
  140. assert_array_almost_equal(y,y2)
  141. y = fftpack.zrfft(x)
  142. assert_array_almost_equal(y,y2)
  143. def test_invalid_sizes(self):
  144. assert_raises(ValueError, fft, [])
  145. assert_raises(ValueError, fft, [[1,1],[2,2]], -5)
  146. def test__is_safe_size(self):
  147. vals = [(0, True), (1, True), (2, True), (3, True), (4, True), (5, True), (6, True), (7, False),
  148. (15, True), (16, True), (17, False), (18, True), (21, False), (25, True), (50, True),
  149. (120, True), (210, False)]
  150. for n, is_safe in vals:
  151. assert_equal(_is_safe_size(n), is_safe)
  152. class TestDoubleFFT(_TestFFTBase):
  153. def setup_method(self):
  154. self.cdt = np.cdouble
  155. self.rdt = np.double
  156. class TestSingleFFT(_TestFFTBase):
  157. def setup_method(self):
  158. self.cdt = np.complex64
  159. self.rdt = np.float32
  160. @pytest.mark.xfail(run=False, reason="single-precision FFT implementation is partially disabled, until accuracy issues with large prime powers are resolved")
  161. def test_notice(self):
  162. pass
  163. class TestFloat16FFT(object):
  164. def test_1_argument_real(self):
  165. x1 = np.array([1, 2, 3, 4], dtype=np.float16)
  166. y = fft(x1, n=4)
  167. assert_equal(y.dtype, np.complex64)
  168. assert_equal(y.shape, (4, ))
  169. assert_array_almost_equal(y, direct_dft(x1.astype(np.float32)))
  170. def test_n_argument_real(self):
  171. x1 = np.array([1, 2, 3, 4], dtype=np.float16)
  172. x2 = np.array([1, 2, 3, 4], dtype=np.float16)
  173. y = fft([x1, x2], n=4)
  174. assert_equal(y.dtype, np.complex64)
  175. assert_equal(y.shape, (2, 4))
  176. assert_array_almost_equal(y[0], direct_dft(x1.astype(np.float32)))
  177. assert_array_almost_equal(y[1], direct_dft(x2.astype(np.float32)))
  178. class _TestIFFTBase(object):
  179. def setup_method(self):
  180. np.random.seed(1234)
  181. def test_definition(self):
  182. x = np.array([1,2,3,4+1j,1,2,3,4+2j], self.cdt)
  183. y = ifft(x)
  184. y1 = direct_idft(x)
  185. assert_equal(y.dtype, self.cdt)
  186. assert_array_almost_equal(y,y1)
  187. x = np.array([1,2,3,4+0j,5], self.cdt)
  188. assert_array_almost_equal(ifft(x),direct_idft(x))
  189. def test_definition_real(self):
  190. x = np.array([1,2,3,4,1,2,3,4], self.rdt)
  191. y = ifft(x)
  192. assert_equal(y.dtype, self.cdt)
  193. y1 = direct_idft(x)
  194. assert_array_almost_equal(y,y1)
  195. x = np.array([1,2,3,4,5], dtype=self.rdt)
  196. assert_equal(y.dtype, self.cdt)
  197. assert_array_almost_equal(ifft(x),direct_idft(x))
  198. def test_djbfft(self):
  199. for i in range(2,14):
  200. n = 2**i
  201. x = list(range(n))
  202. y = fftpack.zfft(x,direction=-1)
  203. y2 = numpy.fft.ifft(x)
  204. assert_array_almost_equal(y,y2)
  205. y = fftpack.zrfft(x,direction=-1)
  206. assert_array_almost_equal(y,y2)
  207. def test_random_complex(self):
  208. for size in [1,51,111,100,200,64,128,256,1024]:
  209. x = random([size]).astype(self.cdt)
  210. x = random([size]).astype(self.cdt) + 1j*x
  211. y1 = ifft(fft(x))
  212. y2 = fft(ifft(x))
  213. assert_equal(y1.dtype, self.cdt)
  214. assert_equal(y2.dtype, self.cdt)
  215. assert_array_almost_equal(y1, x)
  216. assert_array_almost_equal(y2, x)
  217. def test_random_real(self):
  218. for size in [1,51,111,100,200,64,128,256,1024]:
  219. x = random([size]).astype(self.rdt)
  220. y1 = ifft(fft(x))
  221. y2 = fft(ifft(x))
  222. assert_equal(y1.dtype, self.cdt)
  223. assert_equal(y2.dtype, self.cdt)
  224. assert_array_almost_equal(y1, x)
  225. assert_array_almost_equal(y2, x)
  226. def test_size_accuracy(self):
  227. # Sanity check for the accuracy for prime and non-prime sized inputs
  228. if self.rdt == np.float32:
  229. rtol = 1e-5
  230. elif self.rdt == np.float64:
  231. rtol = 1e-10
  232. for size in LARGE_COMPOSITE_SIZES + LARGE_PRIME_SIZES:
  233. np.random.seed(1234)
  234. x = np.random.rand(size).astype(self.rdt)
  235. y = ifft(fft(x))
  236. _assert_close_in_norm(x, y, rtol, size, self.rdt)
  237. y = fft(ifft(x))
  238. _assert_close_in_norm(x, y, rtol, size, self.rdt)
  239. x = (x + 1j*np.random.rand(size)).astype(self.cdt)
  240. y = ifft(fft(x))
  241. _assert_close_in_norm(x, y, rtol, size, self.rdt)
  242. y = fft(ifft(x))
  243. _assert_close_in_norm(x, y, rtol, size, self.rdt)
  244. def test_invalid_sizes(self):
  245. assert_raises(ValueError, ifft, [])
  246. assert_raises(ValueError, ifft, [[1,1],[2,2]], -5)
  247. class TestDoubleIFFT(_TestIFFTBase):
  248. def setup_method(self):
  249. self.cdt = np.cdouble
  250. self.rdt = np.double
  251. class TestSingleIFFT(_TestIFFTBase):
  252. def setup_method(self):
  253. self.cdt = np.complex64
  254. self.rdt = np.float32
  255. class _TestRFFTBase(object):
  256. def setup_method(self):
  257. np.random.seed(1234)
  258. def test_definition(self):
  259. for t in [[1, 2, 3, 4, 1, 2, 3, 4], [1, 2, 3, 4, 1, 2, 3, 4, 5]]:
  260. x = np.array(t, dtype=self.rdt)
  261. y = rfft(x)
  262. y1 = direct_rdft(x)
  263. assert_array_almost_equal(y,y1)
  264. assert_equal(y.dtype, self.rdt)
  265. def test_djbfft(self):
  266. from numpy.fft import fft as numpy_fft
  267. for i in range(2,14):
  268. n = 2**i
  269. x = list(range(n))
  270. y2 = numpy_fft(x)
  271. y1 = zeros((n,),dtype=double)
  272. y1[0] = y2[0].real
  273. y1[-1] = y2[n//2].real
  274. for k in range(1, n//2):
  275. y1[2*k-1] = y2[k].real
  276. y1[2*k] = y2[k].imag
  277. y = fftpack.drfft(x)
  278. assert_array_almost_equal(y,y1)
  279. def test_invalid_sizes(self):
  280. assert_raises(ValueError, rfft, [])
  281. assert_raises(ValueError, rfft, [[1,1],[2,2]], -5)
  282. # See gh-5790
  283. class MockSeries(object):
  284. def __init__(self, data):
  285. self.data = np.asarray(data)
  286. def __getattr__(self, item):
  287. try:
  288. return getattr(self.data, item)
  289. except AttributeError:
  290. raise AttributeError(("'MockSeries' object "
  291. "has no attribute '{attr}'".
  292. format(attr=item)))
  293. def test_non_ndarray_with_dtype(self):
  294. x = np.array([1., 2., 3., 4., 5.])
  295. xs = _TestRFFTBase.MockSeries(x)
  296. expected = [1, 2, 3, 4, 5]
  297. out = rfft(xs)
  298. # Data should not have been overwritten
  299. assert_equal(x, expected)
  300. assert_equal(xs.data, expected)
  301. class TestRFFTDouble(_TestRFFTBase):
  302. def setup_method(self):
  303. self.cdt = np.cdouble
  304. self.rdt = np.double
  305. class TestRFFTSingle(_TestRFFTBase):
  306. def setup_method(self):
  307. self.cdt = np.complex64
  308. self.rdt = np.float32
  309. class _TestIRFFTBase(object):
  310. def setup_method(self):
  311. np.random.seed(1234)
  312. def test_definition(self):
  313. x1 = [1,2,3,4,1,2,3,4]
  314. x1_1 = [1,2+3j,4+1j,2+3j,4,2-3j,4-1j,2-3j]
  315. x2 = [1,2,3,4,1,2,3,4,5]
  316. x2_1 = [1,2+3j,4+1j,2+3j,4+5j,4-5j,2-3j,4-1j,2-3j]
  317. def _test(x, xr):
  318. y = irfft(np.array(x, dtype=self.rdt))
  319. y1 = direct_irdft(x)
  320. assert_equal(y.dtype, self.rdt)
  321. assert_array_almost_equal(y,y1, decimal=self.ndec)
  322. assert_array_almost_equal(y,ifft(xr), decimal=self.ndec)
  323. _test(x1, x1_1)
  324. _test(x2, x2_1)
  325. def test_djbfft(self):
  326. from numpy.fft import ifft as numpy_ifft
  327. for i in range(2,14):
  328. n = 2**i
  329. x = list(range(n))
  330. x1 = zeros((n,),dtype=cdouble)
  331. x1[0] = x[0]
  332. for k in range(1, n//2):
  333. x1[k] = x[2*k-1]+1j*x[2*k]
  334. x1[n-k] = x[2*k-1]-1j*x[2*k]
  335. x1[n//2] = x[-1]
  336. y1 = numpy_ifft(x1)
  337. y = fftpack.drfft(x,direction=-1)
  338. assert_array_almost_equal(y,y1)
  339. def test_random_real(self):
  340. for size in [1,51,111,100,200,64,128,256,1024]:
  341. x = random([size]).astype(self.rdt)
  342. y1 = irfft(rfft(x))
  343. y2 = rfft(irfft(x))
  344. assert_equal(y1.dtype, self.rdt)
  345. assert_equal(y2.dtype, self.rdt)
  346. assert_array_almost_equal(y1, x, decimal=self.ndec,
  347. err_msg="size=%d" % size)
  348. assert_array_almost_equal(y2, x, decimal=self.ndec,
  349. err_msg="size=%d" % size)
  350. def test_size_accuracy(self):
  351. # Sanity check for the accuracy for prime and non-prime sized inputs
  352. if self.rdt == np.float32:
  353. rtol = 1e-5
  354. elif self.rdt == np.float64:
  355. rtol = 1e-10
  356. for size in LARGE_COMPOSITE_SIZES + LARGE_PRIME_SIZES:
  357. np.random.seed(1234)
  358. x = np.random.rand(size).astype(self.rdt)
  359. y = irfft(rfft(x))
  360. _assert_close_in_norm(x, y, rtol, size, self.rdt)
  361. y = rfft(irfft(x))
  362. _assert_close_in_norm(x, y, rtol, size, self.rdt)
  363. def test_invalid_sizes(self):
  364. assert_raises(ValueError, irfft, [])
  365. assert_raises(ValueError, irfft, [[1,1],[2,2]], -5)
  366. # self.ndec is bogus; we should have a assert_array_approx_equal for number of
  367. # significant digits
  368. class TestIRFFTDouble(_TestIRFFTBase):
  369. def setup_method(self):
  370. self.cdt = np.cdouble
  371. self.rdt = np.double
  372. self.ndec = 14
  373. class TestIRFFTSingle(_TestIRFFTBase):
  374. def setup_method(self):
  375. self.cdt = np.complex64
  376. self.rdt = np.float32
  377. self.ndec = 5
  378. class Testfft2(object):
  379. def setup_method(self):
  380. np.random.seed(1234)
  381. def test_regression_244(self):
  382. """FFT returns wrong result with axes parameter."""
  383. # fftn (and hence fft2) used to break when both axes and shape were
  384. # used
  385. x = numpy.ones((4, 4, 2))
  386. y = fft2(x, shape=(8, 8), axes=(-3, -2))
  387. y_r = numpy.fft.fftn(x, s=(8, 8), axes=(-3, -2))
  388. assert_array_almost_equal(y, y_r)
  389. def test_invalid_sizes(self):
  390. assert_raises(ValueError, fft2, [[]])
  391. assert_raises(ValueError, fft2, [[1, 1], [2, 2]], (4, -3))
  392. class TestFftnSingle(object):
  393. def setup_method(self):
  394. np.random.seed(1234)
  395. def test_definition(self):
  396. x = [[1, 2, 3],
  397. [4, 5, 6],
  398. [7, 8, 9]]
  399. y = fftn(np.array(x, np.float32))
  400. assert_(y.dtype == np.complex64,
  401. msg="double precision output with single precision")
  402. y_r = np.array(fftn(x), np.complex64)
  403. assert_array_almost_equal_nulp(y, y_r)
  404. @pytest.mark.parametrize('size', SMALL_COMPOSITE_SIZES + SMALL_PRIME_SIZES)
  405. def test_size_accuracy_small(self, size):
  406. x = np.random.rand(size, size) + 1j*np.random.rand(size, size)
  407. y1 = fftn(x.real.astype(np.float32))
  408. y2 = fftn(x.real.astype(np.float64)).astype(np.complex64)
  409. assert_equal(y1.dtype, np.complex64)
  410. assert_array_almost_equal_nulp(y1, y2, 2000)
  411. @pytest.mark.parametrize('size', LARGE_COMPOSITE_SIZES + LARGE_PRIME_SIZES)
  412. def test_size_accuracy_large(self, size):
  413. x = np.random.rand(size, 3) + 1j*np.random.rand(size, 3)
  414. y1 = fftn(x.real.astype(np.float32))
  415. y2 = fftn(x.real.astype(np.float64)).astype(np.complex64)
  416. assert_equal(y1.dtype, np.complex64)
  417. assert_array_almost_equal_nulp(y1, y2, 2000)
  418. def test_definition_float16(self):
  419. x = [[1, 2, 3],
  420. [4, 5, 6],
  421. [7, 8, 9]]
  422. y = fftn(np.array(x, np.float16))
  423. assert_equal(y.dtype, np.complex64)
  424. y_r = np.array(fftn(x), np.complex64)
  425. assert_array_almost_equal_nulp(y, y_r)
  426. @pytest.mark.parametrize('size', SMALL_COMPOSITE_SIZES + SMALL_PRIME_SIZES)
  427. def test_float16_input_small(self, size):
  428. x = np.random.rand(size, size) + 1j*np.random.rand(size, size)
  429. y1 = fftn(x.real.astype(np.float16))
  430. y2 = fftn(x.real.astype(np.float64)).astype(np.complex64)
  431. assert_equal(y1.dtype, np.complex64)
  432. assert_array_almost_equal_nulp(y1, y2, 5e5)
  433. @pytest.mark.parametrize('size', LARGE_COMPOSITE_SIZES + LARGE_PRIME_SIZES)
  434. def test_float16_input_large(self, size):
  435. x = np.random.rand(size, 3) + 1j*np.random.rand(size, 3)
  436. y1 = fftn(x.real.astype(np.float16))
  437. y2 = fftn(x.real.astype(np.float64)).astype(np.complex64)
  438. assert_equal(y1.dtype, np.complex64)
  439. assert_array_almost_equal_nulp(y1, y2, 2e6)
  440. class TestFftn(object):
  441. def setup_method(self):
  442. np.random.seed(1234)
  443. def test_definition(self):
  444. x = [[1, 2, 3],
  445. [4, 5, 6],
  446. [7, 8, 9]]
  447. y = fftn(x)
  448. assert_array_almost_equal(y, direct_dftn(x))
  449. x = random((20, 26))
  450. assert_array_almost_equal(fftn(x), direct_dftn(x))
  451. x = random((5, 4, 3, 20))
  452. assert_array_almost_equal(fftn(x), direct_dftn(x))
  453. def test_axes_argument(self):
  454. # plane == ji_plane, x== kji_space
  455. plane1 = [[1, 2, 3],
  456. [4, 5, 6],
  457. [7, 8, 9]]
  458. plane2 = [[10, 11, 12],
  459. [13, 14, 15],
  460. [16, 17, 18]]
  461. plane3 = [[19, 20, 21],
  462. [22, 23, 24],
  463. [25, 26, 27]]
  464. ki_plane1 = [[1, 2, 3],
  465. [10, 11, 12],
  466. [19, 20, 21]]
  467. ki_plane2 = [[4, 5, 6],
  468. [13, 14, 15],
  469. [22, 23, 24]]
  470. ki_plane3 = [[7, 8, 9],
  471. [16, 17, 18],
  472. [25, 26, 27]]
  473. jk_plane1 = [[1, 10, 19],
  474. [4, 13, 22],
  475. [7, 16, 25]]
  476. jk_plane2 = [[2, 11, 20],
  477. [5, 14, 23],
  478. [8, 17, 26]]
  479. jk_plane3 = [[3, 12, 21],
  480. [6, 15, 24],
  481. [9, 18, 27]]
  482. kj_plane1 = [[1, 4, 7],
  483. [10, 13, 16], [19, 22, 25]]
  484. kj_plane2 = [[2, 5, 8],
  485. [11, 14, 17], [20, 23, 26]]
  486. kj_plane3 = [[3, 6, 9],
  487. [12, 15, 18], [21, 24, 27]]
  488. ij_plane1 = [[1, 4, 7],
  489. [2, 5, 8],
  490. [3, 6, 9]]
  491. ij_plane2 = [[10, 13, 16],
  492. [11, 14, 17],
  493. [12, 15, 18]]
  494. ij_plane3 = [[19, 22, 25],
  495. [20, 23, 26],
  496. [21, 24, 27]]
  497. ik_plane1 = [[1, 10, 19],
  498. [2, 11, 20],
  499. [3, 12, 21]]
  500. ik_plane2 = [[4, 13, 22],
  501. [5, 14, 23],
  502. [6, 15, 24]]
  503. ik_plane3 = [[7, 16, 25],
  504. [8, 17, 26],
  505. [9, 18, 27]]
  506. ijk_space = [jk_plane1, jk_plane2, jk_plane3]
  507. ikj_space = [kj_plane1, kj_plane2, kj_plane3]
  508. jik_space = [ik_plane1, ik_plane2, ik_plane3]
  509. jki_space = [ki_plane1, ki_plane2, ki_plane3]
  510. kij_space = [ij_plane1, ij_plane2, ij_plane3]
  511. x = array([plane1, plane2, plane3])
  512. assert_array_almost_equal(fftn(x),
  513. fftn(x, axes=(-3, -2, -1))) # kji_space
  514. assert_array_almost_equal(fftn(x), fftn(x, axes=(0, 1, 2)))
  515. assert_array_almost_equal(fftn(x, axes=(0, 2)), fftn(x, axes=(0, -1)))
  516. y = fftn(x, axes=(2, 1, 0)) # ijk_space
  517. assert_array_almost_equal(swapaxes(y, -1, -3), fftn(ijk_space))
  518. y = fftn(x, axes=(2, 0, 1)) # ikj_space
  519. assert_array_almost_equal(swapaxes(swapaxes(y, -1, -3), -1, -2),
  520. fftn(ikj_space))
  521. y = fftn(x, axes=(1, 2, 0)) # jik_space
  522. assert_array_almost_equal(swapaxes(swapaxes(y, -1, -3), -3, -2),
  523. fftn(jik_space))
  524. y = fftn(x, axes=(1, 0, 2)) # jki_space
  525. assert_array_almost_equal(swapaxes(y, -2, -3), fftn(jki_space))
  526. y = fftn(x, axes=(0, 2, 1)) # kij_space
  527. assert_array_almost_equal(swapaxes(y, -2, -1), fftn(kij_space))
  528. y = fftn(x, axes=(-2, -1)) # ji_plane
  529. assert_array_almost_equal(fftn(plane1), y[0])
  530. assert_array_almost_equal(fftn(plane2), y[1])
  531. assert_array_almost_equal(fftn(plane3), y[2])
  532. y = fftn(x, axes=(1, 2)) # ji_plane
  533. assert_array_almost_equal(fftn(plane1), y[0])
  534. assert_array_almost_equal(fftn(plane2), y[1])
  535. assert_array_almost_equal(fftn(plane3), y[2])
  536. y = fftn(x, axes=(-3, -2)) # kj_plane
  537. assert_array_almost_equal(fftn(x[:, :, 0]), y[:, :, 0])
  538. assert_array_almost_equal(fftn(x[:, :, 1]), y[:, :, 1])
  539. assert_array_almost_equal(fftn(x[:, :, 2]), y[:, :, 2])
  540. y = fftn(x, axes=(-3, -1)) # ki_plane
  541. assert_array_almost_equal(fftn(x[:, 0, :]), y[:, 0, :])
  542. assert_array_almost_equal(fftn(x[:, 1, :]), y[:, 1, :])
  543. assert_array_almost_equal(fftn(x[:, 2, :]), y[:, 2, :])
  544. y = fftn(x, axes=(-1, -2)) # ij_plane
  545. assert_array_almost_equal(fftn(ij_plane1), swapaxes(y[0], -2, -1))
  546. assert_array_almost_equal(fftn(ij_plane2), swapaxes(y[1], -2, -1))
  547. assert_array_almost_equal(fftn(ij_plane3), swapaxes(y[2], -2, -1))
  548. y = fftn(x, axes=(-1, -3)) # ik_plane
  549. assert_array_almost_equal(fftn(ik_plane1),
  550. swapaxes(y[:, 0, :], -1, -2))
  551. assert_array_almost_equal(fftn(ik_plane2),
  552. swapaxes(y[:, 1, :], -1, -2))
  553. assert_array_almost_equal(fftn(ik_plane3),
  554. swapaxes(y[:, 2, :], -1, -2))
  555. y = fftn(x, axes=(-2, -3)) # jk_plane
  556. assert_array_almost_equal(fftn(jk_plane1),
  557. swapaxes(y[:, :, 0], -1, -2))
  558. assert_array_almost_equal(fftn(jk_plane2),
  559. swapaxes(y[:, :, 1], -1, -2))
  560. assert_array_almost_equal(fftn(jk_plane3),
  561. swapaxes(y[:, :, 2], -1, -2))
  562. y = fftn(x, axes=(-1,)) # i_line
  563. for i in range(3):
  564. for j in range(3):
  565. assert_array_almost_equal(fft(x[i, j, :]), y[i, j, :])
  566. y = fftn(x, axes=(-2,)) # j_line
  567. for i in range(3):
  568. for j in range(3):
  569. assert_array_almost_equal(fft(x[i, :, j]), y[i, :, j])
  570. y = fftn(x, axes=(0,)) # k_line
  571. for i in range(3):
  572. for j in range(3):
  573. assert_array_almost_equal(fft(x[:, i, j]), y[:, i, j])
  574. y = fftn(x, axes=()) # point
  575. assert_array_almost_equal(y, x)
  576. def test_shape_argument(self):
  577. small_x = [[1, 2, 3],
  578. [4, 5, 6]]
  579. large_x1 = [[1, 2, 3, 0],
  580. [4, 5, 6, 0],
  581. [0, 0, 0, 0],
  582. [0, 0, 0, 0]]
  583. y = fftn(small_x, shape=(4, 4))
  584. assert_array_almost_equal(y, fftn(large_x1))
  585. y = fftn(small_x, shape=(3, 4))
  586. assert_array_almost_equal(y, fftn(large_x1[:-1]))
  587. def test_shape_axes_argument(self):
  588. small_x = [[1, 2, 3],
  589. [4, 5, 6],
  590. [7, 8, 9]]
  591. large_x1 = array([[1, 2, 3, 0],
  592. [4, 5, 6, 0],
  593. [7, 8, 9, 0],
  594. [0, 0, 0, 0]])
  595. y = fftn(small_x, shape=(4, 4), axes=(-2, -1))
  596. assert_array_almost_equal(y, fftn(large_x1))
  597. y = fftn(small_x, shape=(4, 4), axes=(-1, -2))
  598. assert_array_almost_equal(y, swapaxes(
  599. fftn(swapaxes(large_x1, -1, -2)), -1, -2))
  600. def test_shape_axes_argument2(self):
  601. # Change shape of the last axis
  602. x = numpy.random.random((10, 5, 3, 7))
  603. y = fftn(x, axes=(-1,), shape=(8,))
  604. assert_array_almost_equal(y, fft(x, axis=-1, n=8))
  605. # Change shape of an arbitrary axis which is not the last one
  606. x = numpy.random.random((10, 5, 3, 7))
  607. y = fftn(x, axes=(-2,), shape=(8,))
  608. assert_array_almost_equal(y, fft(x, axis=-2, n=8))
  609. # Change shape of axes: cf #244, where shape and axes were mixed up
  610. x = numpy.random.random((4, 4, 2))
  611. y = fftn(x, axes=(-3, -2), shape=(8, 8))
  612. assert_array_almost_equal(y,
  613. numpy.fft.fftn(x, axes=(-3, -2), s=(8, 8)))
  614. def test_shape_argument_more(self):
  615. x = zeros((4, 4, 2))
  616. with assert_raises(ValueError,
  617. match="when given, axes and shape arguments"
  618. " have to be of the same length"):
  619. fftn(x, shape=(8, 8, 2, 1))
  620. def test_invalid_sizes(self):
  621. with assert_raises(ValueError,
  622. match="invalid number of data points"
  623. r" \(\[1 0\]\) specified"):
  624. fftn([[]])
  625. with assert_raises(ValueError,
  626. match="invalid number of data points"
  627. r" \(\[ 4 -3\]\) specified"):
  628. fftn([[1, 1], [2, 2]], (4, -3))
  629. class TestIfftn(object):
  630. dtype = None
  631. cdtype = None
  632. def setup_method(self):
  633. np.random.seed(1234)
  634. @pytest.mark.parametrize('dtype,cdtype,maxnlp',
  635. [(np.float64, np.complex128, 2000),
  636. (np.float32, np.complex64, 3500)])
  637. def test_definition(self, dtype, cdtype, maxnlp):
  638. x = np.array([[1, 2, 3],
  639. [4, 5, 6],
  640. [7, 8, 9]], dtype=dtype)
  641. y = ifftn(x)
  642. assert_equal(y.dtype, cdtype)
  643. assert_array_almost_equal_nulp(y, direct_idftn(x), maxnlp)
  644. x = random((20, 26))
  645. assert_array_almost_equal_nulp(ifftn(x), direct_idftn(x), maxnlp)
  646. x = random((5, 4, 3, 20))
  647. assert_array_almost_equal_nulp(ifftn(x), direct_idftn(x), maxnlp)
  648. @pytest.mark.parametrize('maxnlp', [2000, 3500])
  649. @pytest.mark.parametrize('size', [1, 2, 51, 32, 64, 92])
  650. def test_random_complex(self, maxnlp, size):
  651. x = random([size, size]) + 1j*random([size, size])
  652. assert_array_almost_equal_nulp(ifftn(fftn(x)), x, maxnlp)
  653. assert_array_almost_equal_nulp(fftn(ifftn(x)), x, maxnlp)
  654. def test_invalid_sizes(self):
  655. with assert_raises(ValueError,
  656. match="invalid number of data points"
  657. r" \(\[1 0\]\) specified"):
  658. ifftn([[]])
  659. with assert_raises(ValueError,
  660. match="invalid number of data points"
  661. r" \(\[ 4 -3\]\) specified"):
  662. ifftn([[1, 1], [2, 2]], (4, -3))
  663. class TestLongDoubleFailure(object):
  664. def setup_method(self):
  665. np.random.seed(1234)
  666. def test_complex(self):
  667. if np.dtype(np.longcomplex).itemsize == np.dtype(complex).itemsize:
  668. # longdouble == double; so fft is supported
  669. return
  670. x = np.random.randn(10).astype(np.longdouble) + \
  671. 1j * np.random.randn(10).astype(np.longdouble)
  672. for f in [fft, ifft]:
  673. try:
  674. f(x)
  675. raise AssertionError("Type {0} not supported but does not fail" %
  676. np.longcomplex)
  677. except ValueError:
  678. pass
  679. def test_real(self):
  680. if np.dtype(np.longdouble).itemsize == np.dtype(np.double).itemsize:
  681. # longdouble == double; so fft is supported
  682. return
  683. x = np.random.randn(10).astype(np.longcomplex)
  684. for f in [fft, ifft]:
  685. try:
  686. f(x)
  687. raise AssertionError("Type %r not supported but does not fail" %
  688. np.longcomplex)
  689. except ValueError:
  690. pass
  691. class FakeArray(object):
  692. def __init__(self, data):
  693. self._data = data
  694. self.__array_interface__ = data.__array_interface__
  695. class FakeArray2(object):
  696. def __init__(self, data):
  697. self._data = data
  698. def __array__(self):
  699. return self._data
  700. class TestOverwrite(object):
  701. """Check input overwrite behavior of the FFT functions."""
  702. real_dtypes = [np.float32, np.float64]
  703. dtypes = real_dtypes + [np.complex64, np.complex128]
  704. fftsizes = [8, 16, 32]
  705. def _check(self, x, routine, fftsize, axis, overwrite_x, should_overwrite):
  706. x2 = x.copy()
  707. for fake in [lambda x: x, FakeArray, FakeArray2]:
  708. routine(fake(x2), fftsize, axis, overwrite_x=overwrite_x)
  709. sig = "%s(%s%r, %r, axis=%r, overwrite_x=%r)" % (
  710. routine.__name__, x.dtype, x.shape, fftsize, axis, overwrite_x)
  711. if not should_overwrite:
  712. assert_equal(x2, x, err_msg="spurious overwrite in %s" % sig)
  713. def _check_1d(self, routine, dtype, shape, axis, overwritable_dtypes,
  714. fftsize, overwrite_x):
  715. np.random.seed(1234)
  716. if np.issubdtype(dtype, np.complexfloating):
  717. data = np.random.randn(*shape) + 1j*np.random.randn(*shape)
  718. else:
  719. data = np.random.randn(*shape)
  720. data = data.astype(dtype)
  721. should_overwrite = (overwrite_x
  722. and dtype in overwritable_dtypes
  723. and fftsize <= shape[axis]
  724. and (len(shape) == 1 or
  725. (axis % len(shape) == len(shape)-1
  726. and fftsize == shape[axis])))
  727. self._check(data, routine, fftsize, axis,
  728. overwrite_x=overwrite_x,
  729. should_overwrite=should_overwrite)
  730. @pytest.mark.parametrize('dtype', dtypes)
  731. @pytest.mark.parametrize('fftsize', fftsizes)
  732. @pytest.mark.parametrize('overwrite_x', [True, False])
  733. @pytest.mark.parametrize('shape,axes', [((16,), -1),
  734. ((16, 2), 0),
  735. ((2, 16), 1)])
  736. def test_fft_ifft(self, dtype, fftsize, overwrite_x, shape, axes):
  737. overwritable = (np.complex128, np.complex64)
  738. self._check_1d(fft, dtype, shape, axes, overwritable,
  739. fftsize, overwrite_x)
  740. self._check_1d(ifft, dtype, shape, axes, overwritable,
  741. fftsize, overwrite_x)
  742. @pytest.mark.parametrize('dtype', real_dtypes)
  743. @pytest.mark.parametrize('fftsize', fftsizes)
  744. @pytest.mark.parametrize('overwrite_x', [True, False])
  745. @pytest.mark.parametrize('shape,axes', [((16,), -1),
  746. ((16, 2), 0),
  747. ((2, 16), 1)])
  748. def test_rfft_irfft(self, dtype, fftsize, overwrite_x, shape, axes):
  749. overwritable = self.real_dtypes
  750. self._check_1d(irfft, dtype, shape, axes, overwritable,
  751. fftsize, overwrite_x)
  752. self._check_1d(rfft, dtype, shape, axes, overwritable,
  753. fftsize, overwrite_x)
  754. def _check_nd_one(self, routine, dtype, shape, axes, overwritable_dtypes,
  755. overwrite_x):
  756. np.random.seed(1234)
  757. if np.issubdtype(dtype, np.complexfloating):
  758. data = np.random.randn(*shape) + 1j*np.random.randn(*shape)
  759. else:
  760. data = np.random.randn(*shape)
  761. data = data.astype(dtype)
  762. def fftshape_iter(shp):
  763. if len(shp) <= 0:
  764. yield ()
  765. else:
  766. for j in (shp[0]//2, shp[0], shp[0]*2):
  767. for rest in fftshape_iter(shp[1:]):
  768. yield (j,) + rest
  769. if axes is None:
  770. part_shape = shape
  771. else:
  772. part_shape = tuple(np.take(shape, axes))
  773. for fftshape in fftshape_iter(part_shape):
  774. should_overwrite = (overwrite_x
  775. and data.ndim == 1
  776. and np.all([x < y for x, y in zip(fftshape,
  777. part_shape)])
  778. and dtype in overwritable_dtypes)
  779. self._check(data, routine, fftshape, axes,
  780. overwrite_x=overwrite_x,
  781. should_overwrite=should_overwrite)
  782. if data.ndim > 1:
  783. # check fortran order: it never overwrites
  784. self._check(data.T, routine, fftshape, axes,
  785. overwrite_x=overwrite_x,
  786. should_overwrite=False)
  787. @pytest.mark.parametrize('dtype', dtypes)
  788. @pytest.mark.parametrize('overwrite_x', [True, False])
  789. @pytest.mark.parametrize('shape,axes', [((16,), None),
  790. ((16,), (0,)),
  791. ((16, 2), (0,)),
  792. ((2, 16), (1,)),
  793. ((8, 16), None),
  794. ((8, 16), (0, 1)),
  795. ((8, 16, 2), (0, 1)),
  796. ((8, 16, 2), (1, 2)),
  797. ((8, 16, 2), (0,)),
  798. ((8, 16, 2), (1,)),
  799. ((8, 16, 2), (2,)),
  800. ((8, 16, 2), None),
  801. ((8, 16, 2), (0, 1, 2))])
  802. def test_fftn_ifftn(self, dtype, overwrite_x, shape, axes):
  803. overwritable = (np.complex128, np.complex64)
  804. self._check_nd_one(fftn, dtype, shape, axes, overwritable,
  805. overwrite_x)
  806. self._check_nd_one(ifftn, dtype, shape, axes, overwritable,
  807. overwrite_x)