Performing FFTs in PyCBC
Introduction
Many applications in gravitational-wave analysis rely on Fast Fourier Transforms. These are often the dominant computational cost in analyses. PyCBC needs to balance the requirement that analyses be efficient with ease of use for end users. To meet this requirement PyCBC provides two different APIs to do FFTs:
A function based API, which is easy to use, but not optimized.
A class-based API, which is a little more involved to use, but allows the use of optimized FFT routines.
These APIs offer access to a number of FFT backends. PyCBC knows how to do FFTs using the FFTW, MKL and numpy backends, and will enable these if they are present on your system. By default FFTW will be used, then MKL if FFTW is not and numpy will be used only if the other 2 are not present. However, you can override this and choose a specific backend if multiple are available.
When running on GPUs, PyCBC knows how to do CUDA FFTs through the same interface.
Using the function based API
The PyCBC function based API offers a simple way to Fourier transform an input array into an output array. This is done by:
>>> from pycbc import fft
>>> fft.fft(input_array, output_array)
Or for an inverse FFT:
>>> from pycbc import fft
>>> fft.ifft(input_array, output_array)
To do this requires having the output_array, which would be the Fourier transform of the input, already existing as an array of 0s. This output array must be the correct length, the correct type (complex or real) and the correct precision (double or float precision). It’s also worth noting that fast Fourier transforms are more efficient if their lengths are 2**N where N is some integer. This becomes a little complicated if doing real->complex or complex->real transforms; in those cases the longer array should have a length of 2**N to be efficient.
Here’s a few examples:
>>> import numpy as np
>>> from pycbc import types
>>> from pycbc import fft
>>> inarr = types.Array(np.ones([64], dtype=np.complex64))
>>> outarr = types.Array(np.zeros([64], dtype=np.complex64))
>>> fft.fft(inarr, outarr)
>>> print(outarr)
or (note here the length of the complex array is the length of the real array divided by 2 and then + 1):
>>> import numpy as np
>>> from pycbc import types
>>> from pycbc import fft
>>> inarr = types.Array(np.ones([64], dtype=np.float32))
>>> outarr = types.Array(np.zeros([33], dtype=np.complex64))
>>> fft.fft(inarr, outarr)
>>> print(outarr)
or (this one is an inverse FFT):
>>> import numpy as np
>>> from pycbc import types
>>> from pycbc import fft
>>> inarr = types.Array(np.ones([33], dtype=np.complex64))
>>> outarr = types.Array(np.zeros([64], dtype=np.float32))
>>> fft.ifft(inarr, outarr)
>>> print(outarr)
This will work the pycbc Timeseries and Frequencyseries as well. Except you must FFT a TimeSeries to a FrequencySeries or IFFT a FrequencySeries to a TimeSeries. In this case the time and frequency spacing must also be consistent. For this reason we provide convenience functions that use the function API, but figure out these details, and create the output array, for you. As an example:
>>> import numpy as np
>>> from pycbc import types
>>> inarr = types.TimeSeries(np.ones([64], dtype=np.float64), delta_t=1./64.)
>>> outarr = inarr.to_frequencyseries()
or:
>>> import numpy as np
>>> from pycbc import types
>>> inarr = types.FrequencySeries(np.ones([33], dtype=np.complex128), delta_f=1.)
>>> outarr = inarr.to_timeseries()
Using the class-based API
The PyCBC class-based API should be used if you care about performance. If you are performing FFTs many times, with inputs that are the same size each time, using this will offer significance perfommance improvement.
Here’s how to use this:
>>> from pycbc import fft
>>> fft_class = pycbc.fft.FFT(inarr, outarr)
>>> fft_class.execute()
>>> outarr *= inarr._delta_t # ONLY IF inarr is a TimeSeries
>>> outarr *= inarr._delta_f # ONLY IF inarr is a FrequencySeries
Or for an inverse FFT:
>>> from pycbc import fft
>>> ifft_class = pycbc.fft.IFFT(inarr, outarr)
>>> ifft_class.execute()
>>> outarr *= inarr._delta_t # ONLY IF inarr is a TimeSeries
>>> outarr *= inarr._delta_f # ONLY IF inarr is a FrequencySeries
The idea would be that the fft_class or ifft_class would only be created once and the execute command called many times. You would change the contents of inarr before each call and outarr will update when execute is run. After creating the FFT class do not reassign inarr, but instead set values. So:
>>> fft_class = pycbc.fft.FFT(inarr, outarr)
>>> inarr = types.TimeSeries(np.ones([64], dtype=np.float64), delta_t=1./64.)
>>> fft_class.execute()
would not work! Instead do:
>>> fft_class = pycbc.fft.FFT(inarr, outarr)
>>> inarr[:] = np.ones([64], dtype=np.float64)[:]
>>> fft_class.execute()
Choosing a specific backend
If you want to choose a specific backend, you can see what is available with:
>>> from pycbc.fft import backend_support
>>> backend_support.get_backend_names()
and then do:
>>> from pycbc.fft import backend_support
>>> backend_support.set_backend(['mkl'])
to set a specific backend. Running:
>>> from pycbc.fft import backend_support
>>> backend_support.get_backend()
will tell you what you are currently using. You can also use the MKL Scheme to default to using MKL FFTs, instead of FFTW.