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:
- 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