Discussion : Generic interfaces for the Fourier transforms

[loiseaujc]

Hej everyone,

Given how far we've gone with the linear algebra module, I think it is time to think a bit about what could be the next "big" addition to stdlib. Following some old discussions (see #1, #20, #21, #66, #409 and #540), introducing fast Fourier transforms capabilities seem like the natural next step. Having access to the FFT (and its variants) would allow for numerous additional functionalities, including (but not limited to):

• signal processing: spectral analysis, convolution, auto-correlation, cross-correlation, etc.
• partial differential equations: spectral and pseudo-spectral solvers, fast elliptic (e.g. Poisson) solver, etc.
• linear algebra: circulant matrices, fast matrix-vector product for Toeplitz and Hankel matrices, etc.

There are many FFT libraries available out there, with different licenses and API. From reading the previous discussions, it seems to me that the consensus back then was something along the following lines:

• Have fftpack as a stand-alone package.
• Implement some high-level interfaces (e.g. y = fft(x), p = plan_fft(x)) in stdlib which would allow to switch between different backends (to alleviate licencing issues with FFTW for instance) with fftpack being the default one.
• Have wrappers for the most common backends, e.g. fftpack, FFTW, MKL, and possibly a couple of others.

Note that if fftpack is the default backend, we would need to keep the modernization effort started by @jacobwilliams in order to cover all the different precisions by stdlib.

@certik, @milancurcic, @ivan-pi : is my understanding of the discussions you had a couple of years ago correct ? Do you still stand by this point of view or has it evolved? @perazz, @jvdp1, @jalvesz : Given your experience with developing the linear algebra module, do you have insights that could be useful for the fft endeavor?

[jalvesz]

I'm not familiar with FFTW or the MKL APIs for FFT, but my understanding is that MKL actually proposes a drop-in replacement with the same API signature. So I guess that a similar approach to what was done for BLAS/LAPACK can be done here: integrate fftpack within stdlib as the default backend and have a C-preprocessing macro with something like STDLIB_EXTERNAL_FFT to link against precompiled-optimized versions (FFTW or MKL).

[loiseaujc]

Some heads up: I've continued the modernization of fftpack, most notably mimicking the CMake build (and Github Actions) of stdlib and replacing a whole bunch of nested loops with do concurrent. On the small benchmark-ish examples included in the repo, it leads to a 1.2x to 1.4x speed-up compared to the legacy version. There is still some work I'd like to do (e.g. module/submodule structure and fypp-ification for supporting all intrinsic data types for instance), but I think we could start discussing a bit more seriously about the high-level API before figuring out all of the nitty-gritty details of supporting multiple backends.

Here are some random thoughts in no particular order.

High-level and low-level interfaces

Similarly to what has been done for the linalg module, we should expose two interfaces for the most important functionalities. Looking at what scipy is doing, the high-level interface for the 1D fast Fourier transform could look like

function fft(x, n, norm, plan, err) result(y)
    complex(dp), intent(in) :: x(:)
    !! Input array.
    integer(ilp), optional, intent(in) :: n
    !! Length of the output array. If n < size(x), the input is cropped.
    !! If n > size(x), it is zero-padded.
    character(:), optional, intent(in) :: norm
    !! How the FFT should be normalized (backward, forward, ortho)
    type(fft_plan_type), optional, intent(in) :: plan
    !! Pre-computed plan.
    type(fft_state_type), optional, intent(out) :: err
    complex(dp), allocatable :: y(:)
end function fft

The low-level version wouldn't be too different. Probably something along the lines of

subroutine compute_fft(x, n, norm, plan, err)
    complex(dp), intent(inout) :: x(:)
    integer(ilp), optional, intent(in) :: n
    character(:), optional, intent(in) :: norm
    type(fft_plan_type), optional, intent(in) :: plan
    type(fft_state_type), optional, intent(out) :: err
end subroutine

The fft_state_type would be a new derived-type mimicking the linalg_state_type. The definition of the fft_plan_type would use conditional compilation to ensure it contains anything required by the different FFT backends eventually supported (e.g pre-allocated arrays, informations about the data to be transformed, etc).

Multi-dimensional Fourier transforms

I guess there are no need for discussion about including a 2D FFT (fft2 and ifft2) and a generic multi-dimensional one (fftn and ifftn). Note however that, currently, fftpack does not provide such functionalities.

Fourier-like transforms

Currently, fftpack provides implementations for:

• Fast Fourier transform
• Discrete Sine and Cosine transforms
• Quarter-wave Sine and Cosine transforms (I think these are essentially different types of dct)

If I recall correctly, scipy also has functions for computing the Hilbert transform (which I do use regularly) and the Hankel one (which I've never came across). Another low-hanging fruit that could be useful for people using spectral methods for PDE is the Chebyshev transform (essentially a dct with modified input).

Additional functions

These would be the standard fftfreq, fftshift, etc. Some functions for computing the next fast size of input data (e.g. next_fast_len in scipy). Similarly, a convolve function would be nice.

---

From a practical point of view, I guess that many of these features could be crafted directly in fftpack since they would be useful there independently of using stdlib or not. Once we have something that sort of fits everything we want, we can have stdlib directly re-export these interfaces from fftpack and then figure out how to handle different backends.

Ping: @zoziha since they're the one reviewing my PR in fftpack :]

[jvdp1]

Nice work, @loiseaujc . I am not familiar with fft neither a user of fftpack. However, the proposed high-level and low-level interfaces make sense to me, and aligns well with the linalg module and @perazz work on BLAS and LAPACK. Handling different backends can be done later. I think it makes also sense to implement the new features in the fftpack project directly, especically as it is hosted in fortran-lang. Therefore, I wonder if the code should be copied in stdlib or if the user should rely on cmake or on fpm to load the code from the fftpack project as a dependency.

[loiseaujc]

I wonder if the code should be copied in stdlib or if the user should rely on cmake or on fpm to load the code from the fftpack project as a dependency.

I think the initial idea of @certik, @milancurcic and others was to have fftpack as a standalone package, expose high-level interfaces inside stdlib and use conditional compilation to switch between the standard fftpack backend or other optimized ones (e.g. mkl, fftw, etc). For the sake of development, I believe that it does make sense to keep it as standalone. Even more so given we can partially (and relatively easily) mimic what has been done by @perazz for blas/lapack.

Side note - The current implementation of in the fortran-lang repo is based on fftpack 4.0 (dfftpack 1.0 on netlib). On the other hand, fftpack 5.1 has been released some years ago (see here and here) for instance. One question thus is: should we start crafting the high-level interfaces right away and update to fftpack 5.1 later on or start with updating fftpack and craft the interfaces afterward?

Craft interfaces now and update fttpack later on

• Pros - Fourier analysis capabilities would be available earlier in stdlib.
• Cons - Interfaces may have to be revamped after updating fftpack due to possibly changed definitions.

Update fftpack now and craft interfaces later on

• Pros - The high-level interfaces are directly based on the latest version of fftpack and would thus be less likely to change.
• Cons - Fourier analysis capabilities in stdlib would be delayed by some time.

I have conflicting interests here. As a heavy stdlib user, I wouldn't mind having access to fft as soon as possible from stdlib directly. As the main person working on modernizing the fortran-lang/fftpack implementation at the moment, updating to the latest dev also makes a lot of sense. We would have to check the licence though. fortran-lang/fftpack is released under the MIT Licence while fftpack 5 is under the GPL one I believe.