diff --git a/CHANGELOG.md b/CHANGELOG.md index 5d752a8..ebb5d30 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,13 +1,17 @@ # changelog All notable changes to this project will be documented in this file. -The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), +The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.1.0/), and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). ## [dev] (MM/DD/YY) ### Added -* SciPy interface `mkl_fft.interfaces.scipy_fft` now includes Hermitian FFT functions: `hfft`, `ihfft`, `hfftn`, `ihfftn`, `hfft2`, and `ihfft2` [gh-161](https://github.com/IntelPython/mkl_fft/pull/161) +* Added Hermitian FFT functions to SciPy interface `mkl_fft.interfaces.scipy_fft`: `hfft`, `ihfft`, `hfftn`, `ihfftn`, `hfft2`, and `ihfft2` [gh-161](https://github.com/IntelPython/mkl_fft/pull/161) +* Added support for `out` kwarg to all FFT functions in `mkl_fft` and `mkl_fft.interfaces.numpy_fft` [gh-157](https://github.com/IntelPython/mkl_fft/pull/157) + +### Changed +* NumPy interface `mkl_fft.interfaces.numpy_fft` is aligned with numpy-2.* [gh-139](https://github.com/IntelPython/mkl_fft/pull/139), [gh-157](https://github.com/IntelPython/mkl_fft/pull/157) ## [1.3.14] (04/10/2025) diff --git a/README.md b/README.md index c38285b..c7356ab 100644 --- a/README.md +++ b/README.md @@ -1,11 +1,11 @@ -## `mkl_fft` -- a NumPy-based Python interface to Intel (R) MKL FFT functionality +## `mkl_fft` -- a NumPy-based Python interface to Intel® oneAPI Math Kernel Library (OneMKL) FFT functionality [![Conda package](https://github.com/IntelPython/mkl_fft/actions/workflows/conda-package.yml/badge.svg)](https://github.com/IntelPython/mkl_fft/actions/workflows/conda-package.yml) [![Editable build using pip and pre-release NumPy](https://github.com/IntelPython/mkl_fft/actions/workflows/build_pip.yaml/badge.svg)](https://github.com/IntelPython/mkl_fft/actions/workflows/build_pip.yaml) [![Conda package with conda-forge channel only](https://github.com/IntelPython/mkl_fft/actions/workflows/conda-package-cf.yml/badge.svg)](https://github.com/IntelPython/mkl_fft/actions/workflows/conda-package-cf.yml) [![OpenSSF Scorecard](https://api.securityscorecards.dev/projects/github.com/IntelPython/mkl_fft/badge)](https://securityscorecards.dev/viewer/?uri=github.com/IntelPython/mkl_fft) -`mkl_fft` started as a part of Intel (R) Distribution for Python* optimizations to NumPy, and is now being released -as a stand-alone package. It can be installed into conda environment using +`mkl_fft` started as a part of Intel® Distribution for Python* optimizations to NumPy, and is now being released +as a stand-alone package. It can be installed into conda environment from Intel's channel using: ``` conda install -c https://software.repos.intel.com/python/conda mkl_fft @@ -19,13 +19,13 @@ or from conda-forge channel: --- -To install mkl_fft Pypi package please use following command: +To install `mkl_fft` PyPI package please use following command: ``` python -m pip install --index-url https://software.repos.intel.com/python/pypi --extra-index-url https://pypi.org/simple mkl_fft ``` -If command above installs NumPy package from the Pypi, please use following command to install Intel optimized NumPy wheel package from Intel Pypi Cloud: +If command above installs NumPy package from the PyPI, please use following command to install Intel optimized NumPy wheel package from Intel PyPI Cloud: ``` python -m pip install --index-url https://software.repos.intel.com/python/pypi --extra-index-url https://pypi.org/simple mkl_fft numpy== @@ -35,7 +35,7 @@ Where `` should be the latest version from https://software.repos --- -Since MKL FFT supports performing discrete Fourier transforms over non-contiguously laid out arrays, MKL can be directly +Since MKL FFT supports performing discrete Fourier transforms over non-contiguously laid out arrays, OneMKL can be directly used on any well-behaved floating point array with no internal overlaps for both in-place and not in-place transforms of arrays in single and double floating point precision. @@ -50,38 +50,47 @@ More details can be found in SciPy 2017 conference proceedings: `mkl_fft` implements the following functions: -### Complex transforms, similar to those in `scipy.fft`: +### complex-to-complex (c2c) transforms: -`fft(x, n=None, axis=-1, overwrite_x=False)` +`fft(x, n=None, axis=-1, overwrite_x=False, fwd_scale=1.0, out=out)` - 1D FFT, similar to `scipy.fft.fft` -`ifft(x, n=None, axis=-1, overwrite_x=False)` +`fft2(x, s=None, axes=(-2, -1), overwrite_x=False, fwd_scale=1.0, out=out)` - 2D FFT, similar to `scipy.fft.fft2` -`fft2(x, shape=None, axes=(-2,-1), overwrite_x=False)` +`fftn(x, s=None, axes=None, overwrite_x=False, fwd_scale=1.0, out=out)` - ND FFT, similar to `scipy.fft.fftn` -`ifft2(x, shape=None, axes=(-2,-1), overwrite_x=False)` +and similar inverse FFT (`ifft*`) functions. -`fftn(x, n=None, axes=None, overwrite_x=False)` +### real-to-complex (r2c) and complex-to-real (c2r) transforms: -`ifftn(x, n=None, axes=None, overwrite_x=False)` +`rfft(x, n=None, axis=-1, fwd_scale=1.0, out=out)` - r2c 1D FFT, similar to `numpy.fft.rfft` -### Real transforms +`rfft2(x, s=None, axes=(-2, -1), fwd_scale=1.0, out=out)` - r2c 2D FFT, similar to `numpy.fft.rfft2` -`rfftpack(x, n=None, axis=-1, overwrite_x=False)` - real 1D Fourier transform, like `scipy.fftpack.rfft` +`rfftn(x, s=None, axes=None, fwd_scale=1.0, out=out)` - r2c ND FFT, similar to `numpy.fft.rfftn` -`rfft(x, n=None, axis=-1)` - real 1D Fourier transform, like `numpy.fft.rfft` +and similar inverse c2r FFT (`irfft*`) functions. -`rfft2(x, s=None, axes=(-2,-1))` - real 2D Fourier transform, like `numpy.fft.rfft2` - -`rfftn(x, s=None, axes=None)` - real ND Fourier transform, like `numpy.fft.rfftn` - -... and similar `irfft*` functions. - - -The package also provides `mkl_fft.interfaces.numpy_fft` and `mkl_fft.interfaces.scipy_fft` interfaces which provide drop-in replacements for equivalent functions in NumPy and SciPy respectively. +The package also provides `mkl_fft.interfaces.numpy_fft` and `mkl_fft.interfaces.scipy_fft` interfaces which provide drop-in replacements for equivalent functions in NumPy and SciPy, respectively. --- -To build `mkl_fft` from sources on Linux: - - install a recent version of MKL, if necessary; - - execute `source /path_to_oneapi/mkl/latest/env/vars.sh`; - - execute `python -m pip install .` +To build `mkl_fft` from sources on Linux with Intel® OneMKL: + - create a virtual environment: `python3 -m venv fft_env` + - activate the environment: `source fft_env/bin/activate` + - install a recent version of OneMKL, if necessary + - execute `source /path_to_oneapi/mkl/latest/env/vars.sh` + - `git clone https://github.com/IntelPython/mkl_fft.git mkl_fft` + - `cd mkl_fft` + - `python -m pip install .` + - `cd ..` + - `python -c "import mkl_fft"` + +To build `mkl_fft` from sources on Linux with conda follow these steps: + - `conda create -n fft_env python=3.12 mkl-devel` + - `conda activate fft_env` + - `export MKLROOT=$CONDA_PREFIX` + - `git clone https://github.com/IntelPython/mkl_fft.git mkl_fft` + - `cd mkl_fft` + - `python -m pip install .` + - `cd ..` + - `python -c "import mkl_fft"` diff --git a/conda-recipe-cf/meta.yaml b/conda-recipe-cf/meta.yaml index af71883..73da10d 100644 --- a/conda-recipe-cf/meta.yaml +++ b/conda-recipe-cf/meta.yaml @@ -43,4 +43,4 @@ about: home: http://github.com/IntelPython/mkl_fft license: BSD-3-Clause license_file: LICENSE.txt - summary: NumPy-based implementation of Fast Fourier Transform using Intel (R) Math Kernel Library + summary: NumPy-based implementation of Fast Fourier Transform using Intel® oneAPI Math Kernel Library (OneMKL) diff --git a/conda-recipe/meta.yaml b/conda-recipe/meta.yaml index 3d53717..75f227a 100644 --- a/conda-recipe/meta.yaml +++ b/conda-recipe/meta.yaml @@ -43,4 +43,4 @@ about: home: http://github.com/IntelPython/mkl_fft license: BSD-3-Clause license_file: LICENSE.txt - summary: NumPy-based implementation of Fast Fourier Transform using Intel (R) Math Kernel Library + summary: NumPy-based implementation of Fast Fourier Transform using Intel® oneAPI Math Kernel Library (OneMKL) diff --git a/mkl_fft/_numpy_fft.py b/mkl_fft/_numpy_fft.py index fe7a4c6..9ea1203 100644 --- a/mkl_fft/_numpy_fft.py +++ b/mkl_fft/_numpy_fft.py @@ -23,34 +23,10 @@ # CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, # OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE # OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. -""" -Discrete Fourier Transforms - -Routines in this module: - -fft(a, n=None, axis=-1, norm=None) -ifft(a, n=None, axis=-1, norm=None) -rfft(a, n=None, axis=-1, norm=None) -irfft(a, n=None, axis=-1, norm=None) -hfft(a, n=None, axis=-1, norm=None) -ihfft(a, n=None, axis=-1, norm=None) -fftn(a, s=None, axes=None, norm=None) -ifftn(a, s=None, axes=None, norm=None) -rfftn(a, s=None, axes=None, norm=None) -irfftn(a, s=None, axes=None, norm=None) -fft2(a, s=None, axes=(-2,-1), norm=None) -ifft2(a, s=None, axes=(-2, -1), norm=None) -rfft2(a, s=None, axes=(-2,-1), norm=None) -irfft2(a, s=None, axes=(-2, -1), norm=None) - -i = inverse transform -r = transform of purely real data -h = Hermite transform -n = n-dimensional transform -2 = 2-dimensional transform -(Note: 2D routines are just nD routines with different default -behavior.) +""" +An interface for FFT module of NumPy (`numpy.fft`) that uses OneMKL FFT +in the backend. """ __all__ = [ @@ -93,7 +69,7 @@ def _cook_nd_args(a, s=None, axes=None, invreal=False): shapeless = False s = list(s) if axes is None: - if not shapeless and np.__version__ >= "2.0": + if not shapeless and np.lib.NumpyVersion(np.__version__) >= "2.0.0": msg = ( "`axes` should not be `None` if `s` is not `None` " "(Deprecated in NumPy 2.0). In a future version of NumPy, " @@ -108,7 +84,7 @@ def _cook_nd_args(a, s=None, axes=None, invreal=False): raise ValueError("Shape and axes have different lengths.") if invreal and shapeless: s[-1] = (a.shape[axes[-1]] - 1) * 2 - if None in s and np.__version__ >= "2.0": + if None in s and np.lib.NumpyVersion(np.__version__) >= "2.0.0": msg = ( "Passing an array containing `None` values to `s` is " "deprecated in NumPy 2.0 and will raise an error in " @@ -135,461 +111,76 @@ def trycall(func, args, kwrds): return res -def fft(a, n=None, axis=-1, norm=None): +def fft(a, n=None, axis=-1, norm=None, out=None): """ Compute the one-dimensional discrete Fourier Transform. - This function computes the one-dimensional *n*-point discrete Fourier - Transform (DFT) with the efficient Fast Fourier Transform (FFT) - algorithm [CT]. - - Parameters - ---------- - a : array_like - Input array, can be complex. - n : int, optional - Length of the transformed axis of the output. - If `n` is smaller than the length of the input, the input is cropped. - If it is larger, the input is padded with zeros. If `n` is not given, - the length of the input along the axis specified by `axis` is used. - axis : int, optional - Axis over which to compute the FFT. If not given, the last axis is - used. - norm : {"backward", "ortho", "forward"}, optional - .. versionadded:: 1.10.0 - - Normalization mode (see `numpy.fft`). Default is "backward". - Indicates which direction of the forward/backward pair of transforms - is scaled and with what normalization factor. - - .. versionadded:: 1.20.0 - - The "backward", "forward" values were added. - - Returns - ------- - out : complex ndarray - The truncated or zero-padded input, transformed along the axis - indicated by `axis`, or the last one if `axis` is not specified. - - Raises - ------ - IndexError - if `axes` is larger than the last axis of `a`. - - See Also - -------- - numpy.fft : for definition of the DFT and conventions used. - ifft : The inverse of `fft`. - fft2 : The two-dimensional FFT. - fftn : The *n*-dimensional FFT. - rfftn : The *n*-dimensional FFT of real input. - fftfreq : Frequency bins for given FFT parameters. - - Notes - ----- - FFT (Fast Fourier Transform) refers to a way the discrete Fourier - Transform (DFT) can be calculated efficiently, by using symmetries in the - calculated terms. The symmetry is highest when `n` is a power of 2, and - the transform is therefore most efficient for these sizes. - - The DFT is defined, with the conventions used in this implementation, in - the documentation for the `numpy.fft` module. - - 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. - - Examples - -------- - >>> np.fft.fft(np.exp(2j * np.pi * np.arange(8) / 8)) - array([ -3.44505240e-16 +1.14383329e-17j, - 8.00000000e+00 -5.71092652e-15j, - 2.33482938e-16 +1.22460635e-16j, - 1.64863782e-15 +1.77635684e-15j, - 9.95839695e-17 +2.33482938e-16j, - 0.00000000e+00 +1.66837030e-15j, - 1.14383329e-17 +1.22460635e-16j, - -1.64863782e-15 +1.77635684e-15j]) - - >>> import matplotlib.pyplot as plt - >>> t = np.arange(256) - >>> sp = np.fft.fft(np.sin(t)) - >>> freq = np.fft.fftfreq(t.shape[-1]) - >>> plt.plot(freq, sp.real, freq, sp.imag) - [, ] - >>> plt.show() - - In this example, real input has an FFT which is Hermitian, i.e., symmetric - in the real part and anti-symmetric in the imaginary part, as described in - the `numpy.fft` documentation. + For full documentation refer to `numpy.fft.fft`. """ - x = _downcast_float128_array(a) fsc = _compute_fwd_scale(norm, n, x.shape[axis]) - return trycall(mkl_fft.fft, (x,), {"n": n, "axis": axis, "fwd_scale": fsc}) + return trycall( + mkl_fft.fft, (x,), {"n": n, "axis": axis, "fwd_scale": fsc, "out": out} + ) -def ifft(a, n=None, axis=-1, norm=None): +def ifft(a, n=None, axis=-1, norm=None, out=None): """ Compute the one-dimensional inverse discrete Fourier Transform. - This function computes the inverse of the one-dimensional *n*-point - discrete Fourier transform computed by `fft`. In other words, - ``ifft(fft(a)) == a`` to within numerical accuracy. - For a general description of the algorithm and definitions, - see `numpy.fft`. - - The input should be ordered in the same way as is returned by `fft`, - i.e., - - * ``a[0]`` should contain the zero frequency term, - * ``a[1:n//2]`` should contain the positive-frequency terms, - * ``a[n//2 + 1:]`` should contain the negative-frequency terms, in - increasing order starting from the most negative frequency. - - Parameters - ---------- - a : array_like - Input array, can be complex. - n : int, optional - Length of the transformed axis of the output. - If `n` is smaller than the length of the input, the input is cropped. - If it is larger, the input is padded with zeros. If `n` is not given, - the length of the input along the axis specified by `axis` is used. - See notes about padding issues. - axis : int, optional - Axis over which to compute the inverse DFT. If not given, the last - axis is used. - norm : {"backward", "ortho", "forward"}, optional - .. versionadded:: 1.10.0 - - Normalization mode (see `numpy.fft`). Default is "backward". - Indicates which direction of the forward/backward pair of transforms - is scaled and with what normalization factor. - - .. versionadded:: 1.20.0 - - The "backward", "forward" values were added. - - Returns - ------- - out : complex ndarray - The truncated or zero-padded input, transformed along the axis - indicated by `axis`, or the last one if `axis` is not specified. - - Raises - ------ - IndexError - If `axes` is larger than the last axis of `a`. - - See Also - -------- - numpy.fft : An introduction, with definitions and general explanations. - fft : The one-dimensional (forward) FFT, of which `ifft` is the inverse - ifft2 : The two-dimensional inverse FFT. - ifftn : The n-dimensional inverse FFT. - - Notes - ----- - If the input parameter `n` is larger than the size of the input, the input - is padded by appending zeros at the end. Even though this is the common - approach, it might lead to surprising results. If a different padding is - desired, it must be performed before calling `ifft`. - - Examples - -------- - >>> np.fft.ifft([0, 4, 0, 0]) - array([ 1.+0.j, 0.+1.j, -1.+0.j, 0.-1.j]) - - Create and plot a band-limited signal with random phases: - - >>> import matplotlib.pyplot as plt - >>> t = np.arange(400) - >>> n = np.zeros((400,), dtype=complex) - >>> n[40:60] = np.exp(1j*np.random.uniform(0, 2*np.pi, (20,))) - >>> s = np.fft.ifft(n) - >>> plt.plot(t, s.real, 'b-', t, s.imag, 'r--') - ... - >>> plt.legend(('real', 'imaginary')) - ... - >>> plt.show() + For full documentation refer to `numpy.fft.ifft`. """ - x = _downcast_float128_array(a) fsc = _compute_fwd_scale(norm, n, x.shape[axis]) - return trycall(mkl_fft.ifft, (x,), {"n": n, "axis": axis, "fwd_scale": fsc}) + return trycall( + mkl_fft.ifft, (x,), {"n": n, "axis": axis, "fwd_scale": fsc, "out": out} + ) -def rfft(a, n=None, axis=-1, norm=None): +def rfft(a, n=None, axis=-1, norm=None, out=None): """ Compute the one-dimensional discrete Fourier Transform for real input. - This function computes the one-dimensional *n*-point discrete Fourier - Transform (DFT) of a real-valued array by means of an efficient algorithm - called the Fast Fourier Transform (FFT). - - Parameters - ---------- - a : array_like - Input array - n : int, optional - Number of points along transformation axis in the input to use. - If `n` is smaller than the length of the input, the input is cropped. - If it is larger, the input is padded with zeros. If `n` is not given, - the length of the input along the axis specified by `axis` is used. - axis : int, optional - Axis over which to compute the FFT. If not given, the last axis is - used. - norm : {"backward", "ortho", "forward"}, optional - .. versionadded:: 1.10.0 - - Normalization mode (see `numpy.fft`). Default is "backward". - Indicates which direction of the forward/backward pair of transforms - is scaled and with what normalization factor. - - .. versionadded:: 1.20.0 - - The "backward", "forward" values were added. - - Returns - ------- - out : complex ndarray - The truncated or zero-padded input, transformed along the axis - indicated by `axis`, or the last one if `axis` is not specified. - If `n` is even, the length of the transformed axis is ``(n/2)+1``. - If `n` is odd, the length is ``(n+1)/2``. - - Raises - ------ - IndexError - If `axis` is larger than the last axis of `a`. - - See Also - -------- - numpy.fft : For definition of the DFT and conventions used. - irfft : The inverse of `rfft`. - fft : The one-dimensional FFT of general (complex) input. - fftn : The *n*-dimensional FFT. - rfftn : The *n*-dimensional FFT of real input. - - Notes - ----- - When the DFT is computed for purely real input, the output is - Hermitian-symmetric, i.e. the negative frequency terms are just the complex - conjugates of the corresponding positive-frequency terms, and the - negative-frequency terms are therefore redundant. This function does not - compute the negative frequency terms, and the length of the transformed - axis of the output is therefore ``n//2 + 1``. - - When ``A = rfft(a)`` and fs is the sampling frequency, ``A[0]`` contains - the zero-frequency term 0*fs, which is real due to Hermitian symmetry. - - If `n` is even, ``A[-1]`` contains the term representing both positive - and negative Nyquist frequency (+fs/2 and -fs/2), and must also be purely - real. If `n` is odd, there is no term at fs/2; ``A[-1]`` contains - the largest positive frequency (fs/2*(n-1)/n), and is complex in the - general case. - - If the input `a` contains an imaginary part, it is silently discarded. - - Examples - -------- - >>> np.fft.fft([0, 1, 0, 0]) - array([ 1.+0.j, 0.-1.j, -1.+0.j, 0.+1.j]) - >>> np.fft.rfft([0, 1, 0, 0]) - array([ 1.+0.j, 0.-1.j, -1.+0.j]) - - Notice how the final element of the `fft` output is the complex conjugate - of the second element, for real input. For `rfft`, this symmetry is - exploited to compute only the non-negative frequency terms. + For full documentation refer to `numpy.fft.rfft`. """ - x = _downcast_float128_array(a) fsc = _compute_fwd_scale(norm, n, x.shape[axis]) - return trycall(mkl_fft.rfft, (x,), {"n": n, "axis": axis, "fwd_scale": fsc}) + return trycall( + mkl_fft.rfft, (x,), {"n": n, "axis": axis, "fwd_scale": fsc, "out": out} + ) -def irfft(a, n=None, axis=-1, norm=None): +def irfft(a, n=None, axis=-1, norm=None, out=None): """ - Compute the inverse of the n-point DFT for real input. - - This function computes the inverse of the one-dimensional *n*-point - discrete Fourier Transform of real input computed by `rfft`. - In other words, ``irfft(rfft(a), len(a)) == a`` to within numerical - accuracy. (See Notes below for why ``len(a)`` is necessary here.) - - The input is expected to be in the form returned by `rfft`, i.e. the - real zero-frequency term followed by the complex positive frequency terms - in order of increasing frequency. Since the discrete Fourier Transform of - real input is Hermitian-symmetric, the negative frequency terms are taken - to be the complex conjugates of the corresponding positive frequency terms. - - Parameters - ---------- - a : array_like - The input array. - n : int, optional - Length of the transformed axis of the output. - For `n` output points, ``n//2+1`` input points are necessary. If the - input is longer than this, it is cropped. If it is shorter than this, - it is padded with zeros. If `n` is not given, it is determined from - the length of the input along the axis specified by `axis`. - axis : int, optional - Axis over which to compute the inverse FFT. If not given, the last - axis is used. - norm : {"backward", "ortho", "forward"}, optional - .. versionadded:: 1.10.0 - - Normalization mode (see `numpy.fft`). Default is "backward". - Indicates which direction of the forward/backward pair of transforms - is scaled and with what normalization factor. - - .. versionadded:: 1.20.0 - - The "backward", "forward" values were added. - - Returns - ------- - out : ndarray - The truncated or zero-padded input, transformed along the axis - indicated by `axis`, or the last one if `axis` is not specified. - The length of the transformed axis is `n`, or, if `n` is not given, - ``2*(m-1)`` where ``m`` is the length of the transformed axis of the - input. To get an odd number of output points, `n` must be specified. - - Raises - ------ - IndexError - If `axis` is larger than the last axis of `a`. - - See Also - -------- - numpy.fft : For definition of the DFT and conventions used. - rfft : The one-dimensional FFT of real input, of which `irfft` is inverse. - fft : The one-dimensional FFT. - irfft2 : The inverse of the two-dimensional FFT of real input. - irfftn : The inverse of the *n*-dimensional FFT of real input. - - Notes - ----- - Returns the real valued `n`-point inverse discrete Fourier transform - of `a`, where `a` contains the non-negative frequency terms of a - Hermitian-symmetric sequence. `n` is the length of the result, not the - input. - - If you specify an `n` such that `a` must be zero-padded or truncated, the - extra/removed values will be added/removed at high frequencies. One can - thus resample a series to `m` points via Fourier interpolation by: - ``a_resamp = irfft(rfft(a), m)``. - - Examples - -------- - >>> np.fft.ifft([1, -1j, -1, 1j]) - array([ 0.+0.j, 1.+0.j, 0.+0.j, 0.+0.j]) - >>> np.fft.irfft([1, -1j, -1]) - array([ 0., 1., 0., 0.]) - - Notice how the last term in the input to the ordinary `ifft` is the - complex conjugate of the second term, and the output has zero imaginary - part everywhere. When calling `irfft`, the negative frequencies are not - specified, and the output array is purely real. + Compute the inverse of `rfft`. - """ + For full documentation refer to `numpy.fft.irfft`. + """ x = _downcast_float128_array(a) fsc = _compute_fwd_scale(norm, n, 2 * (x.shape[axis] - 1)) return trycall( - mkl_fft.irfft, (x,), {"n": n, "axis": axis, "fwd_scale": fsc} + mkl_fft.irfft, + (x,), + {"n": n, "axis": axis, "fwd_scale": fsc, "out": out}, ) -def hfft(a, n=None, axis=-1, norm=None): +def hfft(a, n=None, axis=-1, norm=None, out=None): """ - Compute the FFT of a signal which has Hermitian symmetry (real spectrum). - - Parameters - ---------- - a : array_like - The input array. - n : int, optional - Length of the transformed axis of the output. - For `n` output points, ``n//2+1`` input points are necessary. If the - input is longer than this, it is cropped. If it is shorter than this, - it is padded with zeros. If `n` is not given, it is determined from - the length of the input along the axis specified by `axis`. - axis : int, optional - Axis over which to compute the FFT. If not given, the last - axis is used. - norm : {"backward", "ortho", "forward"}, optional - .. versionadded:: 1.10.0 - - Normalization mode (see `numpy.fft`). Default is "backward". - Indicates which direction of the forward/backward pair of transforms - is scaled and with what normalization factor. - - .. versionadded:: 1.20.0 - - The "backward", "forward" values were added. - - Returns - ------- - out : ndarray - The truncated or zero-padded input, transformed along the axis - indicated by `axis`, or the last one if `axis` is not specified. - The length of the transformed axis is `n`, or, if `n` is not given, - ``2*(m-1)`` where ``m`` is the length of the transformed axis of the - input. To get an odd number of output points, `n` must be specified. - - Raises - ------ - IndexError - If `axis` is larger than the last axis of `a`. - - See also - -------- - rfft : Compute the one-dimensional FFT for real input. - ihfft : The inverse of `hfft`. - - Notes - ----- - `hfft`/`ihfft` are a pair analogous to `rfft`/`irfft`, but for the - opposite case: here the signal has Hermitian symmetry in the time domain - and is real in the frequency domain. So here it's `hfft` for which - you must supply the length of the result if it is to be odd: - ``ihfft(hfft(a), len(a)) == a``, within numerical accuracy. - - Examples - -------- - >>> signal = np.array([1, 2, 3, 4, 3, 2]) - >>> np.fft.fft(signal) - array([ 15.+0.j, -4.+0.j, 0.+0.j, -1.-0.j, 0.+0.j, -4.+0.j]) - >>> np.fft.hfft(signal[:4]) # Input first half of signal - array([ 15., -4., 0., -1., 0., -4.]) - >>> np.fft.hfft(signal, 6) # Input entire signal and truncate - array([ 15., -4., 0., -1., 0., -4.]) - - - >>> signal = np.array([[1, 1.j], [-1.j, 2]]) - >>> np.conj(signal.T) - signal # check Hermitian symmetry - array([[ 0.-0.j, 0.+0.j], - [ 0.+0.j, 0.-0.j]]) - >>> freq_spectrum = np.fft.hfft(signal) - >>> freq_spectrum - array([[ 1., 1.], - [ 2., -2.]]) + Compute the FFT of a signal which has Hermitian symmetry, + i.e., a real spectrum.. - """ + For full documentation refer to `numpy.fft.hfft`. + """ norm = _swap_direction(norm) x = _downcast_float128_array(a) x = np.array(x, copy=True) @@ -597,716 +188,118 @@ def hfft(a, n=None, axis=-1, norm=None): fsc = _compute_fwd_scale(norm, n, 2 * (x.shape[axis] - 1)) return trycall( - mkl_fft.irfft, (x,), {"n": n, "axis": axis, "fwd_scale": fsc} + mkl_fft.irfft, + (x,), + {"n": n, "axis": axis, "fwd_scale": fsc, "out": out}, ) -def ihfft(a, n=None, axis=-1, norm=None): +def ihfft(a, n=None, axis=-1, norm=None, out=None): """ Compute the inverse FFT of a signal which has Hermitian symmetry. - Parameters - ---------- - a : array_like - Input array. - n : int, optional - Length of the inverse FFT. - Number of points along transformation axis in the input to use. - If `n` is smaller than the length of the input, the input is cropped. - If it is larger, the input is padded with zeros. If `n` is not given, - the length of the input along the axis specified by `axis` is used. - axis : int, optional - Axis over which to compute the inverse FFT. If not given, the last - axis is used. - norm : {"backward", "ortho", "forward"}, optional - .. versionadded:: 1.10.0 - - Normalization mode (see `numpy.fft`). Default is "backward". - Indicates which direction of the forward/backward pair of transforms - is scaled and with what normalization factor. - - .. versionadded:: 1.20.0 - - The "backward", "forward" values were added. - - Returns - ------- - out : complex ndarray - The truncated or zero-padded input, transformed along the axis - indicated by `axis`, or the last one if `axis` is not specified. - If `n` is even, the length of the transformed axis is ``(n/2)+1``. - If `n` is odd, the length is ``(n+1)/2``. - - See also - -------- - hfft, irfft - - Notes - ----- - `hfft`/`ihfft` are a pair analogous to `rfft`/`irfft`, but for the - opposite case: here the signal has Hermitian symmetry in the time domain - and is real in the frequency domain. So here it's `hfft` for which - you must supply the length of the result if it is to be odd: - ``ihfft(hfft(a), len(a)) == a``, within numerical accuracy. - - Examples - -------- - >>> spectrum = np.array([ 15, -4, 0, -1, 0, -4]) - >>> np.fft.ifft(spectrum) - array([ 1.+0.j, 2.-0.j, 3.+0.j, 4.+0.j, 3.+0.j, 2.-0.j]) - >>> np.fft.ihfft(spectrum) - array([ 1.-0.j, 2.-0.j, 3.-0.j, 4.-0.j]) + For full documentation refer to `numpy.fft.ihfft`. """ - norm = _swap_direction(norm) x = _downcast_float128_array(a) fsc = _compute_fwd_scale(norm, n, x.shape[axis]) output = trycall( - mkl_fft.rfft, (x,), {"n": n, "axis": axis, "fwd_scale": fsc} + mkl_fft.rfft, (x,), {"n": n, "axis": axis, "fwd_scale": fsc, "out": out} ) np.conjugate(output, out=output) return output -def fftn(a, s=None, axes=None, norm=None): +def fftn(a, s=None, axes=None, norm=None, out=None): """ Compute the N-dimensional discrete Fourier Transform. - This function computes the *N*-dimensional discrete Fourier Transform over - any number of axes in an *M*-dimensional array by means of the Fast Fourier - Transform (FFT). - - Parameters - ---------- - a : array_like - Input array, can be complex. - s : sequence of ints, optional - Shape (length of each transformed axis) of the output - (`s[0]` refers to axis 0, `s[1]` to axis 1, etc.). - This corresponds to `n` for `fft(x, n)`. - Along any axis, if the given shape is smaller than that of the input, - the input is cropped. If it is larger, the input is padded with zeros. - if `s` is not given, the shape of the input along the axes specified - by `axes` is used. - axes : sequence of ints, optional - Axes over which to compute the FFT. If not given, the last ``len(s)`` - axes are used, or all axes if `s` is also not specified. - Repeated indices in `axes` means that the transform over that axis is - performed multiple times. - norm : {"backward", "ortho", "forward"}, optional - .. versionadded:: 1.10.0 - - Normalization mode (see `numpy.fft`). Default is "backward". - Indicates which direction of the forward/backward pair of transforms - is scaled and with what normalization factor. - - .. versionadded:: 1.20.0 - - The "backward", "forward" values were added. - - Returns - ------- - out : complex ndarray - The truncated or zero-padded input, transformed along the axes - indicated by `axes`, or by a combination of `s` and `a`, - as explained in the parameters section above. - - Raises - ------ - ValueError - If `s` and `axes` have different length. - IndexError - If an element of `axes` is larger than than the number of axes of `a`. - - See Also - -------- - numpy.fft : Overall view of discrete Fourier transforms, with definitions - and conventions used. - ifftn : The inverse of `fftn`, the inverse *n*-dimensional FFT. - fft : The one-dimensional FFT, with definitions and conventions used. - rfftn : The *n*-dimensional FFT of real input. - fft2 : The two-dimensional FFT. - fftshift : Shifts zero-frequency terms to centre of array - - Notes - ----- - The output, analogously to `fft`, contains the term for zero frequency in - the low-order corner of all axes, the positive frequency terms in the - first half of all axes, the term for the Nyquist frequency in the middle - of all axes and the negative frequency terms in the second half of all - axes, in order of decreasingly negative frequency. - - See `numpy.fft` for details, definitions and conventions used. - - Examples - -------- - >>> a = np.mgrid[:3, :3, :3][0] - >>> np.fft.fftn(a, axes=(1, 2)) - array([[[ 0.+0.j, 0.+0.j, 0.+0.j], - [ 0.+0.j, 0.+0.j, 0.+0.j], - [ 0.+0.j, 0.+0.j, 0.+0.j]], - [[ 9.+0.j, 0.+0.j, 0.+0.j], - [ 0.+0.j, 0.+0.j, 0.+0.j], - [ 0.+0.j, 0.+0.j, 0.+0.j]], - [[ 18.+0.j, 0.+0.j, 0.+0.j], - [ 0.+0.j, 0.+0.j, 0.+0.j], - [ 0.+0.j, 0.+0.j, 0.+0.j]]]) - >>> np.fft.fftn(a, (2, 2), axes=(0, 1)) - array([[[ 2.+0.j, 2.+0.j, 2.+0.j], - [ 0.+0.j, 0.+0.j, 0.+0.j]], - [[-2.+0.j, -2.+0.j, -2.+0.j], - [ 0.+0.j, 0.+0.j, 0.+0.j]]]) - - >>> import matplotlib.pyplot as plt - >>> [X, Y] = np.meshgrid(2 * np.pi * np.arange(200) / 12, - ... 2 * np.pi * np.arange(200) / 34) - >>> S = np.sin(X) + np.cos(Y) + np.random.uniform(0, 1, X.shape) - >>> FS = np.fft.fftn(S) - >>> plt.imshow(np.log(np.abs(np.fft.fftshift(FS))**2)) - - >>> plt.show() + For full documentation refer to `numpy.fft.fftn`. """ - x = _downcast_float128_array(a) s, axes = _cook_nd_args(x, s, axes) fsc = _compute_fwd_scale(norm, s, x.shape) - return trycall(mkl_fft.fftn, (x,), {"s": s, "axes": axes, "fwd_scale": fsc}) + return trycall( + mkl_fft.fftn, (x,), {"s": s, "axes": axes, "fwd_scale": fsc, "out": out} + ) -def ifftn(a, s=None, axes=None, norm=None): +def ifftn(a, s=None, axes=None, norm=None, out=None): """ Compute the N-dimensional inverse discrete Fourier Transform. - This function computes the inverse of the N-dimensional discrete - Fourier Transform over any number of axes in an M-dimensional array by - means of the Fast Fourier Transform (FFT). In other words, - ``ifftn(fftn(a)) == a`` to within numerical accuracy. - For a description of the definitions and conventions used, see `numpy.fft`. - - The input, analogously to `ifft`, should be ordered in the same way as is - returned by `fftn`, i.e. it should have the term for zero frequency - in all axes in the low-order corner, the positive frequency terms in the - first half of all axes, the term for the Nyquist frequency in the middle - of all axes and the negative frequency terms in the second half of all - axes, in order of decreasingly negative frequency. - - Parameters - ---------- - a : array_like - Input array, can be complex. - s : sequence of ints, optional - Shape (length of each transformed axis) of the output - (``s[0]`` refers to axis 0, ``s[1]`` to axis 1, etc.). - This corresponds to ``n`` for ``ifft(x, n)``. - Along any axis, if the given shape is smaller than that of the input, - the input is cropped. If it is larger, the input is padded with zeros. - if `s` is not given, the shape of the input along the axes specified - by `axes` is used. See notes for issue on `ifft` zero padding. - axes : sequence of ints, optional - Axes over which to compute the IFFT. If not given, the last ``len(s)`` - axes are used, or all axes if `s` is also not specified. - Repeated indices in `axes` means that the inverse transform over that - axis is performed multiple times. - norm : {"backward", "ortho", "forward"}, optional - .. versionadded:: 1.10.0 - - Normalization mode (see `numpy.fft`). Default is "backward". - Indicates which direction of the forward/backward pair of transforms - is scaled and with what normalization factor. - - .. versionadded:: 1.20.0 - - The "backward", "forward" values were added. - - Returns - ------- - out : complex ndarray - The truncated or zero-padded input, transformed along the axes - indicated by `axes`, or by a combination of `s` or `a`, - as explained in the parameters section above. - - Raises - ------ - ValueError - If `s` and `axes` have different length. - IndexError - If an element of `axes` is larger than than the number of axes of `a`. - - See Also - -------- - numpy.fft : Overall view of discrete Fourier transforms, with definitions - and conventions used. - fftn : The forward *n*-dimensional FFT, of which `ifftn` is the inverse. - ifft : The one-dimensional inverse FFT. - ifft2 : The two-dimensional inverse FFT. - ifftshift : Undoes `fftshift`, shifts zero-frequency terms to beginning - of array. - - Notes - ----- - See `numpy.fft` for definitions and conventions used. - - Zero-padding, analogously with `ifft`, is performed by appending zeros to - the input along the specified dimension. Although this is the common - approach, it might lead to surprising results. If another form of zero - padding is desired, it must be performed before `ifftn` is called. - - Examples - -------- - >>> a = np.eye(4) - >>> np.fft.ifftn(np.fft.fftn(a, axes=(0,)), axes=(1,)) - array([[ 1.+0.j, 0.+0.j, 0.+0.j, 0.+0.j], - [ 0.+0.j, 1.+0.j, 0.+0.j, 0.+0.j], - [ 0.+0.j, 0.+0.j, 1.+0.j, 0.+0.j], - [ 0.+0.j, 0.+0.j, 0.+0.j, 1.+0.j]]) - - - Create and plot an image with band-limited frequency content: - - >>> import matplotlib.pyplot as plt - >>> n = np.zeros((200,200), dtype=complex) - >>> n[60:80, 20:40] = np.exp(1j*np.random.uniform(0, 2*np.pi, (20, 20))) - >>> im = np.fft.ifftn(n).real - >>> plt.imshow(im) - - >>> plt.show() + For full documentation refer to `numpy.fft.ifftn`. """ - x = _downcast_float128_array(a) s, axes = _cook_nd_args(x, s, axes) fsc = _compute_fwd_scale(norm, s, x.shape) return trycall( - mkl_fft.ifftn, (x,), {"s": s, "axes": axes, "fwd_scale": fsc} + mkl_fft.ifftn, + (x,), + {"s": s, "axes": axes, "fwd_scale": fsc, "out": out}, ) -def fft2(a, s=None, axes=(-2, -1), norm=None): +def fft2(a, s=None, axes=(-2, -1), norm=None, out=None): """ Compute the 2-dimensional discrete Fourier Transform - This function computes the *n*-dimensional discrete Fourier Transform - over any axes in an *M*-dimensional array by means of the - Fast Fourier Transform (FFT). By default, the transform is computed over - the last two axes of the input array, i.e., a 2-dimensional FFT. - - Parameters - ---------- - a : array_like - Input array, can be complex - s : sequence of ints, optional - Shape (length of each transformed axis) of the output - (`s[0]` refers to axis 0, `s[1]` to axis 1, etc.). - This corresponds to `n` for `fft(x, n)`. - Along each axis, if the given shape is smaller than that of the input, - the input is cropped. If it is larger, the input is padded with zeros. - if `s` is not given, the shape of the input along the axes specified - by `axes` is used. - axes : sequence of ints, optional - Axes over which to compute the FFT. If not given, the last two - axes are used. A repeated index in `axes` means the transform over - that axis is performed multiple times. A one-element sequence means - that a one-dimensional FFT is performed. - norm : {"backward", "ortho", "forward"}, optional - .. versionadded:: 1.10.0 - - Normalization mode (see `numpy.fft`). Default is "backward". - Indicates which direction of the forward/backward pair of transforms - is scaled and with what normalization factor. - - .. versionadded:: 1.20.0 - - The "backward", "forward" values were added. - - Returns - ------- - out : complex ndarray - The truncated or zero-padded input, transformed along the axes - indicated by `axes`, or the last two axes if `axes` is not given. - - Raises - ------ - ValueError - If `s` and `axes` have different length, or `axes` not given and - ``len(s) != 2``. - IndexError - If an element of `axes` is larger than than the number of axes of `a`. - - See Also - -------- - numpy.fft : Overall view of discrete Fourier transforms, with definitions - and conventions used. - ifft2 : The inverse two-dimensional FFT. - fft : The one-dimensional FFT. - fftn : The *n*-dimensional FFT. - fftshift : Shifts zero-frequency terms to the center of the array. - For two-dimensional input, swaps first and third quadrants, and second - and fourth quadrants. - - Notes - ----- - `fft2` is just `fftn` with a different default for `axes`. - - The output, analogously to `fft`, contains the term for zero frequency in - the low-order corner of the transformed axes, the positive frequency terms - in the first half of these axes, the term for the Nyquist frequency in the - middle of the axes and the negative frequency terms in the second half of - the axes, in order of decreasingly negative frequency. - - See `fftn` for details and a plotting example, and `numpy.fft` for - definitions and conventions used. - - - Examples - -------- - >>> a = np.mgrid[:5, :5][0] - >>> np.fft.fft2(a) - array([[ 50.0 +0.j , 0.0 +0.j , 0.0 +0.j , - 0.0 +0.j , 0.0 +0.j ], - [-12.5+17.20477401j, 0.0 +0.j , 0.0 +0.j , - 0.0 +0.j , 0.0 +0.j ], - [-12.5 +4.0614962j , 0.0 +0.j , 0.0 +0.j , - 0.0 +0.j , 0.0 +0.j ], - [-12.5 -4.0614962j , 0.0 +0.j , 0.0 +0.j , - 0.0 +0.j , 0.0 +0.j ], - [-12.5-17.20477401j, 0.0 +0.j , 0.0 +0.j , - 0.0 +0.j , 0.0 +0.j ]]) + For full documentation refer to `numpy.fft.fft2`. """ - - return fftn(a, s=s, axes=axes, norm=norm) + return fftn(a, s=s, axes=axes, norm=norm, out=out) -def ifft2(a, s=None, axes=(-2, -1), norm=None): +def ifft2(a, s=None, axes=(-2, -1), norm=None, out=None): """ Compute the 2-dimensional inverse discrete Fourier Transform. - This function computes the inverse of the 2-dimensional discrete Fourier - Transform over any number of axes in an M-dimensional array by means of - the Fast Fourier Transform (FFT). In other words, ``ifft2(fft2(a)) == a`` - to within numerical accuracy. By default, the inverse transform is - computed over the last two axes of the input array. - - The input, analogously to `ifft`, should be ordered in the same way as is - returned by `fft2`, i.e. it should have the term for zero frequency - in the low-order corner of the two axes, the positive frequency terms in - the first half of these axes, the term for the Nyquist frequency in the - middle of the axes and the negative frequency terms in the second half of - both axes, in order of decreasingly negative frequency. - - Parameters - ---------- - a : array_like - Input array, can be complex. - s : sequence of ints, optional - Shape (length of each axis) of the output (``s[0]`` refers to axis 0, - ``s[1]`` to axis 1, etc.). This corresponds to `n` for ``ifft(x, n)``. - Along each axis, if the given shape is smaller than that of the input, - the input is cropped. If it is larger, the input is padded with zeros. - if `s` is not given, the shape of the input along the axes specified - by `axes` is used. See notes for issue on `ifft` zero padding. - axes : sequence of ints, optional - Axes over which to compute the FFT. If not given, the last two - axes are used. A repeated index in `axes` means the transform over - that axis is performed multiple times. A one-element sequence means - that a one-dimensional FFT is performed. - norm : {"backward", "ortho", "forward"}, optional - .. versionadded:: 1.10.0 - - Normalization mode (see `numpy.fft`). Default is "backward". - Indicates which direction of the forward/backward pair of transforms - is scaled and with what normalization factor. - - .. versionadded:: 1.20.0 - - The "backward", "forward" values were added. - - Returns - ------- - out : complex ndarray - The truncated or zero-padded input, transformed along the axes - indicated by `axes`, or the last two axes if `axes` is not given. - - Raises - ------ - ValueError - If `s` and `axes` have different length, or `axes` not given and - ``len(s) != 2``. - IndexError - If an element of `axes` is larger than than the number of axes of `a`. - - See Also - -------- - numpy.fft : Overall view of discrete Fourier transforms, with definitions - and conventions used. - fft2 : The forward 2-dimensional FFT, of which `ifft2` is the inverse. - ifftn : The inverse of the *n*-dimensional FFT. - fft : The one-dimensional FFT. - ifft : The one-dimensional inverse FFT. - - Notes - ----- - `ifft2` is just `ifftn` with a different default for `axes`. - - See `ifftn` for details and a plotting example, and `numpy.fft` for - definition and conventions used. - - Zero-padding, analogously with `ifft`, is performed by appending zeros to - the input along the specified dimension. Although this is the common - approach, it might lead to surprising results. If another form of zero - padding is desired, it must be performed before `ifft2` is called. - - Examples - -------- - >>> a = 4 * np.eye(4) - >>> np.fft.ifft2(a) - array([[ 1.+0.j, 0.+0.j, 0.+0.j, 0.+0.j], - [ 0.+0.j, 0.+0.j, 0.+0.j, 1.+0.j], - [ 0.+0.j, 0.+0.j, 1.+0.j, 0.+0.j], - [ 0.+0.j, 1.+0.j, 0.+0.j, 0.+0.j]]) + For full documentation refer to `numpy.fft.ifft2`. """ + return ifftn(a, s=s, axes=axes, norm=norm, out=out) - return ifftn(a, s=s, axes=axes, norm=norm) - -def rfftn(a, s=None, axes=None, norm=None): +def rfftn(a, s=None, axes=None, norm=None, out=None): """ Compute the N-dimensional discrete Fourier Transform for real input. - This function computes the N-dimensional discrete Fourier Transform over - any number of axes in an M-dimensional real array by means of the Fast - Fourier Transform (FFT). By default, all axes are transformed, with the - real transform performed over the last axis, while the remaining - transforms are complex. - - Parameters - ---------- - a : array_like - Input array, taken to be real. - s : sequence of ints, optional - Shape (length along each transformed axis) to use from the input. - (``s[0]`` refers to axis 0, ``s[1]`` to axis 1, etc.). - The final element of `s` corresponds to `n` for ``rfft(x, n)``, while - for the remaining axes, it corresponds to `n` for ``fft(x, n)``. - Along any axis, if the given shape is smaller than that of the input, - the input is cropped. If it is larger, the input is padded with zeros. - if `s` is not given, the shape of the input along the axes specified - by `axes` is used. - axes : sequence of ints, optional - Axes over which to compute the FFT. If not given, the last ``len(s)`` - axes are used, or all axes if `s` is also not specified. - norm : {"backward", "ortho", "forward"}, optional - .. versionadded:: 1.10.0 - - Normalization mode (see `numpy.fft`). Default is "backward". - Indicates which direction of the forward/backward pair of transforms - is scaled and with what normalization factor. - - .. versionadded:: 1.20.0 - - The "backward", "forward" values were added. - - Returns - ------- - out : complex ndarray - The truncated or zero-padded input, transformed along the axes - indicated by `axes`, or by a combination of `s` and `a`, - as explained in the parameters section above. - The length of the last axis transformed will be ``s[-1]//2+1``, - while the remaining transformed axes will have lengths according to - `s`, or unchanged from the input. - - Raises - ------ - ValueError - If `s` and `axes` have different length. - IndexError - If an element of `axes` is larger than than the number of axes of `a`. - - See Also - -------- - irfftn : The inverse of `rfftn`, i.e. the inverse of the n-dimensional FFT - of real input. - fft : The one-dimensional FFT, with definitions and conventions used. - rfft : The one-dimensional FFT of real input. - fftn : The n-dimensional FFT. - rfft2 : The two-dimensional FFT of real input. - - Notes - ----- - The transform for real input is performed over the last transformation - axis, as by `rfft`, then the transform over the remaining axes is - performed as by `fftn`. The order of the output is as for `rfft` for the - final transformation axis, and as for `fftn` for the remaining - transformation axes. - - See `fft` for details, definitions and conventions used. - - Examples - -------- - >>> a = np.ones((2, 2, 2)) - >>> np.fft.rfftn(a) - array([[[ 8.+0.j, 0.+0.j], - [ 0.+0.j, 0.+0.j]], - [[ 0.+0.j, 0.+0.j], - [ 0.+0.j, 0.+0.j]]]) - - >>> np.fft.rfftn(a, axes=(2, 0)) - array([[[ 4.+0.j, 0.+0.j], - [ 4.+0.j, 0.+0.j]], - [[ 0.+0.j, 0.+0.j], - [ 0.+0.j, 0.+0.j]]]) + For full documentation refer to `numpy.fft.rfftn`. """ - x = _downcast_float128_array(a) s, axes = _cook_nd_args(x, s, axes) fsc = _compute_fwd_scale(norm, s, x.shape) return trycall( - mkl_fft.rfftn, (x,), {"s": s, "axes": axes, "fwd_scale": fsc} + mkl_fft.rfftn, + (x,), + {"s": s, "axes": axes, "fwd_scale": fsc, "out": out}, ) -def rfft2(a, s=None, axes=(-2, -1), norm=None): +def rfft2(a, s=None, axes=(-2, -1), norm=None, out=None): """ Compute the 2-dimensional FFT of a real array. - Parameters - ---------- - a : array - Input array, taken to be real. - s : sequence of ints, optional - Shape of the FFT. - axes : sequence of ints, optional - Axes over which to compute the FFT. - norm : {"backward", "ortho", "forward"}, optional - .. versionadded:: 1.10.0 - - Normalization mode (see `numpy.fft`). Default is "backward". - Indicates which direction of the forward/backward pair of transforms - is scaled and with what normalization factor. - - .. versionadded:: 1.20.0 - - The "backward", "forward" values were added. - - Returns - ------- - out : ndarray - The result of the real 2-D FFT. - - See Also - -------- - rfftn : Compute the N-dimensional discrete Fourier Transform for real - input. - - Notes - ----- - This is really just `rfftn` with different default behavior. - For more details see `rfftn`. + For full documentation refer to `numpy.fft.rfft2`. """ + return rfftn(a, s=s, axes=axes, norm=norm, out=out) - return rfftn(a, s, axes, norm) - -def irfftn(a, s=None, axes=None, norm=None): +def irfftn(a, s=None, axes=None, norm=None, out=None): """ - Compute the inverse of the N-dimensional FFT of real input. - - This function computes the inverse of the N-dimensional discrete - Fourier Transform for real input over any number of axes in an - M-dimensional array by means of the Fast Fourier Transform (FFT). In - other words, ``irfftn(rfftn(a), a.shape) == a`` to within numerical - accuracy. (The ``a.shape`` is necessary like ``len(a)`` is for `irfft`, - and for the same reason.) - - The input should be ordered in the same way as is returned by `rfftn`, - i.e. as for `irfft` for the final transformation axis, and as for `ifftn` - along all the other axes. - - Parameters - ---------- - a : array_like - Input array. - s : sequence of ints, optional - Shape (length of each transformed axis) of the output - (``s[0]`` refers to axis 0, ``s[1]`` to axis 1, etc.). `s` is also the - number of input points used along this axis, except for the last axis, - where ``s[-1]//2+1`` points of the input are used. - Along any axis, if the shape indicated by `s` is smaller than that of - the input, the input is cropped. If it is larger, the input is padded - with zeros. If `s` is not given, the shape of the input along the - axes specified by `axes` is used. - axes : sequence of ints, optional - Axes over which to compute the inverse FFT. If not given, the last - `len(s)` axes are used, or all axes if `s` is also not specified. - Repeated indices in `axes` means that the inverse transform over that - axis is performed multiple times. - norm : {"backward", "ortho", "forward"}, optional - .. versionadded:: 1.10.0 - - Normalization mode (see `numpy.fft`). Default is "backward". - Indicates which direction of the forward/backward pair of transforms - is scaled and with what normalization factor. - - .. versionadded:: 1.20.0 - - The "backward", "forward" values were added. - - Returns - ------- - out : ndarray - The truncated or zero-padded input, transformed along the axes - indicated by `axes`, or by a combination of `s` or `a`, - as explained in the parameters section above. - The length of each transformed axis is as given by the corresponding - element of `s`, or the length of the input in every axis except for the - last one if `s` is not given. In the final transformed axis the length - of the output when `s` is not given is ``2*(m-1)`` where ``m`` is the - length of the final transformed axis of the input. To get an odd - number of output points in the final axis, `s` must be specified. - - Raises - ------ - ValueError - If `s` and `axes` have different length. - IndexError - If an element of `axes` is larger than than the number of axes of `a`. - - See Also - -------- - rfftn : The forward n-dimensional FFT of real input, - of which `ifftn` is the inverse. - fft : The one-dimensional FFT, with definitions and conventions used. - irfft : The inverse of the one-dimensional FFT of real input. - irfft2 : The inverse of the two-dimensional FFT of real input. - - Notes - ----- - See `fft` for definitions and conventions used. - - See `rfft` for definitions and conventions used for real input. - - Examples - -------- - >>> a = np.zeros((3, 2, 2)) - >>> a[0, 0, 0] = 3 * 2 * 2 - >>> np.fft.irfftn(a) - array([[[ 1., 1.], - [ 1., 1.]], - [[ 1., 1.], - [ 1., 1.]], - [[ 1., 1.], - [ 1., 1.]]]) + Compute the inverse of `rfftn`. + + For full documentation refer to `numpy.fft.irfftn`. """ @@ -1315,48 +308,17 @@ def irfftn(a, s=None, axes=None, norm=None): fsc = _compute_fwd_scale(norm, s, x.shape) return trycall( - mkl_fft.irfftn, (x,), {"s": s, "axes": axes, "fwd_scale": fsc} + mkl_fft.irfftn, + (x,), + {"s": s, "axes": axes, "fwd_scale": fsc, "out": out}, ) -def irfft2(a, s=None, axes=(-2, -1), norm=None): +def irfft2(a, s=None, axes=(-2, -1), norm=None, out=None): """ - Compute the 2-dimensional inverse FFT of a real array. - - Parameters - ---------- - a : array_like - The input array - s : sequence of ints, optional - Shape of the inverse FFT. - axes : sequence of ints, optional - The axes over which to compute the inverse fft. - Default is the last two axes. - norm : {"backward", "ortho", "forward"}, optional - .. versionadded:: 1.10.0 - - Normalization mode (see `numpy.fft`). Default is "backward". - Indicates which direction of the forward/backward pair of transforms - is scaled and with what normalization factor. - - .. versionadded:: 1.20.0 - - The "backward", "forward" values were added. - - Returns - ------- - out : ndarray - The result of the inverse real 2-D FFT. - - See Also - -------- - irfftn : Compute the inverse of the N-dimensional FFT of real input. - - Notes - ----- - This is really `irfftn` with different defaults. - For more details see `irfftn`. + Compute the inverse FFT of `rfft2`. - """ + For full documentation refer to `numpy.fft.irfft2`. - return irfftn(a, s, axes, norm) + """ + return irfftn(a, s=s, axes=axes, norm=norm, out=out) diff --git a/mkl_fft/_pydfti.pyx b/mkl_fft/_pydfti.pyx index 3f81f6e..82cebee 100644 --- a/mkl_fft/_pydfti.pyx +++ b/mkl_fft/_pydfti.pyx @@ -135,11 +135,6 @@ cdef extern from "src/mklfft.h": cnp.ndarray, int, int, cnp.ndarray, double, DftiCache* ) - int double_mkl_rfft_in(cnp.ndarray, int, int, double, DftiCache*) - int double_mkl_irfft_in(cnp.ndarray, int, int, double, DftiCache*) - int float_mkl_rfft_in(cnp.ndarray, int, int, double, DftiCache*) - int float_mkl_irfft_in(cnp.ndarray, int, int, double, DftiCache*) - int cdouble_double_mkl_irfft_out( cnp.ndarray, int, int, cnp.ndarray, double, DftiCache* ) @@ -185,15 +180,27 @@ cdef int _datacopied(cnp.ndarray arr, object orig): return 1 if (arr_obj.base is None) else 0 -def fft(x, n=None, axis=-1, overwrite_x=False, fwd_scale=1.0): +def fft(x, n=None, axis=-1, overwrite_x=False, fwd_scale=1.0, out=None): return _fft1d_impl( - x, n=n, axis=axis, overwrite_x=overwrite_x, direction=+1, fsc=fwd_scale + x, + n=n, + axis=axis, + overwrite_x=overwrite_x, + direction=+1, + fsc=fwd_scale, + out=out, ) -def ifft(x, n=None, axis=-1, overwrite_x=False, fwd_scale=1.0): +def ifft(x, n=None, axis=-1, overwrite_x=False, fwd_scale=1.0, out=None): return _fft1d_impl( - x, n=n, axis=axis, overwrite_x=overwrite_x, direction=-1, fsc=fwd_scale + x, + n=n, + axis=axis, + overwrite_x=overwrite_x, + direction=-1, + fsc=fwd_scale, + out=out, ) @@ -345,12 +352,48 @@ cdef cnp.ndarray __allocate_result( return f_arr +def _get_element_strides(array): + """Convert byte strides to element strides.""" + + byte_strides = array.strides + array_itemsize = array.itemsize + return tuple(s // array_itemsize for s in byte_strides) + + +def _validate_out_array(out, x, out_dtype, axis=None, n=None): + """Validate out keyword argument.""" + + if type(out) is not np.ndarray: + raise TypeError("return array must be of ArrayType") + + x_shape = list(x.shape) + if axis is not None: + x_shape[axis] = n + if out.shape != tuple(x_shape): + raise ValueError( + "output array has wrong shape, expected (%s) got (%s)." + % (tuple(x_shape), out.shape) + ) + + if out.dtype != out_dtype: + raise TypeError( + "Cannot cast 'fft' output from dtype(%s) to dtype(%s)." + % (out_dtype, out.dtype) + ) + + # this routine implements complex forward/backward FFT -# Float/double inputs are not cast to complex, but are effectively +# float/double inputs are not cast to complex, but are effectively # treated as complexes with zero imaginary parts. # All other types are cast to complex double. def _fft1d_impl( - x, n=None, axis=-1, overwrite_x=False, direction=+1, double fsc=1.0 + x, + n=None, + axis=-1, + overwrite_x=False, + direction=+1, + double fsc=1.0, + out=None, ): """ Uses MKL to perform 1D FFT on the input array x along the given axis. @@ -370,7 +413,9 @@ def _fft1d_impl( x_type = cnp.PyArray_TYPE(x_arr) - if x_type is cnp.NPY_CFLOAT or x_type is cnp.NPY_CDOUBLE: + if out is not None: + in_place = 0 + elif x_type is cnp.NPY_CFLOAT or x_type is cnp.NPY_CDOUBLE: # we can operate in place if requested. if in_place: if not cnp.PyArray_ISONESEGMENT(x_arr): @@ -436,7 +481,19 @@ def _fft1d_impl( f_type = cnp.NPY_CFLOAT else: f_type = cnp.NPY_CDOUBLE - f_arr = __allocate_result(x_arr, n_, axis_, f_type) + + if out is None: + f_arr = __allocate_result(x_arr, n_, axis_, f_type) + else: + out_dtype = np.dtype(cnp.PyArray_DescrFromType(f_type)) + _validate_out_array(out, x, out_dtype, axis=axis_, n=n_) + # out array that is used in OneMKL c2c FFT must have the exact same + # stride as input array. If not, we need to allocate a new array. + # TODO: check to see if this condition can be relaxed + if _get_element_strides(x) == _get_element_strides(out): + f_arr = out + else: + f_arr = __allocate_result(x_arr, n_, axis_, f_type) # call out-of-place FFT _cache_capsule = _tls_dfti_cache_capsule() @@ -511,7 +568,11 @@ def _fft1d_impl( py_error_msg = c_error_msg raise ValueError("Internal error occurred: {}".format(py_error_msg)) - return f_arr + if out is not None and f_arr is not out: + out[...] = f_arr + return out + else: + return f_arr def rfftpack(x, n=None, axis=-1, overwrite_x=False, fwd_scale=1.0): @@ -784,7 +845,9 @@ def _rr_ifft1d_impl2(x, n=None, axis=-1, overwrite_x=False, double fsc=1.0): # this routine is functionally equivalent to numpy.fft.rfft -def _rc_fft1d_impl(x, n=None, axis=-1, overwrite_x=False, double fsc=1.0): +def _rc_fft1d_impl( + x, n=None, axis=-1, overwrite_x=False, double fsc=1.0, out=None +): """ Uses MKL to perform 1D FFT on the real input array x along the given axis, producing complex output, but giving only half of the harmonics. @@ -831,7 +894,24 @@ def _rc_fft1d_impl(x, n=None, axis=-1, overwrite_x=False, double fsc=1.0): # it can be done only if 2*(n_ // 2 + 1) <= x_arr.shape[axis_] which is not # the common usage f_type = cnp.NPY_CFLOAT if x_type is cnp.NPY_FLOAT else cnp.NPY_CDOUBLE - f_arr = __allocate_result(x_arr, n_ // 2 + 1, axis_, f_type) + f_shape = n_ // 2 + 1 + if out is None: + f_arr = __allocate_result(x_arr, f_shape, axis_, f_type) + else: + out_dtype = np.dtype(cnp.PyArray_DescrFromType(f_type)) + _validate_out_array(out, x, out_dtype, axis=axis_, n=f_shape) + # out array that is used in OneMKL r2c FFT must have comparable strides + # with input array. If not, we need to allocate a new array. + # For r2c, out and input arrays have different size and strides cannot + # be compared directly. + # TODO: currently instead of this condition, we check both input + # and output to be c_contig or f_contig, relax this condition + c_contig = x.flags.c_contiguous and out.flags.c_contiguous + f_contig = x.flags.f_contiguous and out.flags.f_contiguous + if c_contig or f_contig: + f_arr = out + else: + f_arr = __allocate_result(x_arr, f_shape, axis_, f_type) # call out-of-place FFT if x_type is cnp.NPY_FLOAT: @@ -858,7 +938,11 @@ def _rc_fft1d_impl(x, n=None, axis=-1, overwrite_x=False, double fsc=1.0): "Internal error occurred: {}".format(str(py_error_msg)) ) - return f_arr + if out is not None and f_arr is not out: + out[...] = f_arr + return out + else: + return f_arr cdef int _is_integral(object num): @@ -876,7 +960,9 @@ cdef int _is_integral(object num): # this routine is functionally equivalent to numpy.fft.irfft -def _rc_ifft1d_impl(x, n=None, axis=-1, overwrite_x=False, double fsc=1.0): +def _rc_ifft1d_impl( + x, n=None, axis=-1, overwrite_x=False, double fsc=1.0, out=None +): """ Uses MKL to perform 1D FFT on the real input array x along the given axis, producing complex output, but giving only half of the harmonics. @@ -929,7 +1015,23 @@ def _rc_ifft1d_impl(x, n=None, axis=-1, overwrite_x=False, double fsc=1.0): pass else: f_type = cnp.NPY_FLOAT if x_type is cnp.NPY_CFLOAT else cnp.NPY_DOUBLE - f_arr = __allocate_result(x_arr, n_, axis_, f_type) + if out is None: + f_arr = __allocate_result(x_arr, n_, axis_, f_type) + else: + out_dtype = np.dtype(cnp.PyArray_DescrFromType(f_type)) + _validate_out_array(out, x, out_dtype, axis=axis_, n=n_) + # out array that is used in OneMKL c2r FFT must have comparable + # strides with input array. If not, we need to allocate a new + # array. For c2r, out and input arrays have different size and + # strides cannot be compared directly. + # TODO: currently instead of this condition, we check both input + # and output to be c_contig or f_contig, relax this condition + c_contig = x.flags.c_contiguous and out.flags.c_contiguous + f_contig = x.flags.f_contiguous and out.flags.f_contiguous + if c_contig or f_contig: + f_arr = out + else: + f_arr = __allocate_result(x_arr, n_, axis_, f_type) # call out-of-place FFT if x_type is cnp.NPY_CFLOAT: @@ -956,15 +1058,19 @@ def _rc_ifft1d_impl(x, n=None, axis=-1, overwrite_x=False, double fsc=1.0): "Internal error occurred: {}".format(str(py_error_msg)) ) - return f_arr + if out is not None and f_arr is not out: + out[...] = f_arr + return out + else: + return f_arr -def rfft(x, n=None, axis=-1, fwd_scale=1.0): - return _rc_fft1d_impl(x, n=n, axis=axis, fsc=fwd_scale) +def rfft(x, n=None, axis=-1, fwd_scale=1.0, out=None): + return _rc_fft1d_impl(x, n=n, axis=axis, fsc=fwd_scale, out=out) -def irfft(x, n=None, axis=-1, fwd_scale=1.0): - return _rc_ifft1d_impl(x, n=n, axis=axis, fsc=fwd_scale) +def irfft(x, n=None, axis=-1, fwd_scale=1.0, out=None): + return _rc_ifft1d_impl(x, n=n, axis=axis, fsc=fwd_scale, out=out) # ============================== ND ====================================== # @@ -1049,9 +1155,9 @@ def _init_nd_shape_and_axes(x, shape, axes): return shape, axes -def _cook_nd_args(a, s=None, axes=None, invreal=0): +def _cook_nd_args(a, s=None, axes=None, invreal=False): if s is None: - shapeless = 1 + shapeless = True if axes is None: s = list(a.shape) else: @@ -1062,7 +1168,7 @@ def _cook_nd_args(a, s=None, axes=None, invreal=0): s = range(len(axes) + 1) pass else: - shapeless = 0 + shapeless = False s = list(s) if axes is None: axes = list(range(-len(s), 0)) @@ -1081,6 +1187,7 @@ def _iter_fftnd( overwrite_x=False, scale_function=lambda n, ind: 1.0, + out=None, ): a = np.asarray(a) s, axes = _init_nd_shape_and_axes(a, s, axes) @@ -1092,6 +1199,7 @@ def _iter_fftnd( axis = axes[ii], overwrite_x=ovwr, fwd_scale=scale_function(s[ii], ii), + out=out, ) ovwr = True return a @@ -1112,7 +1220,8 @@ def flat_to_multi(ind, shape): def iter_complementary(x, axes, func, kwargs, result): if axes is None: - return func(x, **kwargs) + # s and axes are None, direct N-D FFT + return func(x, **kwargs, out=result) x_shape = x.shape nd = x.ndim r = list(range(nd)) @@ -1134,12 +1243,28 @@ def iter_complementary(x, axes, func, kwargs, result): m_ind = flat_to_multi(ind, sub_shape) for k1, k2 in zip(dual_ind, m_ind): sl[k1] = k2 - np.copyto(result[tuple(sl)], func(x[tuple(sl)], **kwargs)) + if np.issubdtype(x.dtype, np.complexfloating): + func(x[tuple(sl)], **kwargs, out=result[tuple(sl)]) + else: + # For c2c FFT, if the input is real, half of the output is the + # complex conjugate of the other half. Instead of upcasting the + # input to complex and performing c2c FFT, we perform rfft and then + # construct the other half using the first half. + # However, when using the `out` keyword here I encountered a result + # in tests/third_party/scipy/test_basic.py::test_fft_with_order + # that was correct but the elements were not necessarily similar to + # NumPy. For example, an element in the first half of mkl_fft output + # array appeared in the second half of the NumPy output array, + # while the equivalent element in the NumPy array was the conjugate + # of the mkl_fft output array. + np.copyto(result[tuple(sl)], func(x[tuple(sl)], **kwargs)) return result -def _direct_fftnd(x, overwrite_x=False, direction=+1, double fsc=1.0): +def _direct_fftnd( + x, overwrite_x=False, direction=+1, double fsc=1.0, out=None +): """Perform n-dimensional FFT over all axes""" cdef int err cdef cnp.ndarray x_arr "xxnd_arrayObject" @@ -1180,6 +1305,9 @@ def _direct_fftnd(x, overwrite_x=False, direction=+1, double fsc=1.0): assert x_type == cnp.NPY_CDOUBLE in_place = 1 + if out is not None: + in_place = 0 + if in_place: if x_type == cnp.NPY_CDOUBLE or x_type == cnp.NPY_CFLOAT: in_place = 1 @@ -1206,7 +1334,19 @@ def _direct_fftnd(x, overwrite_x=False, direction=+1, double fsc=1.0): f_type = cnp.NPY_CDOUBLE else: f_type = cnp.NPY_CFLOAT - f_arr = __allocate_result(x_arr, -1, 0, f_type) + if out is None: + f_arr = __allocate_result(x_arr, -1, 0, f_type) + else: + out_dtype = np.dtype(cnp.PyArray_DescrFromType(f_type)) + _validate_out_array(out, x, out_dtype) + # out array that is used in OneMKL c2c FFT must have the exact same + # stride as input array. If not, we need to allocate a new array. + # TODO: check to see if this condition can be relaxed + if _get_element_strides(x) == _get_element_strides(out): + f_arr = out + else: + f_arr = __allocate_result(x_arr, -1, 0, f_type) + if x_type == cnp.NPY_CDOUBLE: if dir_ == 1: err = cdouble_cdouble_mkl_fftnd_out(x_arr, f_arr, fsc) @@ -1230,15 +1370,21 @@ def _direct_fftnd(x, overwrite_x=False, direction=+1, double fsc=1.0): else: raise ValueError("An input argument x is not complex type array") - return f_arr + if out is not None and f_arr is not out: + out[...] = f_arr + return out + else: + return f_arr def _check_shapes_for_direct(xs, shape, axes): if len(axes) > 7: # Intel MKL supports up to 7D return False if not (len(xs) == len(shape)): + # full-dimensional transform return False if not (len(set(axes)) == len(axes)): + # repeated axes return False for xsi, ai in zip(xs, axes): try: @@ -1260,11 +1406,18 @@ def _output_dtype(dt): def _fftnd_impl( - x, s=None, axes=None, overwrite_x=False, direction=+1, double fsc=1.0 + x, + s=None, + axes=None, + overwrite_x=False, + direction=+1, + double fsc=1.0, + out=None, ): if direction not in [-1, +1]: raise ValueError("Direction of FFT should +1 or -1") + valid_dtypes = [np.complex64, np.complex128, np.float32, np.float64] # _direct_fftnd requires complex type, and full-dimensional transform if isinstance(x, np.ndarray) and x.size != 0 and x.ndim > 1: _direct = s is None and axes is None @@ -1274,24 +1427,27 @@ def _fftnd_impl( xs, xa = _cook_nd_args(x, s, axes) if _check_shapes_for_direct(xs, x.shape, xa): _direct = True - _direct = ( - _direct - and x.dtype in [np.complex64, np.complex128, np.float32, np.float64] - ) + _direct = _direct and x.dtype in valid_dtypes else: _direct = False if _direct: return _direct_fftnd( - x, overwrite_x=overwrite_x, direction=direction, fsc=fsc + x, + overwrite_x=overwrite_x, + direction=direction, + fsc=fsc, + out=out, ) else: - if ( - s is None - and x.dtype in [np.csingle, np.cdouble, np.single, np.double] - ): + if s is None and x.dtype in valid_dtypes: x = np.asarray(x) - res = np.empty(x.shape, dtype=_output_dtype(x.dtype)) + if out is None: + res = np.empty_like(x, dtype=_output_dtype(x.dtype)) + else: + _validate_out_array(out, x, _output_dtype(x.dtype)) + res = out + return iter_complementary( x, axes, _direct_fftnd, @@ -1303,43 +1459,69 @@ def _fftnd_impl( res ) else: + # perform N-D FFT as a series of 1D FFTs sc = fsc return _iter_fftnd(x, s=s, axes=axes, overwrite_x=overwrite_x, scale_function=lambda n, i: sc if i == 0 else 1., - function=fft if direction == 1 else ifft) + function=fft if direction == 1 else ifft, + out=out) -def fft2(x, s=None, axes=(-2, -1), overwrite_x=False, fwd_scale=1.0): +def fft2(x, s=None, axes=(-2, -1), overwrite_x=False, fwd_scale=1.0, out=None): return _fftnd_impl( - x, s=s, axes=axes, overwrite_x=overwrite_x, direction=+1, fsc=fwd_scale + x, + s=s, + axes=axes, + overwrite_x=overwrite_x, + direction=+1, + fsc=fwd_scale, + out=out, ) -def ifft2(x, s=None, axes=(-2, -1), overwrite_x=False, fwd_scale=1.0): +def ifft2(x, s=None, axes=(-2, -1), overwrite_x=False, fwd_scale=1.0, out=None): return _fftnd_impl( - x, s=s, axes=axes, overwrite_x=overwrite_x, direction=-1, fsc=fwd_scale + x, + s=s, + axes=axes, + overwrite_x=overwrite_x, + direction=-1, + fsc=fwd_scale, + out=out, ) -def fftn(x, s=None, axes=None, overwrite_x=False, fwd_scale=1.0): +def fftn(x, s=None, axes=None, overwrite_x=False, fwd_scale=1.0, out=None): return _fftnd_impl( - x, s=s, axes=axes, overwrite_x=overwrite_x, direction=+1, fsc=fwd_scale + x, + s=s, + axes=axes, + overwrite_x=overwrite_x, + direction=+1, + fsc=fwd_scale, + out=out, ) -def ifftn(x, s=None, axes=None, overwrite_x=False, fwd_scale=1.0): +def ifftn(x, s=None, axes=None, overwrite_x=False, fwd_scale=1.0, out=None): return _fftnd_impl( - x, s=s, axes=axes, overwrite_x=overwrite_x, direction=-1, fsc=fwd_scale + x, + s=s, + axes=axes, + overwrite_x=overwrite_x, + direction=-1, + fsc=fwd_scale, + out=out, ) -def rfft2(x, s=None, axes=(-2, -1), fwd_scale=1.0): - return rfftn(x, s=s, axes=axes, fwd_scale=fwd_scale) +def rfft2(x, s=None, axes=(-2, -1), fwd_scale=1.0, out=None): + return rfftn(x, s=s, axes=axes, fwd_scale=fwd_scale, out=out) -def irfft2(x, s=None, axes=(-2, -1), fwd_scale=1.0): - return irfftn(x, s=s, axes=axes, fwd_scale=fwd_scale) +def irfft2(x, s=None, axes=(-2, -1), fwd_scale=1.0, out=None): + return irfftn(x, s=s, axes=axes, fwd_scale=fwd_scale, out=out) def _remove_axis(s, axes, axis_to_remove): @@ -1394,7 +1576,7 @@ def _fix_dimensions(cnp.ndarray arr, object s, object axes): return np.pad(arr, tuple(pad_widths), "constant") -def rfftn(x, s=None, axes=None, fwd_scale=1.0): +def rfftn(x, s=None, axes=None, fwd_scale=1.0, out=None): a = np.asarray(x) no_trim = (s is None) and (axes is None) s, axes = _cook_nd_args(a, s, axes) @@ -1402,7 +1584,8 @@ def rfftn(x, s=None, axes=None, fwd_scale=1.0): # trim array, so that rfft avoids doing unnecessary computations if not no_trim: a = _trim_array(a, s, axes) - a = rfft(a, n = s[-1], axis=la, fwd_scale=fwd_scale) + # r2c along last axis + a = rfft(a, n = s[-1], axis=la, fwd_scale=fwd_scale, out=out) if len(s) > 1: if not no_trim: ss = list(s) @@ -1410,24 +1593,27 @@ def rfftn(x, s=None, axes=None, fwd_scale=1.0): a = _fix_dimensions(a, tuple(ss), axes) len_axes = len(axes) if len(set(axes)) == len_axes and len_axes == a.ndim and len_axes > 2: + # a series of ND c2c FFTs along last axis ss, aa = _remove_axis(s, axes, -1) ind = [slice(None, None, 1),] * len(s) for ii in range(a.shape[la]): ind[la] = ii tind = tuple(ind) a_inp = a[tind] + res = out[tind] if out is not None else None a_res = _fftnd_impl( a_inp, s=ss, axes=aa, - overwrite_x=True, direction=1) + overwrite_x=True, direction=1, out=res) if a_res is not a_inp: a[tind] = a_res # copy in place else: + # a series of 1D c2c FFTs along all axes except last for ii in range(len(axes) - 2, -1, -1): a = fft(a, s[ii], axes[ii], overwrite_x=True) return a -def irfftn(x, s=None, axes=None, fwd_scale=1.0): +def irfftn(x, s=None, axes=None, fwd_scale=1.0, out=None): a = np.asarray(x) no_trim = (s is None) and (axes is None) s, axes = _cook_nd_args(a, s, axes, invreal=True) @@ -1440,6 +1626,7 @@ def irfftn(x, s=None, axes=None, fwd_scale=1.0): ovr_x = True if _datacopied( a, x) else False len_axes = len(axes) if len(set(axes)) == len_axes and len_axes == a.ndim and len_axes > 2: + # a series of ND c2c FFTs along last axis # due to need to write into a, we must copy if not ovr_x: a = a.copy() @@ -1457,14 +1644,18 @@ def irfftn(x, s=None, axes=None, fwd_scale=1.0): ind[la] = ii tind = tuple(ind) a_inp = a[tind] + # out has real dtype and cannot be used in intermediate steps a_res = _fftnd_impl( a_inp, s=ss, axes=aa, overwrite_x=True, direction=-1) if a_res is not a_inp: a[tind] = a_res # copy in place else: + # a series of 1D c2c FFTs along all axes except last for ii in range(len(axes)-1): + # out has real dtype and cannot be used in intermediate steps a = ifft(a, s[ii], axes[ii], overwrite_x=ovr_x) ovr_x = True - a = irfft(a, n = s[-1], axis=la, fwd_scale=fwd_scale) + # c2r along last axis + a = irfft(a, n = s[-1], axis=la, fwd_scale=fwd_scale, out=out) return a diff --git a/mkl_fft/_scipy_fft.py b/mkl_fft/_scipy_fft.py index 36c165e..7a5f658 100644 --- a/mkl_fft/_scipy_fft.py +++ b/mkl_fft/_scipy_fft.py @@ -25,9 +25,8 @@ # OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. """ -An interface for FFT module of Scipy (`scipy.fft`) that uses OneMKL FFT +An interface for FFT module of SciPy (`scipy.fft`) that uses OneMKL FFT in the backend. - """ import contextlib diff --git a/mkl_fft/tests/__init__.py b/mkl_fft/tests/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/mkl_fft/tests/helper.py b/mkl_fft/tests/helper.py new file mode 100644 index 0000000..6a59664 --- /dev/null +++ b/mkl_fft/tests/helper.py @@ -0,0 +1,34 @@ +#!/usr/bin/env python +# Copyright (c) 2025, Intel Corporation +# +# Redistribution and use in source and binary forms, with or without +# modification, are permitted provided that the following conditions are met: +# +# * Redistributions of source code must retain the above copyright notice, +# this list of conditions and the following disclaimer. +# * Redistributions in binary form must reproduce the above copyright +# notice, this list of conditions and the following disclaimer in the +# documentation and/or other materials provided with the distribution. +# * Neither the name of Intel Corporation nor the names of its contributors +# may be used to endorse or promote products derived from this software +# without specific prior written permission. +# +# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS """AS IS" +# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE +# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE +# FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL +# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER +# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, +# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + + +import numpy as np +import pytest + +requires_numpy_2 = pytest.mark.skipif( + np.lib.NumpyVersion(np.__version__) < "2.0.0", + reason="Requires NumPy >= 2.0.0", +) diff --git a/mkl_fft/tests/test_fft1d.py b/mkl_fft/tests/test_fft1d.py index 1e6fcc9..01631fe 100644 --- a/mkl_fft/tests/test_fft1d.py +++ b/mkl_fft/tests/test_fft1d.py @@ -25,11 +25,14 @@ # OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. import numpy as np +import pytest from numpy import random as rnd from numpy.testing import TestCase, assert_, assert_allclose import mkl_fft +from .helper import requires_numpy_2 + def naive_fft1d(vec): L = len(vec) @@ -396,3 +399,61 @@ def test5(self): f1 = mkl_fft.irfftpack(x, axis=a, overwrite_x=ovwr_x) f2 = mkl_fft.rfftpack(f1, axis=a, overwrite_x=ovwr_x) assert_allclose(f2, self.t3.astype(dt), atol=atol) + + +@requires_numpy_2 +@pytest.mark.parametrize("axis", [0, 1, 2]) +@pytest.mark.parametrize("func", ["fft", "ifft"]) +def test_fft_out_strided(axis, func): + shape = (20, 33, 54) + x = rnd.random(shape) + 1j * rnd.random(shape) + out = np.empty(shape, dtype=x.dtype) + + x = x[::2, ::3, ::4] + out = np.empty(x.shape, dtype=x.dtype) + result = getattr(mkl_fft, func)(x, axis=axis, out=out) + expected = getattr(np.fft, func)(x, axis=axis, out=out) + + assert_allclose(result, expected) + + +@requires_numpy_2 +@pytest.mark.parametrize("axis", [0, 1, 2]) +def test_rfft_out_strided(axis): + shape = (20, 33, 54) + x = rnd.random(shape) + if axis == 0: + out_sh = (12, 33, 54) + elif axis == 1: + out_sh = (20, 18, 54) + else: # axis == 2 + out_sh = (20, 33, 32) + out = np.empty(out_sh, dtype=np.complex128) + + x = x[::2, ::3, ::4] + out = out[::2, ::3, ::4] + result = mkl_fft.rfft(x, axis=axis, out=out) + expected = np.fft.rfft(x, axis=axis, out=out) + + assert_allclose(result, expected) + + +@requires_numpy_2 +@pytest.mark.parametrize("axis", [0, 1, 2]) +def test_irfft_out_strided(axis): + shape = (20, 33, 54) + x = rnd.random(shape) + 1j * rnd.random(shape) + if axis == 0: + out_sh = (36, 33, 54) + elif axis == 1: + out_sh = (20, 60, 54) + else: # axis == 2 + out_sh = (20, 33, 104) + out = np.empty(out_sh, dtype=np.float64) + + x = x[::2, ::3, ::4] + out = out[::2, ::3, ::4] + result = mkl_fft.irfft(x, axis=axis, out=out) + expected = np.fft.irfft(x, axis=axis, out=out) + + assert_allclose(result, expected) diff --git a/mkl_fft/tests/test_fftnd.py b/mkl_fft/tests/test_fftnd.py index 9f32c95..f5a744d 100644 --- a/mkl_fft/tests/test_fftnd.py +++ b/mkl_fft/tests/test_fftnd.py @@ -31,6 +31,8 @@ import mkl_fft +from .helper import requires_numpy_2 + reps_64 = (2**11) * np.finfo(np.float64).eps reps_32 = (2**11) * np.finfo(np.float32).eps atol_64 = (2**9) * np.finfo(np.float64).eps @@ -322,3 +324,19 @@ def test_repeated_axes(dtype, axes, func): rtol, atol = _get_rtol_atol(x) assert_allclose(r1, r2, rtol=rtol, atol=atol) + + +@requires_numpy_2 +@pytest.mark.parametrize("axes", [None, (0, 1), (0, 2), (1, 2)]) +@pytest.mark.parametrize("func", ["fftn", "ifftn"]) +def test_out_strided(axes, func): + shape = (20, 30, 40) + x = rnd.random(shape) + 1j * rnd.random(shape) + out = np.empty(shape, dtype=x.dtype) + + x = x[::2, ::3, ::4] + out = out[::2, ::3, ::4] + result = getattr(mkl_fft, func)(x, axes=axes, out=out) + expected = getattr(np.fft, func)(x, axes=axes, out=out) + + assert_allclose(result, expected) diff --git a/mkl_fft/tests/third_party/numpy/test_fft.py b/mkl_fft/tests/third_party/numpy/test_fft.py index 6d5659e..81b7a5f 100644 --- a/mkl_fft/tests/third_party/numpy/test_fft.py +++ b/mkl_fft/tests/third_party/numpy/test_fft.py @@ -1,9 +1,6 @@ # This file includes tests from numpy.fft module: # https://github.com/numpy/numpy/blob/main/numpy/fft/tests/test_pocketfft.py -# TODO: remove pylint disable when out kwarg is added -# pylint: disable=unexpected-keyword-arg - import queue import threading @@ -18,10 +15,7 @@ ) import mkl_fft.interfaces.numpy_fft as mkl_fft - -requires_numpy_2 = pytest.mark.skipif( - np.__version__ < "2.0", reason="Requires NumPy >= 2.0" -) +from mkl_fft.tests.helper import requires_numpy_2 def fft1(x): @@ -83,7 +77,7 @@ def test_identity_long_short_reversed(self, dtype): for i in range(1, maxlen * 2): check_via_c = mkl_fft.fft(mkl_fft.ifft(x, n=i), n=i) assert check_via_c.dtype == x.dtype - assert_allclose(check_via_c, xx[0:i], atol=atol, rtol=0) + assert_allclose(check_via_c, xx[0:i], atol=atol, rtol=1e-14) # For irfft, we can neither recover the imaginary part of # the first element, nor the imaginary part of the last # element if npts is even. So, set to 0 for the comparison. @@ -108,7 +102,7 @@ def test_fft(self): fft1(x) / 30.0, mkl_fft.fft(x, norm="forward"), atol=1e-6 ) - @pytest.mark.skip("out is not supported") + @requires_numpy_2 @pytest.mark.parametrize("axis", (0, 1)) @pytest.mark.parametrize("dtype", (complex, float)) @pytest.mark.parametrize("transpose", (True, False)) @@ -139,7 +133,7 @@ def zeros_like(x): assert result2 is out2 assert_array_equal(result2, expected2) - @pytest.mark.skip("out is not supported") + @requires_numpy_2 @pytest.mark.parametrize("axis", [0, 1]) def test_fft_inplace_out(self, axis): # Test some weirder in-place combinations @@ -200,7 +194,7 @@ def test_fft_inplace_out(self, axis): assert result6 is out6 assert_array_equal(result6, expected1) - @pytest.mark.skip("out is not supported") + @requires_numpy_2 def test_fft_bad_out(self): x = np.arange(30.0) with pytest.raises(TypeError, match="must be of ArrayType"): @@ -566,7 +560,7 @@ def test_all_1d_norm_preserving(self): tmp = back(tmp, n=n, norm=norm) assert_allclose(x_norm, np.linalg.norm(tmp), atol=1e-6) - @pytest.mark.skip("out is not supported") + @requires_numpy_2 @pytest.mark.parametrize("axes", [(0, 1), (0, 2), None]) @pytest.mark.parametrize("dtype", (complex, float)) @pytest.mark.parametrize("transpose", (True, False)) @@ -597,7 +591,7 @@ def zeros_like(x): assert result2 is out2 assert_array_equal(result2, expected2) - @pytest.mark.skip("out is not supported") + @requires_numpy_2 @pytest.mark.parametrize( "fft", [mkl_fft.fftn, mkl_fft.ifftn, mkl_fft.rfftn] ) @@ -617,7 +611,7 @@ def test_fftn_out_and_s_interaction(self, fft): assert result is out assert_array_equal(result, expected) - @pytest.mark.skip("out is not supported") + @requires_numpy_2 @pytest.mark.parametrize("s", [(9, 5, 5), (3, 3, 3)]) def test_irfftn_out_and_s_interaction(self, s): # Since for irfftn, the output is real and thus cannot be used for diff --git a/setup.py b/setup.py index 667fc88..2125391 100644 --- a/setup.py +++ b/setup.py @@ -1,4 +1,3 @@ -#!/usr/bin/env python # Copyright (c) 2017, Intel Corporation # # Redistribution and use in source and binary forms, with or without @@ -29,7 +28,7 @@ from os.path import join import Cython.Build -import numpy as np +import numpy from setuptools import Extension, setup sys.path.insert(0, os.path.dirname(__file__)) # Ensures local imports work @@ -72,7 +71,7 @@ def extensions(): join("mkl_fft", "src", "mklfft.h"), join("mkl_fft", "src", "multi_iter.h"), ], - include_dirs=[join("mkl_fft", "src"), np.get_include()] + include_dirs=[join("mkl_fft", "src"), numpy.get_include()] + mkl_include_dirs, libraries=mkl_libraries, library_dirs=mkl_library_dirs,