NFFT Package¶

The nfft package is a lightweight implementation of the non-equispaced fast Fourier transform (NFFT), based on numpy and scipy and released under an MIT license.

The NFFT is described in Using NFFT 3 – a software library for various nonequispaced fast Fourier transforms (pdf), which describes a C library that computes the NFFT and several variants and extensions.

Included Algorithms¶

The nfft package currently implements only a few of the algorithms described in the above paper, in particular:

The one-dimensional forward NDFT and NFFT¶

The forward transform is given by

$$ f_j = \sum_{k=-N/2}^{N/2-1} \hat{f}_k e^{-2\pi i k x_j} $$

for complex amplitudes $\{f_k\}$ specified at the range of integer wavenumbers $k$ in the range $-N/2 \le k < N$, evaluated at points $\{x_j\}$ satisfying $-1/2 \le x_j < 1/2$.

This can be computed via the nfft.ndft() and nfft.nfft() functions, respectively.

The one-dimensional adjoint NDFT and NFFT¶

The adjoint transform is given by

$$ \hat{f}_k = \sum_{j=0}^{M-1} f_j e^{2\pi i k x_j} $$

for complex values $\{f_j\}$ at points $\{x_j\}$ satisfying $-1/2 \le x_j < 1/2$, and for the range of integer wavenumbers $k$ in the range $-N/2 \le k < N$.

This can be computed via the nfft.ndft_adjoint() and nfft.nfft_adjoint() functions, respectively.

Complexity¶

The computational complexity of both the forward and adjoint algorithm is approximately

$$ \mathcal{O}[N\log(N) + M\log(1 / \epsilon)] $$

where $\epsilon$ is the desired tolerance of the result. In the current implementation, the memory requirements are approximately

$$ \mathcal{O}[N + M\log(1 / \epsilon)] $$

Comparison to pynfft¶

Another option for computing the NFFT in Python is to use the pynfft package, which wraps the C library referenced in the above paper. The advantage of pynfft is that it provides a more complete set of routines, including multi-dimensional NFFTs and various computing strategies.

The disadvantage is that pynfft is GPL-licensed, and has a more complicated set of dependencies.

Performance-wise, nfft and pynfft are comparable, with the nfft package discussed here being up to a factor of 2 faster in most cases of interest (see Benchmarks.ipynb for some benchmarks).

Installation¶

The nfft package can be installded directly from the Python Package Index:

$ pip install nfft

Dependencies are numpy, scipy, and pytest.

Testing¶

Unit tests can be run using pytest:

$ pytest --pyargs nfft

Usage¶

Here are some examples of computing the NFFT and its adjoint, using both a direct method and the fast method:

Forward Transform¶

$$ f_j = \sum_{k=-N/2}^{N/2-1} \hat{f}_k e^{-2\pi i k x_j} $$

In [1]:

import numpy as np
import nfft

In [2]:

# define evaluation points
x = -0.5 + np.random.rand(1000)

# define Fourier coefficients
N = 10000
k = N // 2 + np.arange(N)
f_k = np.random.randn(N)

In [3]:

# direct Fourier transform
%time f_x_direct = nfft.ndft(x, f_k)

CPU times: user 401 ms, sys: 120 ms, total: 521 ms
Wall time: 505 ms

In [4]:

# fast Fourier transform
%time f_x_fast = nfft.nfft(x, f_k)

CPU times: user 10 ms, sys: 2.11 ms, total: 12.1 ms
Wall time: 7.02 ms

In [5]:

# Compare the results
np.allclose(f_x_direct, f_x_fast)

Out[5]:

True

Adjoint Transform¶

$$ \hat{f}_k = \sum_{j=0}^{M-1} f_j e^{2\pi i k x_j}, $$

In [6]:

# define observations

f = np.random.rand(len(x))

In [7]:

# direct adjoint transform
%time f_k_direct = nfft.ndft_adjoint(x, f, N)

CPU times: user 473 ms, sys: 116 ms, total: 590 ms
Wall time: 450 ms

In [8]:

# fast adjoint transform
%time f_k_fast = nfft.nfft_adjoint(x, f, N)

CPU times: user 7.75 ms, sys: 1.56 ms, total: 9.3 ms
Wall time: 4.94 ms

In [9]:

# Compare the results
np.allclose(f_k_direct, f_k_fast)

Out[9]:

True