Disagreement w/ FFT in limit of uniform grid #6

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()

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