Initial commit
This commit is contained in:
Untriex Programming
@@ -0,0 +1,212 @@
|
||||
"""
|
||||
Discrete Fourier Transform (:mod:`numpy.fft`)
|
||||
=============================================
|
||||
|
||||
.. currentmodule:: numpy.fft
|
||||
|
||||
The SciPy module `scipy.fft` is a more comprehensive superset
|
||||
of ``numpy.fft``, which includes only a basic set of routines.
|
||||
|
||||
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`.
|
||||
|
||||
Type Promotion
|
||||
--------------
|
||||
|
||||
`numpy.fft` promotes ``float32`` and ``complex64`` arrays to ``float64`` and
|
||||
``complex128`` arrays respectively. For an FFT implementation that does not
|
||||
promote input arrays, see `scipy.fftpack`.
|
||||
|
||||
Normalization
|
||||
-------------
|
||||
|
||||
The argument ``norm`` indicates which direction of the pair of direct/inverse
|
||||
transforms is scaled and with what normalization factor.
|
||||
The default normalization (``"backward"``) has the direct (forward) transforms
|
||||
unscaled and the inverse (backward) transforms scaled by :math:`1/n`. It is
|
||||
possible to obtain unitary transforms by setting the keyword argument ``norm``
|
||||
to ``"ortho"`` so that both direct and inverse transforms are scaled by
|
||||
:math:`1/\\sqrt{n}`. Finally, setting the keyword argument ``norm`` to
|
||||
``"forward"`` has the direct transforms scaled by :math:`1/n` and the inverse
|
||||
transforms unscaled (i.e. exactly opposite to the default ``"backward"``).
|
||||
`None` is an alias of the default option ``"backward"`` for backward
|
||||
compatibility.
|
||||
|
||||
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 . import _pocketfft, helper
|
||||
from ._pocketfft import *
|
||||
from .helper import *
|
||||
|
||||
__all__ = _pocketfft.__all__.copy()
|
||||
__all__ += helper.__all__
|
||||
|
||||
from numpy._pytesttester import PytestTester
|
||||
test = PytestTester(__name__)
|
||||
del PytestTester
|
||||
@@ -0,0 +1,22 @@
|
||||
from typing import Any, List
|
||||
|
||||
__all__: List[str]
|
||||
|
||||
def fft(a, n=..., axis=..., norm=...): ...
|
||||
def ifft(a, n=..., axis=..., norm=...): ...
|
||||
def rfft(a, n=..., axis=..., norm=...): ...
|
||||
def irfft(a, n=..., axis=..., norm=...): ...
|
||||
def hfft(a, n=..., axis=..., norm=...): ...
|
||||
def ihfft(a, n=..., axis=..., norm=...): ...
|
||||
def fftn(a, s=..., axes=..., norm=...): ...
|
||||
def ifftn(a, s=..., axes=..., norm=...): ...
|
||||
def rfftn(a, s=..., axes=..., norm=...): ...
|
||||
def irfftn(a, s=..., axes=..., norm=...): ...
|
||||
def fft2(a, s=..., axes=..., norm=...): ...
|
||||
def ifft2(a, s=..., axes=..., norm=...): ...
|
||||
def rfft2(a, s=..., axes=..., norm=...): ...
|
||||
def irfft2(a, s=..., axes=..., norm=...): ...
|
||||
def fftshift(x, axes=...): ...
|
||||
def ifftshift(x, axes=...): ...
|
||||
def fftfreq(n, d=...): ...
|
||||
def rfftfreq(n, d=...): ...
|
||||
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
File diff suppressed because it is too large
Load Diff
Executable
BIN
Binary file not shown.
@@ -0,0 +1,221 @@
|
||||
"""
|
||||
Discrete Fourier Transforms - helper.py
|
||||
|
||||
"""
|
||||
from numpy.core import integer, empty, arange, asarray, roll
|
||||
from numpy.core.overrides import array_function_dispatch, set_module
|
||||
|
||||
# Created by Pearu Peterson, September 2002
|
||||
|
||||
__all__ = ['fftshift', 'ifftshift', 'fftfreq', 'rfftfreq']
|
||||
|
||||
integer_types = (int, integer)
|
||||
|
||||
|
||||
def _fftshift_dispatcher(x, axes=None):
|
||||
return (x,)
|
||||
|
||||
|
||||
@array_function_dispatch(_fftshift_dispatcher, module='numpy.fft')
|
||||
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., -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.]])
|
||||
|
||||
"""
|
||||
x = asarray(x)
|
||||
if axes is None:
|
||||
axes = tuple(range(x.ndim))
|
||||
shift = [dim // 2 for dim in x.shape]
|
||||
elif isinstance(axes, integer_types):
|
||||
shift = x.shape[axes] // 2
|
||||
else:
|
||||
shift = [x.shape[ax] // 2 for ax in axes]
|
||||
|
||||
return roll(x, shift, axes)
|
||||
|
||||
|
||||
@array_function_dispatch(_fftshift_dispatcher, module='numpy.fft')
|
||||
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.]])
|
||||
|
||||
"""
|
||||
x = asarray(x)
|
||||
if axes is None:
|
||||
axes = tuple(range(x.ndim))
|
||||
shift = [-(dim // 2) for dim in x.shape]
|
||||
elif isinstance(axes, integer_types):
|
||||
shift = -(x.shape[axes] // 2)
|
||||
else:
|
||||
shift = [-(x.shape[ax] // 2) for ax in axes]
|
||||
|
||||
return roll(x, shift, axes)
|
||||
|
||||
|
||||
@set_module('numpy.fft')
|
||||
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, -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
|
||||
|
||||
|
||||
@set_module('numpy.fft')
|
||||
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., -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
|
||||
@@ -0,0 +1,22 @@
|
||||
import sys
|
||||
|
||||
def configuration(parent_package='',top_path=None):
|
||||
from numpy.distutils.misc_util import Configuration
|
||||
config = Configuration('fft', parent_package, top_path)
|
||||
|
||||
config.add_subpackage('tests')
|
||||
|
||||
# AIX needs to be told to use large file support - at all times
|
||||
defs = [('_LARGE_FILES', None)] if sys.platform[:3] == "aix" else []
|
||||
# Configure pocketfft_internal
|
||||
config.add_extension('_pocketfft_internal',
|
||||
sources=['_pocketfft.c'],
|
||||
define_macros=defs,
|
||||
)
|
||||
|
||||
config.add_data_files('*.pyi')
|
||||
return config
|
||||
|
||||
if __name__ == '__main__':
|
||||
from numpy.distutils.core import setup
|
||||
setup(configuration=configuration)
|
||||
Binary file not shown.
BIN
Binary file not shown.
BIN
Binary file not shown.
@@ -0,0 +1,167 @@
|
||||
"""Test functions for fftpack.helper module
|
||||
|
||||
Copied from fftpack.helper by Pearu Peterson, October 2005
|
||||
|
||||
"""
|
||||
import numpy as np
|
||||
from numpy.testing import assert_array_almost_equal
|
||||
from numpy import fft, pi
|
||||
|
||||
|
||||
class TestFFTShift:
|
||||
|
||||
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,)))
|
||||
|
||||
assert_array_almost_equal(fft.fftshift(freqs), shifted)
|
||||
assert_array_almost_equal(fft.ifftshift(shifted), freqs)
|
||||
|
||||
def test_uneven_dims(self):
|
||||
""" Test 2D input, which has uneven dimension sizes """
|
||||
freqs = [
|
||||
[0, 1],
|
||||
[2, 3],
|
||||
[4, 5]
|
||||
]
|
||||
|
||||
# shift in dimension 0
|
||||
shift_dim0 = [
|
||||
[4, 5],
|
||||
[0, 1],
|
||||
[2, 3]
|
||||
]
|
||||
assert_array_almost_equal(fft.fftshift(freqs, axes=0), shift_dim0)
|
||||
assert_array_almost_equal(fft.ifftshift(shift_dim0, axes=0), freqs)
|
||||
assert_array_almost_equal(fft.fftshift(freqs, axes=(0,)), shift_dim0)
|
||||
assert_array_almost_equal(fft.ifftshift(shift_dim0, axes=[0]), freqs)
|
||||
|
||||
# shift in dimension 1
|
||||
shift_dim1 = [
|
||||
[1, 0],
|
||||
[3, 2],
|
||||
[5, 4]
|
||||
]
|
||||
assert_array_almost_equal(fft.fftshift(freqs, axes=1), shift_dim1)
|
||||
assert_array_almost_equal(fft.ifftshift(shift_dim1, axes=1), freqs)
|
||||
|
||||
# shift in both dimensions
|
||||
shift_dim_both = [
|
||||
[5, 4],
|
||||
[1, 0],
|
||||
[3, 2]
|
||||
]
|
||||
assert_array_almost_equal(fft.fftshift(freqs, axes=(0, 1)), shift_dim_both)
|
||||
assert_array_almost_equal(fft.ifftshift(shift_dim_both, axes=(0, 1)), freqs)
|
||||
assert_array_almost_equal(fft.fftshift(freqs, axes=[0, 1]), shift_dim_both)
|
||||
assert_array_almost_equal(fft.ifftshift(shift_dim_both, axes=[0, 1]), freqs)
|
||||
|
||||
# axes=None (default) shift in all dimensions
|
||||
assert_array_almost_equal(fft.fftshift(freqs, axes=None), shift_dim_both)
|
||||
assert_array_almost_equal(fft.ifftshift(shift_dim_both, axes=None), freqs)
|
||||
assert_array_almost_equal(fft.fftshift(freqs), shift_dim_both)
|
||||
assert_array_almost_equal(fft.ifftshift(shift_dim_both), freqs)
|
||||
|
||||
def test_equal_to_original(self):
|
||||
""" Test that the new (>=v1.15) implementation (see #10073) is equal to the original (<=v1.14) """
|
||||
from numpy.core import asarray, concatenate, arange, take
|
||||
|
||||
def original_fftshift(x, axes=None):
|
||||
""" How fftshift was implemented in v1.14"""
|
||||
tmp = asarray(x)
|
||||
ndim = tmp.ndim
|
||||
if axes is None:
|
||||
axes = list(range(ndim))
|
||||
elif isinstance(axes, int):
|
||||
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 original_ifftshift(x, axes=None):
|
||||
""" How ifftshift was implemented in v1.14 """
|
||||
tmp = asarray(x)
|
||||
ndim = tmp.ndim
|
||||
if axes is None:
|
||||
axes = list(range(ndim))
|
||||
elif isinstance(axes, int):
|
||||
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
|
||||
|
||||
# create possible 2d array combinations and try all possible keywords
|
||||
# compare output to original functions
|
||||
for i in range(16):
|
||||
for j in range(16):
|
||||
for axes_keyword in [0, 1, None, (0,), (0, 1)]:
|
||||
inp = np.random.rand(i, j)
|
||||
|
||||
assert_array_almost_equal(fft.fftshift(inp, axes_keyword),
|
||||
original_fftshift(inp, axes_keyword))
|
||||
|
||||
assert_array_almost_equal(fft.ifftshift(inp, axes_keyword),
|
||||
original_ifftshift(inp, axes_keyword))
|
||||
|
||||
|
||||
class TestFFTFreq:
|
||||
|
||||
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:
|
||||
|
||||
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:
|
||||
|
||||
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)
|
||||
@@ -0,0 +1,307 @@
|
||||
import numpy as np
|
||||
import pytest
|
||||
from numpy.random import random
|
||||
from numpy.testing import (
|
||||
assert_array_equal, assert_raises, assert_allclose
|
||||
)
|
||||
import threading
|
||||
import 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:
|
||||
|
||||
def test_fft_n(self):
|
||||
assert_raises(ValueError, np.fft.fft, [1, 2, 3], 0)
|
||||
|
||||
|
||||
class TestFFT1D:
|
||||
|
||||
def test_identity(self):
|
||||
maxlen = 512
|
||||
x = random(maxlen) + 1j*random(maxlen)
|
||||
xr = random(maxlen)
|
||||
for i in range(1, maxlen):
|
||||
assert_allclose(np.fft.ifft(np.fft.fft(x[0:i])), x[0:i],
|
||||
atol=1e-12)
|
||||
assert_allclose(np.fft.irfft(np.fft.rfft(xr[0:i]), i),
|
||||
xr[0:i], atol=1e-12)
|
||||
|
||||
def test_fft(self):
|
||||
x = random(30) + 1j*random(30)
|
||||
assert_allclose(fft1(x), np.fft.fft(x), atol=1e-6)
|
||||
assert_allclose(fft1(x), np.fft.fft(x, norm="backward"), atol=1e-6)
|
||||
assert_allclose(fft1(x) / np.sqrt(30),
|
||||
np.fft.fft(x, norm="ortho"), atol=1e-6)
|
||||
assert_allclose(fft1(x) / 30.,
|
||||
np.fft.fft(x, norm="forward"), atol=1e-6)
|
||||
|
||||
@pytest.mark.parametrize('norm', (None, 'backward', 'ortho', 'forward'))
|
||||
def test_ifft(self, norm):
|
||||
x = random(30) + 1j*random(30)
|
||||
assert_allclose(
|
||||
x, np.fft.ifft(np.fft.fft(x, norm=norm), norm=norm),
|
||||
atol=1e-6)
|
||||
# Ensure we get the correct error message
|
||||
with pytest.raises(ValueError,
|
||||
match='Invalid number of FFT data points'):
|
||||
np.fft.ifft([], norm=norm)
|
||||
|
||||
def test_fft2(self):
|
||||
x = random((30, 20)) + 1j*random((30, 20))
|
||||
assert_allclose(np.fft.fft(np.fft.fft(x, axis=1), axis=0),
|
||||
np.fft.fft2(x), atol=1e-6)
|
||||
assert_allclose(np.fft.fft2(x),
|
||||
np.fft.fft2(x, norm="backward"), atol=1e-6)
|
||||
assert_allclose(np.fft.fft2(x) / np.sqrt(30 * 20),
|
||||
np.fft.fft2(x, norm="ortho"), atol=1e-6)
|
||||
assert_allclose(np.fft.fft2(x) / (30. * 20.),
|
||||
np.fft.fft2(x, norm="forward"), atol=1e-6)
|
||||
|
||||
def test_ifft2(self):
|
||||
x = random((30, 20)) + 1j*random((30, 20))
|
||||
assert_allclose(np.fft.ifft(np.fft.ifft(x, axis=1), axis=0),
|
||||
np.fft.ifft2(x), atol=1e-6)
|
||||
assert_allclose(np.fft.ifft2(x),
|
||||
np.fft.ifft2(x, norm="backward"), atol=1e-6)
|
||||
assert_allclose(np.fft.ifft2(x) * np.sqrt(30 * 20),
|
||||
np.fft.ifft2(x, norm="ortho"), atol=1e-6)
|
||||
assert_allclose(np.fft.ifft2(x) * (30. * 20.),
|
||||
np.fft.ifft2(x, norm="forward"), atol=1e-6)
|
||||
|
||||
def test_fftn(self):
|
||||
x = random((30, 20, 10)) + 1j*random((30, 20, 10))
|
||||
assert_allclose(
|
||||
np.fft.fft(np.fft.fft(np.fft.fft(x, axis=2), axis=1), axis=0),
|
||||
np.fft.fftn(x), atol=1e-6)
|
||||
assert_allclose(np.fft.fftn(x),
|
||||
np.fft.fftn(x, norm="backward"), atol=1e-6)
|
||||
assert_allclose(np.fft.fftn(x) / np.sqrt(30 * 20 * 10),
|
||||
np.fft.fftn(x, norm="ortho"), atol=1e-6)
|
||||
assert_allclose(np.fft.fftn(x) / (30. * 20. * 10.),
|
||||
np.fft.fftn(x, norm="forward"), atol=1e-6)
|
||||
|
||||
def test_ifftn(self):
|
||||
x = random((30, 20, 10)) + 1j*random((30, 20, 10))
|
||||
assert_allclose(
|
||||
np.fft.ifft(np.fft.ifft(np.fft.ifft(x, axis=2), axis=1), axis=0),
|
||||
np.fft.ifftn(x), atol=1e-6)
|
||||
assert_allclose(np.fft.ifftn(x),
|
||||
np.fft.ifftn(x, norm="backward"), atol=1e-6)
|
||||
assert_allclose(np.fft.ifftn(x) * np.sqrt(30 * 20 * 10),
|
||||
np.fft.ifftn(x, norm="ortho"), atol=1e-6)
|
||||
assert_allclose(np.fft.ifftn(x) * (30. * 20. * 10.),
|
||||
np.fft.ifftn(x, norm="forward"), atol=1e-6)
|
||||
|
||||
def test_rfft(self):
|
||||
x = random(30)
|
||||
for n in [x.size, 2*x.size]:
|
||||
for norm in [None, 'backward', 'ortho', 'forward']:
|
||||
assert_allclose(
|
||||
np.fft.fft(x, n=n, norm=norm)[:(n//2 + 1)],
|
||||
np.fft.rfft(x, n=n, norm=norm), atol=1e-6)
|
||||
assert_allclose(
|
||||
np.fft.rfft(x, n=n),
|
||||
np.fft.rfft(x, n=n, norm="backward"), atol=1e-6)
|
||||
assert_allclose(
|
||||
np.fft.rfft(x, n=n) / np.sqrt(n),
|
||||
np.fft.rfft(x, n=n, norm="ortho"), atol=1e-6)
|
||||
assert_allclose(
|
||||
np.fft.rfft(x, n=n) / n,
|
||||
np.fft.rfft(x, n=n, norm="forward"), atol=1e-6)
|
||||
|
||||
def test_irfft(self):
|
||||
x = random(30)
|
||||
assert_allclose(x, np.fft.irfft(np.fft.rfft(x)), atol=1e-6)
|
||||
assert_allclose(x, np.fft.irfft(np.fft.rfft(x, norm="backward"),
|
||||
norm="backward"), atol=1e-6)
|
||||
assert_allclose(x, np.fft.irfft(np.fft.rfft(x, norm="ortho"),
|
||||
norm="ortho"), atol=1e-6)
|
||||
assert_allclose(x, np.fft.irfft(np.fft.rfft(x, norm="forward"),
|
||||
norm="forward"), atol=1e-6)
|
||||
|
||||
def test_rfft2(self):
|
||||
x = random((30, 20))
|
||||
assert_allclose(np.fft.fft2(x)[:, :11], np.fft.rfft2(x), atol=1e-6)
|
||||
assert_allclose(np.fft.rfft2(x),
|
||||
np.fft.rfft2(x, norm="backward"), atol=1e-6)
|
||||
assert_allclose(np.fft.rfft2(x) / np.sqrt(30 * 20),
|
||||
np.fft.rfft2(x, norm="ortho"), atol=1e-6)
|
||||
assert_allclose(np.fft.rfft2(x) / (30. * 20.),
|
||||
np.fft.rfft2(x, norm="forward"), atol=1e-6)
|
||||
|
||||
def test_irfft2(self):
|
||||
x = random((30, 20))
|
||||
assert_allclose(x, np.fft.irfft2(np.fft.rfft2(x)), atol=1e-6)
|
||||
assert_allclose(x, np.fft.irfft2(np.fft.rfft2(x, norm="backward"),
|
||||
norm="backward"), atol=1e-6)
|
||||
assert_allclose(x, np.fft.irfft2(np.fft.rfft2(x, norm="ortho"),
|
||||
norm="ortho"), atol=1e-6)
|
||||
assert_allclose(x, np.fft.irfft2(np.fft.rfft2(x, norm="forward"),
|
||||
norm="forward"), atol=1e-6)
|
||||
|
||||
def test_rfftn(self):
|
||||
x = random((30, 20, 10))
|
||||
assert_allclose(np.fft.fftn(x)[:, :, :6], np.fft.rfftn(x), atol=1e-6)
|
||||
assert_allclose(np.fft.rfftn(x),
|
||||
np.fft.rfftn(x, norm="backward"), atol=1e-6)
|
||||
assert_allclose(np.fft.rfftn(x) / np.sqrt(30 * 20 * 10),
|
||||
np.fft.rfftn(x, norm="ortho"), atol=1e-6)
|
||||
assert_allclose(np.fft.rfftn(x) / (30. * 20. * 10.),
|
||||
np.fft.rfftn(x, norm="forward"), atol=1e-6)
|
||||
|
||||
def test_irfftn(self):
|
||||
x = random((30, 20, 10))
|
||||
assert_allclose(x, np.fft.irfftn(np.fft.rfftn(x)), atol=1e-6)
|
||||
assert_allclose(x, np.fft.irfftn(np.fft.rfftn(x, norm="backward"),
|
||||
norm="backward"), atol=1e-6)
|
||||
assert_allclose(x, np.fft.irfftn(np.fft.rfftn(x, norm="ortho"),
|
||||
norm="ortho"), atol=1e-6)
|
||||
assert_allclose(x, np.fft.irfftn(np.fft.rfftn(x, norm="forward"),
|
||||
norm="forward"), atol=1e-6)
|
||||
|
||||
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_allclose(np.fft.fft(x), np.fft.hfft(x_herm), atol=1e-6)
|
||||
assert_allclose(np.fft.hfft(x_herm),
|
||||
np.fft.hfft(x_herm, norm="backward"), atol=1e-6)
|
||||
assert_allclose(np.fft.hfft(x_herm) / np.sqrt(30),
|
||||
np.fft.hfft(x_herm, norm="ortho"), atol=1e-6)
|
||||
assert_allclose(np.fft.hfft(x_herm) / 30.,
|
||||
np.fft.hfft(x_herm, norm="forward"), atol=1e-6)
|
||||
|
||||
def test_ihfft(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_allclose(x_herm, np.fft.ihfft(np.fft.hfft(x_herm)), atol=1e-6)
|
||||
assert_allclose(x_herm, np.fft.ihfft(np.fft.hfft(x_herm,
|
||||
norm="backward"), norm="backward"), atol=1e-6)
|
||||
assert_allclose(x_herm, np.fft.ihfft(np.fft.hfft(x_herm,
|
||||
norm="ortho"), norm="ortho"), atol=1e-6)
|
||||
assert_allclose(x_herm, np.fft.ihfft(np.fft.hfft(x_herm,
|
||||
norm="forward"), norm="forward"), atol=1e-6)
|
||||
|
||||
@pytest.mark.parametrize("op", [np.fft.fftn, np.fft.ifftn,
|
||||
np.fft.rfftn, np.fft.irfftn])
|
||||
def test_axes(self, op):
|
||||
x = random((30, 20, 10))
|
||||
axes = [(0, 1, 2), (0, 2, 1), (1, 0, 2), (1, 2, 0), (2, 0, 1), (2, 1, 0)]
|
||||
for a in axes:
|
||||
op_tr = op(np.transpose(x, a))
|
||||
tr_op = np.transpose(op(x, axes=a), a)
|
||||
assert_allclose(op_tr, tr_op, atol=1e-6)
|
||||
|
||||
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, 'backward', 'ortho', 'forward']:
|
||||
tmp = forw(x, n=n, norm=norm)
|
||||
tmp = back(tmp, n=n, norm=norm)
|
||||
assert_allclose(x_norm,
|
||||
np.linalg.norm(tmp), atol=1e-6)
|
||||
|
||||
@pytest.mark.parametrize("dtype", [np.half, np.single, np.double,
|
||||
np.longdouble])
|
||||
def test_dtypes(self, dtype):
|
||||
# make sure that all input precisions are accepted and internally
|
||||
# converted to 64bit
|
||||
x = random(30).astype(dtype)
|
||||
assert_allclose(np.fft.ifft(np.fft.fft(x)), x, atol=1e-6)
|
||||
assert_allclose(np.fft.irfft(np.fft.rfft(x)), x, atol=1e-6)
|
||||
|
||||
|
||||
@pytest.mark.parametrize(
|
||||
"dtype",
|
||||
[np.float32, np.float64, np.complex64, np.complex128])
|
||||
@pytest.mark.parametrize("order", ["F", 'non-contiguous'])
|
||||
@pytest.mark.parametrize(
|
||||
"fft",
|
||||
[np.fft.fft, np.fft.fft2, np.fft.fftn,
|
||||
np.fft.ifft, np.fft.ifft2, np.fft.ifftn])
|
||||
def test_fft_with_order(dtype, order, fft):
|
||||
# Check that FFT/IFFT produces identical results for C, Fortran and
|
||||
# non contiguous arrays
|
||||
rng = np.random.RandomState(42)
|
||||
X = rng.rand(8, 7, 13).astype(dtype, copy=False)
|
||||
# See discussion in pull/14178
|
||||
_tol = 8.0 * np.sqrt(np.log2(X.size)) * np.finfo(X.dtype).eps
|
||||
if order == 'F':
|
||||
Y = np.asfortranarray(X)
|
||||
else:
|
||||
# Make a non contiguous array
|
||||
Y = X[::-1]
|
||||
X = np.ascontiguousarray(X[::-1])
|
||||
|
||||
if fft.__name__.endswith('fft'):
|
||||
for axis in range(3):
|
||||
X_res = fft(X, axis=axis)
|
||||
Y_res = fft(Y, axis=axis)
|
||||
assert_allclose(X_res, Y_res, atol=_tol, rtol=_tol)
|
||||
elif fft.__name__.endswith(('fft2', 'fftn')):
|
||||
axes = [(0, 1), (1, 2), (0, 2)]
|
||||
if fft.__name__.endswith('fftn'):
|
||||
axes.extend([(0,), (1,), (2,), None])
|
||||
for ax in axes:
|
||||
X_res = fft(X, axes=ax)
|
||||
Y_res = fft(Y, axes=ax)
|
||||
assert_allclose(X_res, Y_res, atol=_tol, rtol=_tol)
|
||||
else:
|
||||
raise ValueError()
|
||||
|
||||
|
||||
class TestFFTThreadSafe:
|
||||
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)
|
||||
Reference in New Issue
Block a user