mmfutils.performance.fft

FFTW wrappers for high-performance computing.

This module requires you to have installed the fftw libraries and pyfftw. Note that you must build the fftw with all precisions using something like:

PREFIX=/data/apps/fftw
VER=3.3.4
for opt in " " "--enable-sse2 --enable-single" \
           "--enable-long-double" "--enable-quad-precision"; do
  ./configure --prefix="${PREFIX}/${VER}"\
              --enable-threads\
              --enable-shared\
              $opt
  make -j8 install
done

Note: The FFTW library does not work with negative indices for axis. Indices should first be normalized by inds % len(shape).

For best performance, you should call the appropriate get_*fft*(). By default, this uses _PLANNER_EFFORT = “FFTW_MEASURE” with _THREADS. For optimal performance, you need to adjust _THREADS or use _PLANNER_EFFORT = “FFTW_PATIENT”, but this can be extremely slow.

The methods fft, ifft, fftn`, and ifftn are provided for convenience. They call an appropriate builder upon first invocation, check that the returned function is faster than the NumPy version, and then cache this. The returned function checks the global flag _COPY_OUTPUT ensures that the resulting array has flags[‘OWNDATA’], otherwise a copy is made. If the fft function will not be called before the array is copied, you might gain some performance improvement by setting this to False.

mmfutils.performance.fft.fft(a, n=None, axis=-1)[source]

High-performance replacement for numpy.fft.fft()

mmfutils.performance.fft.fftfreq(n, d=1.0, device=None)

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.

  • device (str, optional) –

    The device on which to place the created array. Default: None. For Array-API interoperability only, so must be "cpu" if passed.

    Added in version 2.0.0.

Returns:

f – Array of length n containing the sample frequencies.

Return type:

ndarray

Examples

>>> import numpy as np
>>> 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])
mmfutils.performance.fft.fftn(a, s=None, axes=None)[source]

High-performance replacement for numpy.fft.fftn()

mmfutils.performance.fft.get_fft(a, n=None, axis=-1, repeat=3, **kw)[source]

Return the fastest function to compute the fft.

Parameters:

repeat (int) – Uses timeit.repeat() to compare versions (if pyfftw is defined).

Other arguments are as for numpy.fft.fft().

mmfutils.performance.fft.get_fftn(a, s=None, axes=None, repeat=3, **kw)[source]

Return the fastest function to compute the fftn.

Parameters:

repeat (int) – Uses timeit.repeat() to compare versions (if pyfftw is defined).

Other arguments are as for numpy.fft.fftn().

mmfutils.performance.fft.get_ifft(a, n=None, axis=-1, repeat=3, **kw)[source]

Return the fastest function to compute the ifft.

Parameters:

repeat (int) – Uses timeit.repeat() to compare versions (if pyfftw is defined).

Other arguments are as for numpy.fft.ifft().

mmfutils.performance.fft.get_ifftn(a, s=None, axes=None, repeat=3, **kw)[source]

Return the fastest function to compute the ifftn.

Parameters:

repeat (int) – Uses timeit.repeat() to compare versions (if pyfftw is defined).

Other arguments are as for numpy.fft.ifftn().

mmfutils.performance.fft.ifft(a, n=None, axis=-1)[source]

High-performance replacement for numpy.fft.ifft()

mmfutils.performance.fft.ifftn(a, s=None, axes=None)[source]

High-performance replacement for numpy.fft.ifftn()

mmfutils.performance.fft.resample(f, N)[source]

Resample f to a new grid of size N.

This uses the FFT to resample the function f on a new grid with N points. Note: this assumes that the function f is periodic. Resampling non-periodic functions to finer lattices may introduce aliasing artifacts.

Parameters:
  • f (array) – The function to be resampled. May be n-dimensional

  • N (int or array) – The number of lattice points in the new array. If this is an integer, then all dimensions of the output array will have this length.

Examples

>>> def f(x, y):
...     "Function with only low frequencies"
...     return (np.sin(2*np.pi*x)-np.cos(4*np.pi*y))
>>> L = 1.0
>>> Nx, Ny = 16, 13   # Small grid
>>> NX, NY = 31, 24   # Large grid
>>> dx, dy = L/Nx, L/Ny
>>> dX, dY = L/NX, L/NY
>>> x = (np.arange(Nx)*dx - L/2)[:, None]
>>> y = (np.arange(Ny)*dy - L/2)[None, :]
>>> X = (np.arange(NX)*dX - L/2)[:, None]
>>> Y = (np.arange(NY)*dY - L/2)[None, :]
>>> f_XY = resample(f(x,y), (NX, NY))
>>> np.allclose(f_XY, f(X,Y))                      # To larger grid
True
>>> np.allclose(resample(f_XY, (Nx, Ny)), f(x,y))  # Back down
True