The following packages extend the functionality provided by AbstractFFTs:
- FFTW.jl: Bindings for the FFTW library. This also used to be part of Base Julia.
- FastTransforms.jl: Pure-Julia implementation of FFT, with support for arbitrary AbstractFloat types.
To define a new FFT implementation in your own module, you should
-
Define a new subtype (e.g.
MyPlan) ofAbstractFFTs.Plan{T}for FFTs and related transforms on arrays ofT. This must have apinv::Planfield, initially undefined when aMyPlanis created, that is used for caching the inverse plan. -
Define a new method
AbstractFFTs.plan_fft(x, region; kws...)that returns aMyPlanfor at least some types ofxand some set of dimensionsregion. Theregion(or a copy thereof) should be accessible viafftdims(p::MyPlan)(which defaults top.region), and the input sizesize(x)should be accessible viasize(p::MyPlan). -
Define a method of
LinearAlgebra.mul!(y, p::MyPlan, x)that computes the transformpofxand stores the result iny. -
Define a method of
*(p::MyPlan, x), which can simply call yourmul!method. This is not defined generically in this package due to subtleties that arise for in-place and real-input FFTs. -
If the inverse transform is implemented, you should also define
plan_inv(p::MyPlan), which should construct the inverse plan top, andplan_bfft(x, region; kws...)for an unnormalized inverse ("backwards") transform ofx. Implementations only need to provide the unnormalized backwards FFT, similar to FFTW, and we do the scaling generically to get the inverse FFT. -
You can also define similar methods of
plan_rfftandplan_brfftfor real-input FFTs. -
To support adjoints in a new plan, define the trait
AbstractFFTs.AdjointStyle.AbstractFFTsimplements the following adjoint styles:AbstractFFTs.FFTAdjointStyle,AbstractFFTs.RFFTAdjointStyle,AbstractFFTs.IRFFTAdjointStyle, andAbstractFFTs.UnitaryAdjointStyle. To define a new adjoint style, define the methodsAbstractFFTs.adjoint_mulandAbstractFFTs.output_size.
The normalization convention for your FFT should be that it computes y_k = \sum_j x_j \exp(-2\pi i j k/n) for a transform of
length n, and the "backwards" (unnormalized inverse) transform computes the same thing but with \exp(+2\pi i jk/n).