Disagreement w/ FFT in limit of uniform grid

[emd]

In your implementation walkthrough, you show that your nfft implementation agrees with your ndft implementation. However, you never show that either agree with the FFT in the limit of a uniform grid, and, having never considered the NFFT before, I figured this would be a natural place for me to start before I began using non-uniform grids. I am particularly interested in transforming measurements made on a non-uniform grid into Fourier space, so my discussion below will be focused on the adjoint transform.

Using the below MWE, I find that, in the limit of a uniform grid, the ndft_adjoint implementation disagrees with the FFT:

import numpy as np
import matplotlib.pyplot as plt

# Construct *uniform* timebase to benchmark `ndft_adjoint` against `fft`.
# (Note: `nfft` docs state that algorithms expect -0.5 <= t < 0.5.)
Npts = 128
dt = 1. / Npts
t = np.arange(-0.5, 0.5, dt)

# Construct signal. (As we will see, the problem is in the phase
# of the adjoint NDFT computation. To ensure that the phase
# is not computed from numerical noise, we should select a
# signal frequency such that
#
#   (signal frequency) * (observation window) != integer
#
# This ensures that there is significant spectral leakage
# to compute a sensible phase outside of the spectral peaks.
f0 = 5.1
x = np.sin(2 * np.pi * f0 * t)

# Compute FFT, shifting zero-frequency component of spectrum to center.
f_fft = np.fft.fftshift(np.fft.fftfreq(len(x), d=dt))
xhat_fft = np.fft.fftshift(np.fft.fft(x))

# Compute adjoint NDFT on *uniform* grid -- should be equal to FFT
xhat_ndft_adj = ndft_adjoint(t, x, len(f_fft))

# ... but `fft` basis is exp(-2j * pi * f * t), whereas the
# `ndft_adjoint` basis is exp(+2j * pi * f * t). Thus, benchmarking
# against the FFT requires the reverse indexing of `xhat_ndft_adj.
xhat_ndft_adj = xhat_ndft_adj[::-1]

# The frequencies corresponding to `ndft_adjoint` must also be
# appropriately reversed (and negated) for comparison to the FFT.
f_ndft_adj = -f_fft[::-1]

# Finally, plot magnitude and phase.
fig, axes = plt.subplots(2, 1, sharex=True)

# Magnitudes are in agreement...
axes[0].semilogy(f_fft, np.abs(xhat_fft), '-o', label='fft')
axes[0].semilogy(f_ndft_adj, np.abs(xhat_ndft_adj), label='ndft_adjoint')
axes[0].set_ylabel('magnitude [arbitrary]')
axes[0].legend()

# ... but phase is way off for default `nfft` NDFT adjoint algorithm
axes[1].plot(f_fft, np.unwrap(np.angle(xhat_fft)), '-o')
axes[1].plot(f_ndft_adj, np.unwrap(np.angle(xhat_ndft_adj)))
axes[1].set_xlabel('frequency')
axes[1].set_ylabel('phase [radians]')

plt.show()

Making a simple modification to your ndft_adjoint code, as detailed below, and re-running the above MWE produces agreement with the FFT.

def ndft_adjoint(x, f, N):
    """Compute the adjoint non-equispaced direct Fourier transform

    \hat{f}_k = \sum_{0 \le j < N} f_j \exp(2 \pi i k x_j)

    where k = range(-N/2, N/2)

    Parameters
    ----------
    x : array_like, shape=(M,)
        The locations of the data points.
    f : array_like, shape=(M,)
        The amplitudes at each location x
    N : int
        The number of frequencies at which to evaluate the result

    Returns
    -------
    f_hat : ndarray, shape=(N,)
        The amplitudes corresponding to each wave number k = range(-N/2, N/2)

    See Also
    --------
    nfft_adjoint : adjoint non-equispaced fast Fourier transform
    ndft : non-equispaced direct Fourier transform
    nfft : non-equispaced fast Fourier transform
    """
    x, f = np.broadcast_arrays(x, f)
    assert x.ndim == 1

    N = int(N)
    assert N % 2 == 0

    k = -(N // 2) + np.arange(N)

    # Compute phase of the exponential for each value of `k` and `x`.
    # In a similar vein to the FFT, the phase corresponding to `k[0]`
    # and `x[0]` should vanish, the phase corresponding to `k[-1]` and
    # `x[-1]` should be maximal, and so on and so forth.
    ph = 2j * np.pi * (k - k[0]) * ((x - x[0])[:, None])

    # Compute transform. With the above FFT-like ordering in the phase `ph`,
    # we must also remember to shift the zero-frequency component of the
    # spectrum to the center, just as we normally need to do with the FFT.
    fhat = np.dot(f, np.exp(ph))
    fhat = np.fft.fftshift(fhat)

    return fhat

I hope the commentary around my modifications is clear, but please let me know if more details would be helpful or if you have further questions.

Just as you test the equivalence of your ndft and nfft implementations in your walkthrough, I've also verified the equivalence of your ndft_adjoint and nfft_adjoint implementations. Thus, the above phase-computation error also plagues nfft_adjoint. I've glanced through the Keiner, Kunis, & Potts paper, but I can't even pretend that I understand it, so I'm not sure if the fix for nfft_adjoint is as simple as what I've suggested above for ndft_adjoint...

Thanks for starting this package, and I look forward to seeing where it goes!

[jakevdp]

Thanks!

The FFT uses a different frequency convention than the standard form of the DFT, which manifests either as a frequency shift, or a phase shift, or both.

Because this package gives a fast approximation of the DFT, I think sticking with the current form makes the most sense. It might be useful to add some documentation about this difference when it comes to the FFT, though, because it can be a bit confusing if you're not used to the convention used by the FFT.