Take the 2-minute tour ×
Computational Science Stack Exchange is a question and answer site for scientists using computers to solve scientific problems. It's 100% free, no registration required.

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.

share|improve this question
1  
Now that I've asked this I wonder if it may be more appropriate for Signal Processing... but I'll let the community figure it out. –  David Z Jan 6 '12 at 23:21
add comment

4 Answers

up vote 22 down vote accepted

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.

share|improve this answer
add comment

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.

share|improve this answer
add comment

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.

share|improve this answer
1  
As far as I know, PNFFT is the fastest library for parallel 3d non-uniform FFT's. –  Jack Poulson Jan 17 '12 at 22:59
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.