pyfftw
- The core¶
The core of pyfftw
consists of the FFTW
class,
wisdom functions and a couple of
utility functions for dealing with aligned
arrays.
This module represents the full interface to the underlying FFTW
library. However, users may find it easier to
use the helper routines provided in pyfftw.builders
. Default values
used by the helper routines can be controlled as via
configuration variables.
FFTW Class¶
- class pyfftw.FFTW(input_array, output_array, axes=(- 1,), direction='FFTW_FORWARD', flags=('FFTW_MEASURE',), threads=1, planning_timelimit=None)¶
FFTW is a class for computing a variety of discrete Fourier transforms of multidimensional, strided arrays using the FFTW library. The interface is designed to be somewhat pythonic, with the correct transform being inferred from the dtypes of the passed arrays.
The exact scheme may be either directly specified with the
direction
parameter or inferred from the dtypes and relative shapes of the input arrays. Information on which shapes and dtypes imply which transformations is available in the FFTW schemes. If a match is found, the plan corresponding to that scheme is created, operating on the arrays that are passed in. If no scheme can be created then aValueError
is raised.The actual transformation is performed by calling the
execute()
method.The arrays can be updated by calling the
update_arrays()
method.The created instance of the class is itself callable, and can perform the execution of the FFT, both with or without array updates, returning the result of the FFT. Unlike calling the
execute()
method, calling the class instance will also optionally normalise the output as necessary. Additionally, calling with an input array update will also coerce that array to be the correct dtype.See the documentation on the
__call__()
method for more information.Arguments:
input_array
andoutput_array
should be numpy arrays. The contents of these arrays will be destroyed by the planning process during initialisation. Information on supported dtypes for the arrays is given below.axes
describes along which axes the DFT should be taken. This should be a valid list of axes. Repeated axes are only transformed once. Invalid axes will raise anIndexError
exception. This argument is equivalent to the same argument innumpy.fft.fftn()
, except for the fact that the behaviour of repeated axes is different (numpy.fft
will happily take the fft of the same axis if it is repeated in theaxes
argument). Rudimentary testing has suggested this is down to the underlying FFTW library and so unlikely to be fixed in these wrappers.The
direction
parameter describes what sort of transformation the object should compute. This parameter is poorly named for historical reasons: older versions of pyFFTW only supported forward and backward transformations, for which this name made sense. Since then pyFFTW has been expanded to support real to real transforms as well and the name is not quite as descriptive.direction
should either be a string, or, in the case of multiple real transforms, a list of strings. The two values corresponding to the DFT are'FFTW_FORWARD'
, which is the forward discrete Fourier transform, and'FFTW_BACKWARD'
, which is the backward discrete Fourier transform.
Note that, for the two above options, only the Complex schemes allow a free choice for
direction
. The direction must agree with the the table below if a Real scheme is used, otherwise aValueError
is raised.Alternatively, if you are interested in one of the real to real transforms, then pyFFTW supports four different discrete cosine transforms:
'FFTW_REDFT00'
,'FFTW_REDFT01'
,'FFTW_REDFT10'
, and'FFTW_REDFT01'
,
and four discrete sine transforms:
'FFTW_RODFT00'
,'FFTW_RODFT01'
,'FFTW_RODFT10'
, and'FFTW_RODFT01'
.
pyFFTW uses the same naming convention for these flags as FFTW: the
'REDFT'
part of the name is an acronym for ‘real even discrete Fourier transform, and, similarly,'RODFT'
stands for ‘real odd discrete Fourier transform’. The trailing'0'
is notation for even data (in terms of symmetry) and the trailing'1'
is for odd data.Unlike the plain discrete Fourier transform, one may specify a different real to real transformation over each axis: for example,
will create a transformation across the first and last axes with a discrete cosine transform over the first and a discrete sine transform over the last.
Unfortunately, since this class is ultimately just a wrapper for various transforms implemented in FFTW, one cannot combine real transformations with real to complex transformations in a single object.
flags
is a list of strings and is a subset of the flags that FFTW allows for the planners:'FFTW_ESTIMATE'
,'FFTW_MEASURE'
,'FFTW_PATIENT'
and'FFTW_EXHAUSTIVE'
are supported. These describe the increasing amount of effort spent during the planning stage to create the fastest possible transform. Usually'FFTW_MEASURE'
is a good compromise. If no flag is passed, the default'FFTW_MEASURE'
is used.'FFTW_UNALIGNED'
is supported. This tells FFTW not to assume anything about the alignment of the data and disabling any SIMD capability (see below).'FFTW_DESTROY_INPUT'
is supported. This tells FFTW that the input array can be destroyed during the transform, sometimes allowing a faster algorithm to be used. The default behaviour is, if possible, to preserve the input. In the case of the 1D Backwards Real transform, this may result in a performance hit. In the case of a backwards real transform for greater than one dimension, it is not possible to preserve the input, making this flag implicit in that case. A little more on this is given below.'FFTW_WISDOM_ONLY'
is supported. This tells FFTW to raise an error if no plan for this transform and data type is already in the wisdom. It thus provides a method to determine whether planning would require additional effort or the cached wisdom can be used. This flag should be combined with the various planning-effort flags ('FFTW_ESTIMATE'
,'FFTW_MEASURE'
, etc.); if so, then an error will be raised if wisdom derived from that level of planning effort (or higher) is not present. If no planning-effort flag is used, the default of'FFTW_ESTIMATE'
is assumed. Note that wisdom is specific to all the parameters, including the data alignment. That is, if wisdom was generated with input/output arrays with one specific alignment, using'FFTW_WISDOM_ONLY'
to create a plan for arrays with any different alignment will cause the'FFTW_WISDOM_ONLY'
planning to fail. Thus it is important to specifically control the data alignment to make the best use of'FFTW_WISDOM_ONLY'
.
The FFTW planner flags documentation has more information about the various flags and their impact. Note that only the flags documented here are supported.
threads
tells the wrapper how many threads to use when invoking FFTW, with a default of 1. If the number of threads is greater than 1, then the GIL is released by necessity.planning_timelimit
is a floating point number that indicates to the underlying FFTW planner the maximum number of seconds it should spend planning the FFT. This is a rough estimate and corresponds to calling offftw_set_timelimit()
(or an equivalent dependent on type) in the underlying FFTW library. IfNone
is set, the planner will run indefinitely until all the planning modes allowed by the flags have been tried. See the FFTW planner flags page for more information on this.
Schemes
The currently supported full (so not discrete sine or discrete cosine) DFT schemes are as follows:
Type
input_array.dtype
output_array.dtype
Direction
Complex
complex64
complex64
Both
Complex
complex128
complex128
Both
Complex
clongdouble
clongdouble
Both
Real
float32
complex64
Forwards
Real
float64
complex128
Forwards
Real
longdouble
clongdouble
Forwards
Real^{1}
complex64
float32
Backwards
Real^{1}
complex128
float64
Backwards
Real^{1}
clongdouble
longdouble
Backwards
^{1} Note that the Backwards Real transform for the case in which the dimensionality of the transform is greater than 1 will destroy the input array. This is inherent to FFTW and the only general work-around for this is to copy the array prior to performing the transform. In the case where the dimensionality of the transform is 1, the default is to preserve the input array. This is different from the default in the underlying library, and some speed gain may be achieved by allowing the input array to be destroyed by passing the
'FFTW_DESTROY_INPUT'
flag.The discrete sine and discrete cosine transforms are supported for all three real types.
clongdouble
typically maps directly tocomplex256
orcomplex192
, andlongdouble
tofloat128
orfloat96
, dependent on platform.The relative shapes of the arrays should be as follows:
For a Complex transform,
output_array.shape == input_array.shape
For a Real transform in the Forwards direction, both the following should be true:
output_array.shape[axes][-1] == input_array.shape[axes][-1]//2 + 1
All the other axes should be equal in length.
For a Real transform in the Backwards direction, both the following should be true:
input_array.shape[axes][-1] == output_array.shape[axes][-1]//2 + 1
All the other axes should be equal in length.
In the above expressions for the Real transform, the
axes
arguments denotes the unique set of axes on which we are taking the FFT, in the order passed. It is the last of these axes that is subject to the special case shown.The shapes for the real transforms corresponds to those stipulated by the FFTW library. Further information can be found in the FFTW documentation on the real DFT.
The actual arrangement in memory is arbitrary and the scheme can be planned for any set of strides on either the input or the output. The user should not have to worry about this and any valid numpy array should work just fine.
What is calculated is exactly what FFTW calculates. Notably, this is an unnormalized transform so should be scaled as necessary (fft followed by ifft will scale the input by N, the product of the dimensions along which the DFT is taken). For further information, see the FFTW documentation.
The FFTW library benefits greatly from the beginning of each DFT axes being aligned on the correct byte boundary, enabling SIMD instructions. By default, if the data begins on such a boundary, then FFTW will be allowed to try and enable SIMD instructions. This means that all future changes to the data arrays will be checked for similar alignment. SIMD instructions can be explicitly disabled by setting the FFTW_UNALIGNED flags, to allow for updates with unaligned data.
byte_align()
andempty_aligned()
are two methods included with this module for producing aligned arrays.The optimum alignment for the running platform is provided by
pyfftw.simd_alignment
, though a different alignment may still result in some performance improvement. For example, if the processor supports AVX (requiring 32-byte alignment) as well as SSE (requiring 16-byte alignment), then if the array is 16-byte aligned, SSE will still be used.It’s worth noting that just being aligned may not be sufficient to create the fastest possible transform. For example, if the array is not contiguous (i.e. certain axes are displaced in memory), it may be faster to plan a transform for a contiguous array, and then rely on the array being copied in before the transform (which
pyfftw.FFTW
will handle for you when accessed through__call__()
).- N¶
The product of the lengths of the DFT over all DFT axes. 1/N is the normalisation constant. For any input array A, and for any set of axes, 1/N * ifft(fft(A)) = A
- simd_aligned¶
Return whether or not this FFTW object requires simd aligned input and output data.
- input_alignment¶
Returns the byte alignment of the input arrays for which the
FFTW
object was created.Input array updates with arrays that are not aligned on this byte boundary will result in a ValueError being raised, or a copy being made if the
__call__()
interface is used.
- output_alignment¶
Returns the byte alignment of the output arrays for which the
FFTW
object was created.Output array updates with arrays that are not aligned on this byte boundary will result in a ValueError being raised.
- flags¶
Return which flags were used to construct the FFTW object.
This includes flags that were added during initialisation.
- input_array¶
Return the input array that is associated with the FFTW instance.
- output_array¶
Return the output array that is associated with the FFTW instance.
- input_shape¶
Return the shape of the input array for which the FFT is planned.
- output_shape¶
Return the shape of the output array for which the FFT is planned.
- input_strides¶
Return the strides of the input array for which the FFT is planned.
- output_strides¶
Return the strides of the output array for which the FFT is planned.
- input_dtype¶
Return the dtype of the input array for which the FFT is planned.
- output_dtype¶
Return the shape of the output array for which the FFT is planned.
- direction¶
Return the planned FFT direction. Either ‘FFTW_FORWARD’, ‘FFTW_BACKWARD’, or a list of real transform codes of the form [‘FFTW_R*DFT**’].
- axes¶
Return the axes for the planned FFT in canonical form. That is, as a tuple of positive integers. The order in which they were passed is maintained.
- ortho¶
If
ortho=True
both the forward and inverse transforms are scaled by 1/sqrt(N).
- normalise_idft¶
If
normalise_idft=True
, the inverse transform is scaled by 1/N.
- __call__()¶
- __call__(input_array=None, output_array=None, normalise_idft=True,
ortho=False)
Calling the class instance (optionally) updates the arrays, then calls
execute()
, before optionally normalising the output and returning the output array.It has some built-in helpers to make life simpler for the calling functions (as distinct from manually updating the arrays and calling
execute()
).If
normalise_idft
isTrue
(the default), then the output from an inverse DFT (i.e. when the direction flag is'FFTW_BACKWARD'
) is scaled by 1/N, where N is the product of the lengths of input array on which the FFT is taken. If the direction is'FFTW_FORWARD'
, this flag makes no difference to the output array.If
ortho
isTrue
, then the output of both forward and inverse DFT operations is scaled by 1/sqrt(N), where N is the product of the lengths of input array on which the FFT is taken. This ensures that the DFT is a unitary operation, meaning that it satisfies Parseval’s theorem (the sum of the squared values of the transform output is equal to the sum of the squared values of the input). In other words, the energy of the signal is preserved.If either
normalise_idft
orortho
areTrue
, then ifft(fft(A)) = A.When
input_array
is something other than None, then the passed in array is coerced to be the same dtype as the input array used when the class was instantiated, the byte-alignment of the passed in array is made consistent with the expected byte-alignment and the striding is made consistent with the expected striding. All this may, but not necessarily, require a copy to be made.As noted in the scheme table, if the FFTW instance describes a backwards real transform of more than one dimension, the contents of the input array will be destroyed. It is up to the calling function to make a copy if it is necessary to maintain the input array.
output_array
is always used as-is if possible. If the dtype, the alignment or the striding is incorrect for the FFTW object, then aValueError
is raised.The coerced input array and the output array (as appropriate) are then passed as arguments to
update_arrays()
, after whichexecute()
is called, and then normalisation is applied to the output array if that is desired.Note that it is possible to pass some data structure that can be converted to an array, such as a list, so long as it fits the data requirements of the class instance, such as array shape.
Other than the dtype and the alignment of the passed in arrays, the rest of the requirements on the arrays mandated by
update_arrays()
are enforced.A
None
argument to either keyword means that that array is not updated.The result of the FFT is returned. This is the same array that is used internally and will be overwritten again on subsequent calls. If you need the data to persist longer than a subsequent call, you should copy the returned array.
- update_arrays(new_input_array, new_output_array)¶
Update the arrays upon which the DFT is taken.
The new arrays should be of the same dtypes as the originals, the same shapes as the originals and should have the same strides between axes. If the original data was aligned so as to allow SIMD instructions (e.g. by being aligned on a 16-byte boundary), then the new array must also be aligned so as to allow SIMD instructions (assuming, of course, that the
FFTW_UNALIGNED
flag was not enabled).The byte alignment requirement extends to requiring natural alignment in the non-SIMD cases as well, but this is much less stringent as it simply means avoiding arrays shifted by, say, a single byte (which invariably takes some effort to achieve!).
If all these conditions are not met, a
ValueError
will be raised and the data will not be updated (though the object will still be in a sane state).
- execute()¶
Execute the planned operation, taking the correct kind of FFT of the input array (i.e.
FFTW.input_array
), and putting the result in the output array (i.e.FFTW.output_array
).
- get_input_array()¶
Return the input array that is associated with the FFTW instance.
Deprecated since 0.10. Consider using the
FFTW.input_array
property instead.
- get_output_array()¶
Return the output array that is associated with the FFTW instance.
Deprecated since 0.10. Consider using the
FFTW.output_array
property instead.
Wisdom Functions¶
Functions for dealing with FFTW’s ability to export and restore plans, referred to as wisdom. For further information, refer to the FFTW wisdom documentation.
- pyfftw.export_wisdom()¶
Return the FFTW wisdom as a tuple of strings.
The first string in the tuple is the string for the double precision wisdom, the second is for single precision, and the third for long double precision. If any of the precisions is not supported in the build, the string is empty.
The tuple that is returned from this function can be used as the argument to
import_wisdom()
.
- pyfftw.import_wisdom(wisdom)¶
Function that imports wisdom from the passed tuple of strings.
The first string in the tuple is the string for the double precision wisdom. The second string in the tuple is the string for the single precision wisdom. The third string in the tuple is the string for the long double precision wisdom.
The tuple that is returned from
export_wisdom()
can be used as the argument to this function.This function returns a tuple of boolean values indicating the success of loading each of the wisdom types (double, float and long double, in that order).
- pyfftw.forget_wisdom()¶
Forget all the accumulated wisdom.
Utility Functions¶
- pyfftw.simd_alignment¶
An integer giving the optimum SIMD alignment in bytes, found by inspecting the CPU (e.g. if AVX is supported, its value will be 32).
This can be used as
n
in the arguments forbyte_align()
,empty_aligned()
,zeros_aligned()
, andones_aligned()
to create optimally aligned arrays for the running platform.
- pyfftw.byte_align(array, n=None, dtype=None)¶
Function that takes a numpy array and checks it is aligned on an n-byte boundary, where
n
is an optional parameter. Ifn
is not provided then this function will inspect the CPU to determine alignment. If the array is aligned then it is returned without further ado. If it is not aligned then a new array is created and the data copied in, but aligned on the n-byte boundary.dtype
is an optional argument that forces the resultant array to be of that dtype.
- pyfftw.empty_aligned(shape, dtype='float64', order='C', n=None)¶
Function that returns an empty numpy array that is n-byte aligned, where
n
is determined by inspecting the CPU if it is not provided.The alignment is given by the final optional argument,
n
. Ifn
is not provided then this function will inspect the CPU to determine alignment. The rest of the arguments are as pernumpy.empty()
.
- pyfftw.zeros_aligned(shape, dtype='float64', order='C', n=None)¶
Function that returns a numpy array of zeros that is n-byte aligned, where
n
is determined by inspecting the CPU if it is not provided.The alignment is given by the final optional argument,
n
. Ifn
is not provided then this function will inspect the CPU to determine alignment. The rest of the arguments are as pernumpy.zeros()
.
- pyfftw.ones_aligned(shape, dtype='float64', order='C', n=None)¶
Function that returns a numpy array of ones that is n-byte aligned, where
n
is determined by inspecting the CPU if it is not provided.The alignment is given by the final optional argument,
n
. Ifn
is not provided then this function will inspect the CPU to determine alignment. The rest of the arguments are as pernumpy.ones()
.
- pyfftw.is_byte_aligned()¶
is_n_byte_aligned(array, n=None)
Function that takes a numpy array and checks it is aligned on an n-byte boundary, where
n
is an optional parameter, returningTrue
if it is, andFalse
if it is not. Ifn
is not provided then this function will inspect the CPU to determine alignment.
- pyfftw.n_byte_align(array, n, dtype=None)¶
This function is deprecated:
byte_align
should be used instead.Function that takes a numpy array and checks it is aligned on an n-byte boundary, where
n
is an optional parameter. Ifn
is not provided then this function will inspect the CPU to determine alignment. If the array is aligned then it is returned without further ado. If it is not aligned then a new array is created and the data copied in, but aligned on the n-byte boundary.dtype
is an optional argument that forces the resultant array to be of that dtype.
- pyfftw.n_byte_align_empty(shape, n, dtype='float64', order='C')¶
This function is deprecated:
empty_aligned
should be used instead.Function that returns an empty numpy array that is n-byte aligned.
The alignment is given by the first optional argument,
n
. Ifn
is not provided then this function will inspect the CPU to determine alignment. The rest of the arguments are as pernumpy.empty()
.
- pyfftw.is_n_byte_aligned(array, n)¶
This function is deprecated:
is_byte_aligned
should be used instead.Function that takes a numpy array and checks it is aligned on an n-byte boundary, where
n
is a passed parameter, returningTrue
if it is, andFalse
if it is not.
- pyfftw.next_fast_len(target)¶
Find the next fast transform length for FFTW.
FFTW has efficient functions for transforms of length 2**a * 3**b * 5**c * 7**d * 11**e * 13**f, where e + f is either 0 or 1.
- Parameters
target (int) – Length to start searching from. Must be a positive integer.
- Returns
out – The first fast length greater than or equal to target.
- Return type
Examples
On a particular machine, an FFT of prime length takes 2.1 ms:
>>> from pyfftw.interfaces import scipy_fftpack >>> min_len = 10007 # prime length is worst case for speed >>> a = numpy.random.randn(min_len) >>> b = scipy_fftpack.fft(a)
Zero-padding to the next fast length reduces computation time to 406 us, a speedup of ~5 times:
>>> next_fast_len(min_len) 10080 >>> b = scipy_fftpack.fft(a, 10080)
Rounding up to the next power of 2 is not optimal, taking 598 us to compute, 1.5 times as long as the size selected by next_fast_len.
>>> b = fftpack.fft(a, 16384)
Similar speedups will occur for pre-planned FFTs as generated via pyfftw.builders.
FFTW Configuration¶
- pyfftw.config.NUM_THREADS¶
This variable controls the default number of threads used by the functions in
pyfftw.builders
andpyfftw.interfaces
.The default value is read from the environment variable
PYFFTW_NUM_THREADS
. If this variable is undefined and the user’s underlying FFTW library was built using OpenMP threading, the number of threads will be read from the environment variableOMP_NUM_THREADS
instead. If neither environment variable is defined, the default value is 1.If the specified value is
<= 0
, the library will usemultiprocessing.cpu_count()
to determine the number of threads.The user can modify the value at run time by assigning to this variable.
- pyfftw.config.PLANNER_EFFORT¶
This variable controls the default planning effort used by the functions in
pyfftw.builders
andpyfftw.interfaces
.The default value of is determined by reading from the environment variable
PYFFTW_PLANNER_EFFORT
. If this environment variable is undefined, it defaults to'FFTW_ESTIMATE'
.The user can modify the value at run time by assigning to this variable.