first commit

This commit is contained in:
2018-07-24 15:40:59 +02:00
commit d0d03faed9
1884 changed files with 366061 additions and 0 deletions

View File

@@ -0,0 +1,11 @@
from __future__ import division, absolute_import, print_function
# To get sub-modules
from .info import __doc__
from .fftpack import *
from .helper import *
from numpy.testing import _numpy_tester
test = _numpy_tester().test
bench = _numpy_tester().bench

File diff suppressed because it is too large Load Diff

View File

@@ -0,0 +1,323 @@
"""
Discrete Fourier Transforms - helper.py
"""
from __future__ import division, absolute_import, print_function
import collections
import threading
from numpy.compat import integer_types
from numpy.core import (
asarray, concatenate, arange, take, integer, empty
)
# Created by Pearu Peterson, September 2002
__all__ = ['fftshift', 'ifftshift', 'fftfreq', 'rfftfreq']
integer_types = integer_types + (integer,)
def fftshift(x, axes=None):
"""
Shift the zero-frequency component to the center of the spectrum.
This function swaps half-spaces for all axes listed (defaults to all).
Note that ``y[0]`` is the Nyquist component only if ``len(x)`` is even.
Parameters
----------
x : array_like
Input array.
axes : int or shape tuple, optional
Axes over which to shift. Default is None, which shifts all axes.
Returns
-------
y : ndarray
The shifted array.
See Also
--------
ifftshift : The inverse of `fftshift`.
Examples
--------
>>> freqs = np.fft.fftfreq(10, 0.1)
>>> freqs
array([ 0., 1., 2., 3., 4., -5., -4., -3., -2., -1.])
>>> np.fft.fftshift(freqs)
array([-5., -4., -3., -2., -1., 0., 1., 2., 3., 4.])
Shift the zero-frequency component only along the second axis:
>>> freqs = np.fft.fftfreq(9, d=1./9).reshape(3, 3)
>>> freqs
array([[ 0., 1., 2.],
[ 3., 4., -4.],
[-3., -2., -1.]])
>>> np.fft.fftshift(freqs, axes=(1,))
array([[ 2., 0., 1.],
[-4., 3., 4.],
[-1., -3., -2.]])
"""
tmp = asarray(x)
ndim = tmp.ndim
if axes is None:
axes = list(range(ndim))
elif isinstance(axes, integer_types):
axes = (axes,)
y = tmp
for k in axes:
n = tmp.shape[k]
p2 = (n+1)//2
mylist = concatenate((arange(p2, n), arange(p2)))
y = take(y, mylist, k)
return y
def ifftshift(x, axes=None):
"""
The inverse of `fftshift`. Although identical for even-length `x`, the
functions differ by one sample for odd-length `x`.
Parameters
----------
x : array_like
Input array.
axes : int or shape tuple, optional
Axes over which to calculate. Defaults to None, which shifts all axes.
Returns
-------
y : ndarray
The shifted array.
See Also
--------
fftshift : Shift zero-frequency component to the center of the spectrum.
Examples
--------
>>> freqs = np.fft.fftfreq(9, d=1./9).reshape(3, 3)
>>> freqs
array([[ 0., 1., 2.],
[ 3., 4., -4.],
[-3., -2., -1.]])
>>> np.fft.ifftshift(np.fft.fftshift(freqs))
array([[ 0., 1., 2.],
[ 3., 4., -4.],
[-3., -2., -1.]])
"""
tmp = asarray(x)
ndim = tmp.ndim
if axes is None:
axes = list(range(ndim))
elif isinstance(axes, integer_types):
axes = (axes,)
y = tmp
for k in axes:
n = tmp.shape[k]
p2 = n-(n+1)//2
mylist = concatenate((arange(p2, n), arange(p2)))
y = take(y, mylist, k)
return y
def fftfreq(n, d=1.0):
"""
Return the Discrete Fourier Transform sample frequencies.
The returned float array `f` contains the frequency bin centers in cycles
per unit of the sample spacing (with zero at the start). For instance, if
the sample spacing is in seconds, then the frequency unit is cycles/second.
Given a window length `n` and a sample spacing `d`::
f = [0, 1, ..., n/2-1, -n/2, ..., -1] / (d*n) if n is even
f = [0, 1, ..., (n-1)/2, -(n-1)/2, ..., -1] / (d*n) if n is odd
Parameters
----------
n : int
Window length.
d : scalar, optional
Sample spacing (inverse of the sampling rate). Defaults to 1.
Returns
-------
f : ndarray
Array of length `n` containing the sample frequencies.
Examples
--------
>>> signal = np.array([-2, 8, 6, 4, 1, 0, 3, 5], dtype=float)
>>> fourier = np.fft.fft(signal)
>>> n = signal.size
>>> timestep = 0.1
>>> freq = np.fft.fftfreq(n, d=timestep)
>>> freq
array([ 0. , 1.25, 2.5 , 3.75, -5. , -3.75, -2.5 , -1.25])
"""
if not isinstance(n, integer_types):
raise ValueError("n should be an integer")
val = 1.0 / (n * d)
results = empty(n, int)
N = (n-1)//2 + 1
p1 = arange(0, N, dtype=int)
results[:N] = p1
p2 = arange(-(n//2), 0, dtype=int)
results[N:] = p2
return results * val
#return hstack((arange(0,(n-1)/2 + 1), arange(-(n/2),0))) / (n*d)
def rfftfreq(n, d=1.0):
"""
Return the Discrete Fourier Transform sample frequencies
(for usage with rfft, irfft).
The returned float array `f` contains the frequency bin centers in cycles
per unit of the sample spacing (with zero at the start). For instance, if
the sample spacing is in seconds, then the frequency unit is cycles/second.
Given a window length `n` and a sample spacing `d`::
f = [0, 1, ..., n/2-1, n/2] / (d*n) if n is even
f = [0, 1, ..., (n-1)/2-1, (n-1)/2] / (d*n) if n is odd
Unlike `fftfreq` (but like `scipy.fftpack.rfftfreq`)
the Nyquist frequency component is considered to be positive.
Parameters
----------
n : int
Window length.
d : scalar, optional
Sample spacing (inverse of the sampling rate). Defaults to 1.
Returns
-------
f : ndarray
Array of length ``n//2 + 1`` containing the sample frequencies.
Examples
--------
>>> signal = np.array([-2, 8, 6, 4, 1, 0, 3, 5, -3, 4], dtype=float)
>>> fourier = np.fft.rfft(signal)
>>> n = signal.size
>>> sample_rate = 100
>>> freq = np.fft.fftfreq(n, d=1./sample_rate)
>>> freq
array([ 0., 10., 20., 30., 40., -50., -40., -30., -20., -10.])
>>> freq = np.fft.rfftfreq(n, d=1./sample_rate)
>>> freq
array([ 0., 10., 20., 30., 40., 50.])
"""
if not isinstance(n, integer_types):
raise ValueError("n should be an integer")
val = 1.0/(n*d)
N = n//2 + 1
results = arange(0, N, dtype=int)
return results * val
class _FFTCache(object):
"""
Cache for the FFT twiddle factors as an LRU (least recently used) cache.
Parameters
----------
max_size_in_mb : int
Maximum memory usage of the cache before items are being evicted.
max_item_count : int
Maximum item count of the cache before items are being evicted.
Notes
-----
Items will be evicted if either limit has been reached upon getting and
setting. The maximum memory usages is not strictly the given
``max_size_in_mb`` but rather
``max(max_size_in_mb, 1.5 * size_of_largest_item)``. Thus the cache will
never be completely cleared - at least one item will remain and a single
large item can cause the cache to retain several smaller items even if the
given maximum cache size has been exceeded.
"""
def __init__(self, max_size_in_mb, max_item_count):
self._max_size_in_bytes = max_size_in_mb * 1024 ** 2
self._max_item_count = max_item_count
self._dict = collections.OrderedDict()
self._lock = threading.Lock()
def put_twiddle_factors(self, n, factors):
"""
Store twiddle factors for an FFT of length n in the cache.
Putting multiple twiddle factors for a certain n will store it multiple
times.
Parameters
----------
n : int
Data length for the FFT.
factors : ndarray
The actual twiddle values.
"""
with self._lock:
# Pop + later add to move it to the end for LRU behavior.
# Internally everything is stored in a dictionary whose values are
# lists.
try:
value = self._dict.pop(n)
except KeyError:
value = []
value.append(factors)
self._dict[n] = value
self._prune_cache()
def pop_twiddle_factors(self, n):
"""
Pop twiddle factors for an FFT of length n from the cache.
Will return None if the requested twiddle factors are not available in
the cache.
Parameters
----------
n : int
Data length for the FFT.
Returns
-------
out : ndarray or None
The retrieved twiddle factors if available, else None.
"""
with self._lock:
if n not in self._dict or not self._dict[n]:
return None
# Pop + later add to move it to the end for LRU behavior.
all_values = self._dict.pop(n)
value = all_values.pop()
# Only put pack if there are still some arrays left in the list.
if all_values:
self._dict[n] = all_values
return value
def _prune_cache(self):
# Always keep at least one item.
while len(self._dict) > 1 and (
len(self._dict) > self._max_item_count or self._check_size()):
self._dict.popitem(last=False)
def _check_size(self):
item_sizes = [sum(_j.nbytes for _j in _i)
for _i in self._dict.values() if _i]
if not item_sizes:
return False
max_size = max(self._max_size_in_bytes, 1.5 * max(item_sizes))
return sum(item_sizes) > max_size

View File

@@ -0,0 +1,187 @@
"""
Discrete Fourier Transform (:mod:`numpy.fft`)
=============================================
.. currentmodule:: numpy.fft
Standard FFTs
-------------
.. autosummary::
:toctree: generated/
fft Discrete Fourier transform.
ifft Inverse discrete Fourier transform.
fft2 Discrete Fourier transform in two dimensions.
ifft2 Inverse discrete Fourier transform in two dimensions.
fftn Discrete Fourier transform in N-dimensions.
ifftn Inverse discrete Fourier transform in N dimensions.
Real FFTs
---------
.. autosummary::
:toctree: generated/
rfft Real discrete Fourier transform.
irfft Inverse real discrete Fourier transform.
rfft2 Real discrete Fourier transform in two dimensions.
irfft2 Inverse real discrete Fourier transform in two dimensions.
rfftn Real discrete Fourier transform in N dimensions.
irfftn Inverse real discrete Fourier transform in N dimensions.
Hermitian FFTs
--------------
.. autosummary::
:toctree: generated/
hfft Hermitian discrete Fourier transform.
ihfft Inverse Hermitian discrete Fourier transform.
Helper routines
---------------
.. autosummary::
:toctree: generated/
fftfreq Discrete Fourier Transform sample frequencies.
rfftfreq DFT sample frequencies (for usage with rfft, irfft).
fftshift Shift zero-frequency component to center of spectrum.
ifftshift Inverse of fftshift.
Background information
----------------------
Fourier analysis is fundamentally a method for expressing a function as a
sum of periodic components, and for recovering the function from those
components. When both the function and its Fourier transform are
replaced with discretized counterparts, it is called the discrete Fourier
transform (DFT). The DFT has become a mainstay of numerical computing in
part because of a very fast algorithm for computing it, called the Fast
Fourier Transform (FFT), which was known to Gauss (1805) and was brought
to light in its current form by Cooley and Tukey [CT]_. Press et al. [NR]_
provide an accessible introduction to Fourier analysis and its
applications.
Because the discrete Fourier transform separates its input into
components that contribute at discrete frequencies, it has a great number
of applications in digital signal processing, e.g., for filtering, and in
this context the discretized input to the transform is customarily
referred to as a *signal*, which exists in the *time domain*. The output
is called a *spectrum* or *transform* and exists in the *frequency
domain*.
Implementation details
----------------------
There are many ways to define the DFT, varying in the sign of the
exponent, normalization, etc. In this implementation, the DFT is defined
as
.. math::
A_k = \\sum_{m=0}^{n-1} a_m \\exp\\left\\{-2\\pi i{mk \\over n}\\right\\}
\\qquad k = 0,\\ldots,n-1.
The DFT is in general defined for complex inputs and outputs, and a
single-frequency component at linear frequency :math:`f` is
represented by a complex exponential
:math:`a_m = \\exp\\{2\\pi i\\,f m\\Delta t\\}`, where :math:`\\Delta t`
is the sampling interval.
The values in the result follow so-called "standard" order: If ``A =
fft(a, n)``, then ``A[0]`` contains the zero-frequency term (the sum of
the signal), which is always purely real for real inputs. Then ``A[1:n/2]``
contains the positive-frequency terms, and ``A[n/2+1:]`` contains the
negative-frequency terms, in order of decreasingly negative frequency.
For an even number of input points, ``A[n/2]`` represents both positive and
negative Nyquist frequency, and is also purely real for real input. For
an odd number of input points, ``A[(n-1)/2]`` contains the largest positive
frequency, while ``A[(n+1)/2]`` contains the largest negative frequency.
The routine ``np.fft.fftfreq(n)`` returns an array giving the frequencies
of corresponding elements in the output. The routine
``np.fft.fftshift(A)`` shifts transforms and their frequencies to put the
zero-frequency components in the middle, and ``np.fft.ifftshift(A)`` undoes
that shift.
When the input `a` is a time-domain signal and ``A = fft(a)``, ``np.abs(A)``
is its amplitude spectrum and ``np.abs(A)**2`` is its power spectrum.
The phase spectrum is obtained by ``np.angle(A)``.
The inverse DFT is defined as
.. math::
a_m = \\frac{1}{n}\\sum_{k=0}^{n-1}A_k\\exp\\left\\{2\\pi i{mk\\over n}\\right\\}
\\qquad m = 0,\\ldots,n-1.
It differs from the forward transform by the sign of the exponential
argument and the default normalization by :math:`1/n`.
Normalization
-------------
The default normalization has the direct transforms unscaled and the inverse
transforms are scaled by :math:`1/n`. It is possible to obtain unitary
transforms by setting the keyword argument ``norm`` to ``"ortho"`` (default is
`None`) so that both direct and inverse transforms will be scaled by
:math:`1/\\sqrt{n}`.
Real and Hermitian transforms
-----------------------------
When the input is purely real, its transform is Hermitian, i.e., the
component at frequency :math:`f_k` is the complex conjugate of the
component at frequency :math:`-f_k`, which means that for real
inputs there is no information in the negative frequency components that
is not already available from the positive frequency components.
The family of `rfft` functions is
designed to operate on real inputs, and exploits this symmetry by
computing only the positive frequency components, up to and including the
Nyquist frequency. Thus, ``n`` input points produce ``n/2+1`` complex
output points. The inverses of this family assumes the same symmetry of
its input, and for an output of ``n`` points uses ``n/2+1`` input points.
Correspondingly, when the spectrum is purely real, the signal is
Hermitian. The `hfft` family of functions exploits this symmetry by
using ``n/2+1`` complex points in the input (time) domain for ``n`` real
points in the frequency domain.
In higher dimensions, FFTs are used, e.g., for image analysis and
filtering. The computational efficiency of the FFT means that it can
also be a faster way to compute large convolutions, using the property
that a convolution in the time domain is equivalent to a point-by-point
multiplication in the frequency domain.
Higher dimensions
-----------------
In two dimensions, the DFT is defined as
.. math::
A_{kl} = \\sum_{m=0}^{M-1} \\sum_{n=0}^{N-1}
a_{mn}\\exp\\left\\{-2\\pi i \\left({mk\\over M}+{nl\\over N}\\right)\\right\\}
\\qquad k = 0, \\ldots, M-1;\\quad l = 0, \\ldots, N-1,
which extends in the obvious way to higher dimensions, and the inverses
in higher dimensions also extend in the same way.
References
----------
.. [CT] Cooley, James W., and John W. Tukey, 1965, "An algorithm for the
machine calculation of complex Fourier series," *Math. Comput.*
19: 297-301.
.. [NR] Press, W., Teukolsky, S., Vetterline, W.T., and Flannery, B.P.,
2007, *Numerical Recipes: The Art of Scientific Computing*, ch.
12-13. Cambridge Univ. Press, Cambridge, UK.
Examples
--------
For examples, see the various functions.
"""
from __future__ import division, absolute_import, print_function
depends = ['core']

View File

@@ -0,0 +1,19 @@
from __future__ import division, print_function
def configuration(parent_package='',top_path=None):
from numpy.distutils.misc_util import Configuration
config = Configuration('fft', parent_package, top_path)
config.add_data_dir('tests')
# Configure fftpack_lite
config.add_extension('fftpack_lite',
sources=['fftpack_litemodule.c', 'fftpack.c']
)
return config
if __name__ == '__main__':
from numpy.distutils.core import setup
setup(configuration=configuration)

View File

@@ -0,0 +1,190 @@
from __future__ import division, absolute_import, print_function
import numpy as np
from numpy.random import random
from numpy.testing import (
run_module_suite, assert_array_almost_equal, assert_array_equal,
assert_raises,
)
import threading
import sys
if sys.version_info[0] >= 3:
import queue
else:
import Queue as queue
def fft1(x):
L = len(x)
phase = -2j*np.pi*(np.arange(L)/float(L))
phase = np.arange(L).reshape(-1, 1) * phase
return np.sum(x*np.exp(phase), axis=1)
class TestFFTShift(object):
def test_fft_n(self):
assert_raises(ValueError, np.fft.fft, [1, 2, 3], 0)
class TestFFT1D(object):
def test_fft(self):
x = random(30) + 1j*random(30)
assert_array_almost_equal(fft1(x), np.fft.fft(x))
assert_array_almost_equal(fft1(x) / np.sqrt(30),
np.fft.fft(x, norm="ortho"))
def test_ifft(self):
x = random(30) + 1j*random(30)
assert_array_almost_equal(x, np.fft.ifft(np.fft.fft(x)))
assert_array_almost_equal(
x, np.fft.ifft(np.fft.fft(x, norm="ortho"), norm="ortho"))
def test_fft2(self):
x = random((30, 20)) + 1j*random((30, 20))
assert_array_almost_equal(np.fft.fft(np.fft.fft(x, axis=1), axis=0),
np.fft.fft2(x))
assert_array_almost_equal(np.fft.fft2(x) / np.sqrt(30 * 20),
np.fft.fft2(x, norm="ortho"))
def test_ifft2(self):
x = random((30, 20)) + 1j*random((30, 20))
assert_array_almost_equal(np.fft.ifft(np.fft.ifft(x, axis=1), axis=0),
np.fft.ifft2(x))
assert_array_almost_equal(np.fft.ifft2(x) * np.sqrt(30 * 20),
np.fft.ifft2(x, norm="ortho"))
def test_fftn(self):
x = random((30, 20, 10)) + 1j*random((30, 20, 10))
assert_array_almost_equal(
np.fft.fft(np.fft.fft(np.fft.fft(x, axis=2), axis=1), axis=0),
np.fft.fftn(x))
assert_array_almost_equal(np.fft.fftn(x) / np.sqrt(30 * 20 * 10),
np.fft.fftn(x, norm="ortho"))
def test_ifftn(self):
x = random((30, 20, 10)) + 1j*random((30, 20, 10))
assert_array_almost_equal(
np.fft.ifft(np.fft.ifft(np.fft.ifft(x, axis=2), axis=1), axis=0),
np.fft.ifftn(x))
assert_array_almost_equal(np.fft.ifftn(x) * np.sqrt(30 * 20 * 10),
np.fft.ifftn(x, norm="ortho"))
def test_rfft(self):
x = random(30)
for n in [x.size, 2*x.size]:
for norm in [None, 'ortho']:
assert_array_almost_equal(
np.fft.fft(x, n=n, norm=norm)[:(n//2 + 1)],
np.fft.rfft(x, n=n, norm=norm))
assert_array_almost_equal(np.fft.rfft(x, n=n) / np.sqrt(n),
np.fft.rfft(x, n=n, norm="ortho"))
def test_irfft(self):
x = random(30)
assert_array_almost_equal(x, np.fft.irfft(np.fft.rfft(x)))
assert_array_almost_equal(
x, np.fft.irfft(np.fft.rfft(x, norm="ortho"), norm="ortho"))
def test_rfft2(self):
x = random((30, 20))
assert_array_almost_equal(np.fft.fft2(x)[:, :11], np.fft.rfft2(x))
assert_array_almost_equal(np.fft.rfft2(x) / np.sqrt(30 * 20),
np.fft.rfft2(x, norm="ortho"))
def test_irfft2(self):
x = random((30, 20))
assert_array_almost_equal(x, np.fft.irfft2(np.fft.rfft2(x)))
assert_array_almost_equal(
x, np.fft.irfft2(np.fft.rfft2(x, norm="ortho"), norm="ortho"))
def test_rfftn(self):
x = random((30, 20, 10))
assert_array_almost_equal(np.fft.fftn(x)[:, :, :6], np.fft.rfftn(x))
assert_array_almost_equal(np.fft.rfftn(x) / np.sqrt(30 * 20 * 10),
np.fft.rfftn(x, norm="ortho"))
def test_irfftn(self):
x = random((30, 20, 10))
assert_array_almost_equal(x, np.fft.irfftn(np.fft.rfftn(x)))
assert_array_almost_equal(
x, np.fft.irfftn(np.fft.rfftn(x, norm="ortho"), norm="ortho"))
def test_hfft(self):
x = random(14) + 1j*random(14)
x_herm = np.concatenate((random(1), x, random(1)))
x = np.concatenate((x_herm, x[::-1].conj()))
assert_array_almost_equal(np.fft.fft(x), np.fft.hfft(x_herm))
assert_array_almost_equal(np.fft.hfft(x_herm) / np.sqrt(30),
np.fft.hfft(x_herm, norm="ortho"))
def test_ihttf(self):
x = random(14) + 1j*random(14)
x_herm = np.concatenate((random(1), x, random(1)))
x = np.concatenate((x_herm, x[::-1].conj()))
assert_array_almost_equal(x_herm, np.fft.ihfft(np.fft.hfft(x_herm)))
assert_array_almost_equal(
x_herm, np.fft.ihfft(np.fft.hfft(x_herm, norm="ortho"),
norm="ortho"))
def test_all_1d_norm_preserving(self):
# verify that round-trip transforms are norm-preserving
x = random(30)
x_norm = np.linalg.norm(x)
n = x.size * 2
func_pairs = [(np.fft.fft, np.fft.ifft),
(np.fft.rfft, np.fft.irfft),
# hfft: order so the first function takes x.size samples
# (necessary for comparison to x_norm above)
(np.fft.ihfft, np.fft.hfft),
]
for forw, back in func_pairs:
for n in [x.size, 2*x.size]:
for norm in [None, 'ortho']:
tmp = forw(x, n=n, norm=norm)
tmp = back(tmp, n=n, norm=norm)
assert_array_almost_equal(x_norm,
np.linalg.norm(tmp))
class TestFFTThreadSafe(object):
threads = 16
input_shape = (800, 200)
def _test_mtsame(self, func, *args):
def worker(args, q):
q.put(func(*args))
q = queue.Queue()
expected = func(*args)
# Spin off a bunch of threads to call the same function simultaneously
t = [threading.Thread(target=worker, args=(args, q))
for i in range(self.threads)]
[x.start() for x in t]
[x.join() for x in t]
# Make sure all threads returned the correct value
for i in range(self.threads):
assert_array_equal(q.get(timeout=5), expected,
'Function returned wrong value in multithreaded context')
def test_fft(self):
a = np.ones(self.input_shape) * 1+0j
self._test_mtsame(np.fft.fft, a)
def test_ifft(self):
a = np.ones(self.input_shape) * 1+0j
self._test_mtsame(np.fft.ifft, a)
def test_rfft(self):
a = np.ones(self.input_shape)
self._test_mtsame(np.fft.rfft, a)
def test_irfft(self):
a = np.ones(self.input_shape) * 1+0j
self._test_mtsame(np.fft.irfft, a)
if __name__ == "__main__":
run_module_suite()

View File

@@ -0,0 +1,158 @@
"""Test functions for fftpack.helper module
Copied from fftpack.helper by Pearu Peterson, October 2005
"""
from __future__ import division, absolute_import, print_function
import numpy as np
from numpy.testing import (
run_module_suite, assert_array_almost_equal, assert_equal,
)
from numpy import fft
from numpy import pi
from numpy.fft.helper import _FFTCache
class TestFFTShift(object):
def test_definition(self):
x = [0, 1, 2, 3, 4, -4, -3, -2, -1]
y = [-4, -3, -2, -1, 0, 1, 2, 3, 4]
assert_array_almost_equal(fft.fftshift(x), y)
assert_array_almost_equal(fft.ifftshift(y), x)
x = [0, 1, 2, 3, 4, -5, -4, -3, -2, -1]
y = [-5, -4, -3, -2, -1, 0, 1, 2, 3, 4]
assert_array_almost_equal(fft.fftshift(x), y)
assert_array_almost_equal(fft.ifftshift(y), x)
def test_inverse(self):
for n in [1, 4, 9, 100, 211]:
x = np.random.random((n,))
assert_array_almost_equal(fft.ifftshift(fft.fftshift(x)), x)
def test_axes_keyword(self):
freqs = [[0, 1, 2], [3, 4, -4], [-3, -2, -1]]
shifted = [[-1, -3, -2], [2, 0, 1], [-4, 3, 4]]
assert_array_almost_equal(fft.fftshift(freqs, axes=(0, 1)), shifted)
assert_array_almost_equal(fft.fftshift(freqs, axes=0),
fft.fftshift(freqs, axes=(0,)))
assert_array_almost_equal(fft.ifftshift(shifted, axes=(0, 1)), freqs)
assert_array_almost_equal(fft.ifftshift(shifted, axes=0),
fft.ifftshift(shifted, axes=(0,)))
class TestFFTFreq(object):
def test_definition(self):
x = [0, 1, 2, 3, 4, -4, -3, -2, -1]
assert_array_almost_equal(9*fft.fftfreq(9), x)
assert_array_almost_equal(9*pi*fft.fftfreq(9, pi), x)
x = [0, 1, 2, 3, 4, -5, -4, -3, -2, -1]
assert_array_almost_equal(10*fft.fftfreq(10), x)
assert_array_almost_equal(10*pi*fft.fftfreq(10, pi), x)
class TestRFFTFreq(object):
def test_definition(self):
x = [0, 1, 2, 3, 4]
assert_array_almost_equal(9*fft.rfftfreq(9), x)
assert_array_almost_equal(9*pi*fft.rfftfreq(9, pi), x)
x = [0, 1, 2, 3, 4, 5]
assert_array_almost_equal(10*fft.rfftfreq(10), x)
assert_array_almost_equal(10*pi*fft.rfftfreq(10, pi), x)
class TestIRFFTN(object):
def test_not_last_axis_success(self):
ar, ai = np.random.random((2, 16, 8, 32))
a = ar + 1j*ai
axes = (-2,)
# Should not raise error
fft.irfftn(a, axes=axes)
class TestFFTCache(object):
def test_basic_behaviour(self):
c = _FFTCache(max_size_in_mb=1, max_item_count=4)
# Put
c.put_twiddle_factors(1, np.ones(2, dtype=np.float32))
c.put_twiddle_factors(2, np.zeros(2, dtype=np.float32))
# Get
assert_array_almost_equal(c.pop_twiddle_factors(1),
np.ones(2, dtype=np.float32))
assert_array_almost_equal(c.pop_twiddle_factors(2),
np.zeros(2, dtype=np.float32))
# Nothing should be left.
assert_equal(len(c._dict), 0)
# Now put everything in twice so it can be retrieved once and each will
# still have one item left.
for _ in range(2):
c.put_twiddle_factors(1, np.ones(2, dtype=np.float32))
c.put_twiddle_factors(2, np.zeros(2, dtype=np.float32))
assert_array_almost_equal(c.pop_twiddle_factors(1),
np.ones(2, dtype=np.float32))
assert_array_almost_equal(c.pop_twiddle_factors(2),
np.zeros(2, dtype=np.float32))
assert_equal(len(c._dict), 2)
def test_automatic_pruning(self):
# That's around 2600 single precision samples.
c = _FFTCache(max_size_in_mb=0.01, max_item_count=4)
c.put_twiddle_factors(1, np.ones(200, dtype=np.float32))
c.put_twiddle_factors(2, np.ones(200, dtype=np.float32))
assert_equal(list(c._dict.keys()), [1, 2])
# This is larger than the limit but should still be kept.
c.put_twiddle_factors(3, np.ones(3000, dtype=np.float32))
assert_equal(list(c._dict.keys()), [1, 2, 3])
# Add one more.
c.put_twiddle_factors(4, np.ones(3000, dtype=np.float32))
# The other three should no longer exist.
assert_equal(list(c._dict.keys()), [4])
# Now test the max item count pruning.
c = _FFTCache(max_size_in_mb=0.01, max_item_count=2)
c.put_twiddle_factors(2, np.empty(2))
c.put_twiddle_factors(1, np.empty(2))
# Can still be accessed.
assert_equal(list(c._dict.keys()), [2, 1])
c.put_twiddle_factors(3, np.empty(2))
# 1 and 3 can still be accessed - c[2] has been touched least recently
# and is thus evicted.
assert_equal(list(c._dict.keys()), [1, 3])
# One last test. We will add a single large item that is slightly
# bigger then the cache size. Some small items can still be added.
c = _FFTCache(max_size_in_mb=0.01, max_item_count=5)
c.put_twiddle_factors(1, np.ones(3000, dtype=np.float32))
c.put_twiddle_factors(2, np.ones(2, dtype=np.float32))
c.put_twiddle_factors(3, np.ones(2, dtype=np.float32))
c.put_twiddle_factors(4, np.ones(2, dtype=np.float32))
assert_equal(list(c._dict.keys()), [1, 2, 3, 4])
# One more big item. This time it is 6 smaller ones but they are
# counted as one big item.
for _ in range(6):
c.put_twiddle_factors(5, np.ones(500, dtype=np.float32))
# '1' no longer in the cache. Rest still in the cache.
assert_equal(list(c._dict.keys()), [2, 3, 4, 5])
# Another big item - should now be the only item in the cache.
c.put_twiddle_factors(6, np.ones(4000, dtype=np.float32))
assert_equal(list(c._dict.keys()), [6])
if __name__ == "__main__":
run_module_suite()