How do I take the FFT of unevenly spaced data?

Question

The Fast Fourier Transform algorithm computes a Fourier decomposition under the assumption that its input points are equally spaced in the time domain, $t_k = kT$. What if they're not? Is there another algorithm I could use, or some way I could modify the FFT, to account for what is effectively a variable sampling rate?

If the solution depends on how the samples are distributed, there are two particular situations I'm most interested in:

• Constant sampling rate with jitter: $t_k = kT + \delta t_k$ where $\delta t_k$ is a randomly distributed variable. Suppose it's safe to say $|\delta t_k| < T/2$.
• Dropped samples from an otherwise constant sampling rate: $t_k = n_k T$ where $n_k \in\mathbb{Z}\ge k$

Motivation: first of all, this was one of the higher voted questions on the proposal for this site. But in addition, a while ago I got involved in a discussion about FFT usage (prompted by a question on Stack Overflow) in which some input data with unevenly sampled points came up. It turned out that the timestamps on the data were wrong, but it got me thinking about how one could tackle this problem.

Answer 1

There is a wide variety of techniques for non-uniform FFT, and the most efficient ones are all meant for exactly your case: quasi-uniform samples. The basic idea is to smear the unevenly sampled sources onto a slightly finer ("oversampled") uniform grid though local convolutions against Gaussians. A standard FFT can then be run on the oversampled uniform grid, and then the convolution against the Gaussians can be undone. Good implementations are something like $C^d$ times more expensive than a standard FFT in $d$ dimensions, where $C$ is something close to 4 or 5.

I recommend reading Accelerating the Nonuniform Fast Fourier Transform by Greengard and Lee.

There also exist fast, i.e., $O(N^d \log N)$ or faster, techniques when the sources and/or evaluation points are sparse, and there are also generalizations to more general integral operators, e.g., Fourier Integral Operators. If you are interested in these techniques, I recommend Sparse Fourier transform via butterfly algorithm and A fast butterfly algorithm for the computation of Fourier Integral Operators. The price paid in these techniques versus standard FFT's is a much higher coefficient. Disclaimer: My advisor wrote/cowrote those two papers, and I have spent a decent amount of time parallelizing those techniques.

An important point is that all of the above techniques are approximations that can be made arbitrarily accurate at the expense of longer runtimes, whereas the standard FFT algorithm is exact.

Answer 2

In signal processing, aliasing is avoided by sending a signal through a low pass filter before sampling. Jack Poulson already explained one technique for non-uniform FFT using truncated Gaussians as low pass filters. One inconvenient feature of truncated Gaussians is that even after you have decided on the grid spacing for the FFT (=the sampling rate in signal processing), you still have two free parameters: The width of the Gaussian and the truncation radius.

I therefore prefer the "hat" function with a width of two grid cells as low pass filter. This has the effect that the zeroth Fourier order is exact, and that the lower Fourier orders will converge quadratically. The Fourier transform of the "hat" function is easy to compute (it is the square of the sinc function), which simplifies undoing the convolution after the FFT. Note that the "hat" function is the convolution of the characteristic function of the (centered) unit cell with itself. Any desired convergence rate can be achieved, by convoluting the unit cell more than once with itself, and using the resulting function instead of the "hat" function.

Answer 3

If you are interested in software I can recommend the NFFT library (in C with an interface to MATLAB) which can be found here. Note that there is also a PFFT library for parallel FFT computation and even a PNFFT library for parallel non-equispaces FFTs by the same developers.

Answer 4

Addition to the accepted answer. Here is a link to an open source implementation of Greengard's and Lee's method: https://finufft.readthedocs.io/en/latest/ It has wrappers to C, fortran, MATLAB, octave, and python. I believe the FINUFFT is written in C++.

It is maintained and used at NYU Courant institute, SFU, Flatiron institute (obviously), University of Texas Austin and Florida state university. At least these are the ones I know of.

I am myself using an older version, because I am lazy. See: https://cims.nyu.edu/cmcl/nufft/nufft.html

Answer 5

Of interest may be the Date Compensated Discrete Fourier Transform:

Ferraz-Mello, S., 1981, Estimation of periods from unequally spaced observations, The Astronomical Journal, 302:757-763.