Accelerating the Nonuniform Fast Fourier Transform

TL;DR: This paper observes that one of the standard interpolation or "gridding" schemes, based on Gaussians, can be accelerated by a significant factor without precomputation and storage of the interpolation weights, of particular value in two- and three- dimensional settings.

Abstract: The nonequispaced Fourier transform arises in a variety of application areas, from medical imaging to radio astronomy to the numerical solution of partial differential equations. In a typical problem, one is given an irregular sampling of N data in the frequency domain and one is interested in reconstructing the corresponding function in the physical domain. When the sampling is uniform, the fast Fourier transform (FFT) allows this calculation to be computed in O(N log N ) operations rather than O(N 2 ) operations. Unfortunately, when the sampling is nonuniform, the FFT does not apply. Over the last few years, a number of algorithms have been developed to overcome this limitation and are often referred to as nonuniform FFTs (NUFFTs). These rely on a mixture of interpolation and the judicious use of the FFT on an oversampled grid (A. Dutt and V. Rokhlin, SIAM J. Sci. Comput., 14 (1993), pp. 1368-1383). In this paper, we observe that one of the standard interpolation or "gridding" schemes, based on Gaussians, can be accelerated by a significant factor without precomputation and storage of the interpolation weights. This is of particular value in two- and three- dimensional settings, saving either 10 d N in storage in d dimensions or a factor of about 5-10 in CPUtime (independent of dimension).

Summary

1. Introduction.

• The authors describe an extremely simple and efficient implementation of the nonuniform fast Fourier transform .
• Let us begin, however, with a more precise description of the computational task.
• There is some confusion in the literature about the use of NUFFTs in this context (see Remark 1 below).
• Not all schemes for reconstructing Fourier integrals of the type (3) can be represented formally as a quadrature of the type (1) .

3. Fast Gaussian Gridding.

• Following standard practice, the authors will refer to these processes as gridding and the M r -point mesh as the oversampled mesh.
• This cost (in either storage or CPU time or both) becomes a significant burden in two, three, and higher dimensions.
• It is sometimes called the curse of dimensionality; in the absence of a separable coordinate system, interpolation-type processes have costs that grow exponentially with dimension.
• This expression for f τ looks much more expensive than it actually is.
• An elementary calculation shows that EQUATION ) Careful organization of the loop shows that, for each source point, two exponential evaluations are required, followed by two multiplications at each of 2M sp regular mesh points.

Example 2 (Fast Gridding Compared to Gridding).

• Both algorithms use the standard FFT on the oversampled mesh, and the time for this step is indicated in Figure 3 of the math coprocessor, cache size, etc.
• Note that the type-1 and type-2 transforms are very similar in terms of floating point operations; the differences in CPU time are due mainly to memory caching issues.
• In any case, the speed-up of the fast gridding algorithm is significant in two dimensions and would be even more significant in the three-dimensional case.
• The authors have not carried out such fine-tuning.

Example 3 (Comparison with Direct Method).

• The authors compare the computational performance of their fast gridding algorithm with direct summation and the standard FFT.
• The direct code was implemented and compiled with the same options: the gcc-2.95 compiler with -O2 optimization on a 450MHz Sparc Ultra-60.
• Table 2 shows how the computational cost for gridding grows with precision.
• In summary, the NUFFT is about 4 times more expensive in one dimension than a traditional M -point FFT for single precision accuracy.
• The MRI hardware is able to acquire the Fourier transform of a particular tissue property at selected points in the frequency domain.

Figures

Fig. 5 Reconstructed image from Fourier data sampled on a radial (polar coordinate) grid. The curves at the right show the reconstructed function along the lines y = 0.25 and y = −0.609.

Table 1

Fig. 2 Gridding error in the NUFFT. In each figure, the x-axis indicates the spreading distance defined by parameter Msp and the y-axis indicates the computed error. The leftmost figure corresponds to a one-dimensional example with 1024 points. The middle example uses a uniform set of data points in two dimensions and the rightmost figure uses a random distribution of data points in two dimensions. Dash-dotted lines show the error for an oversampled grid with Mr = 1.5M , solid lines for Mr = 2M , and dashed lines for Mr = 3M .

Fig. 3 Computing time. Dotted lines represent the cost for FFT, dashed lines for fast gridding, and dash-dotted lines for naive gridding. Heavy dots (·) are used for type-1 and plus marks (+) for type-2 transformation.

Fig. 1 In the version of the NUFFT described here, each delta function source at a point such as xj in (7) is replaced by a Gaussian. This smears the source strength to nearby regular grid points. The regular grid must be fine enough to resolve the smeared function fτ in (8). Note that we include 2π-periodic images of the sources in the definition of the heat kernel gτ . They decay sufficiently rapidly that all but the nearest ones can be ignored.

Table 2

Fig. 4 CPU requirements of the direct code (solid lines with dots), the fast gridding code (four solid lines with tolerance = 10−3, 10−6, 10−9, 10−12 from bottom to top), and the standard FFT for the oversampled mesh with Mr = 2M (dotted lines).

Full text

SIAM REVIEW
c

2004 Society for Industrial and Applied Mathematics
Vol. 46, No. 3, pp. 443–454
Accelerating the Nonuniform
Fast Fourier Transform
∗
Leslie Greengard
†
June-Yub Lee
‡
Abstract. The nonequispaced Fourier transform arises in a variety of application areas, from medical
imaging to radio astronomy to the numerical solution of partial differential equations. In
a typical problem, one is given an irregular sampling of N data in the frequency domain
and one is interested in reconstructing the corresponding function in the physical domain.
When the sampling is uniform, the fast Fourier transform (FFT) allows this calculation
to be computed in O(N log N) operations rather than O(N
2
) operations. Unfortunately,
when the sampling is nonuniform, the FFT does not apply. Over the last few years,
a number of algorithms have been developed to overcome this limitation and are often
referred to as nonuniform FFTs (NUFFTs). These rely on a mixture of interpolation and
the judicious use of the FFT on an oversampled grid [A. Dutt and V. Rokhlin, SIAM J.
Sci. Comput., 14 (1993), pp. 1368–1383].
In this paper, we observe that one of the standard interpolation or “gridding” schemes,
based on Gaussians, can be accelerated by a significant factor without precomputation
and storage of the interpolation weights. This is of particular value in two- and three-
dimensional settings, saving either 10
d
N in storage in d dimensions or a factor of about
5–10 in CPU time (independent of dimension).
Key words. nonuniform fast Fourier transform, fast gridding, FFT, image reconstruction
AMS subject classifications. 42A38, 44A35, 65T50, 65R10
DOI. 10.1137/S003614450343200X
1. Introduction. In this note, we describe an extremely simple and efficient
implementation of the nonuniform fast Fourier transform (NUFFT). There are a
host of applications of such algorithms, and we refer the reader to the references
[2, 6, 8, 11, 13, 14, 17] for examples. We restrict our attention here to one: function
(or image) reconstruction from Fourier data as discussed in [6, 8, 11, 14]. Let us
begin, however, with a more precise description of the computational task. In two
dimensions, we define the nonuniform discrete Fourier transform of types 1 and 2
according to the formulae
F (k
1
,k
2
)=
1
N
N−1

j=0
f
j
e
−i(k
1
,k
2
)·x
j
,(1)
∗
Received by the editors July 23, 2003; accepted for publication (in revised form) December 1,
2003; published electronically July 30, 2004.
http://www.siam.org/journals/sirev/46-3/43200.html
†
Courant Institute of Mathematical Sciences, New York University, New York, NY 10012
(greengard@cims.nyu.edu). The work of this author was supported by the Applied Mathematical
Sciences Program of the U.S. Department of Energy under contract DEFGO288ER25053.
‡
Department of Mathematics, Ewha Womans University, Seoul, 120-750, Korea (jylee@math.
ewha.ac.kr). The work of this author was supported by the Korea Research Foundation under grant
2002-015-CP0044.
443

444 LESLIE GREENGARD AND JUNE-YUB LEE
f(x
j
)=

k
1

k
2
F (k
1
,k
2
) e
i(k
1
,k
2
)·x
j
,(2)
respectively, where x
j
∈ [0, 2π] × [0, 2π] and −
M
2
≤ k
1
,k
2
<
M
2
.
It is, perhaps, convenient to think of (1) as a discretization of the Fourier integral
F (k
1
,k
2
)=
1
(2π)
2

2π
0

2π
0
f(x) e
−i(k
1
,k
2
)·x
dx(3)
with {x
j
} serving as the discretization points. If we let w
j
denote the quadrature
weight corresponding to {x
j
}, then we obtain (1) by setting f
j
= f(x
j
) w
j
. Equation
(2), of course, is simply the evaluation of a finite Fourier series
f(x)=

k
1

k
2
F (k
1
,k
2
) e
i(k
1
,k
2
)·x
(4)
at an arbitrary set of targets.
Nonuniform FFTs of the types discussed here are based, in essence, on combining
some interpolation scheme with the standard FFT. Oddly enough, it was a number
of years after their use in applications before a rigorous analysis of such schemes was
introduced by Dutt and Rokhlin [5]. Subsequent papers, such as [1, 3, 9, 10], described
variants based on alternative interpolation/approximation approaches.
Before discussing the algorithm itself, we would like to comment briefly on appli-
cations that involve evaluation of (3) or the Fourier integral
H(s
1
,s
2
)=

∞
−∞

∞
−∞
h(x) e
−i(s
1
,s
2
)·x
dx(5)
when the transform data h(x) is known at a scatter of points rather than on a regular
Cartesian mesh. There is some confusion in the literature about the use of NUFFTs in
this context (see Remark 1 below). There are three separate issues involved: acquisi-
tion of data h(x
j
) in the Fourier domain, the choice of a quadrature scheme {x
j
,w
j
},
and the availability of a fast algorithm for computing the discrete transform itself.
They are often blended together when they should not be. Dutt and Rokhlin [5] ap-
pear to have been the first to try to isolate one of these problems; they addressed the
algorithmic question and showed that sums of the form (1) or (2) can be computed
in O(N log N) time with complete control of precision. While they concentrated on
the one-dimensional case, higher dimensional versions have been considered by a va-
riety of authors [3, 6, 14]. The first rigorous two-dimensional version can be found
in a paper by Strain [15], which uses the NUFFT to solve a class of elliptic partial
differential equations.
Remark 1. Not all schemes for reconstructing Fourier integrals of the type (3)
can be represented formally as a quadrature of the type (1). Different schemes for
interpolating f(x
j
) to a uniform mesh do give rise to different reconstructed functions
F . However, if the decision has been made to use a quadrature approach such as (1),
then the remaining task is entirely computational. The best algorithm is the one that
evaluates the relevant sums as quickly and as accurately as possible.
2. The NUFFT. Let us first consider the one-dimensional analogs of the summa-
tion problems (1) and (2). With x
j
∈ [0, 2π], the type-1 NUFFT is defined by the

ACCELERATING THE NONUNIFORM FAST FOURIER TRANSFORM 445
calculation of
F (k)=
1
N
N−1

j=0
f
j
e
−ikx
j
for k = −
M
2
,...,
M
2
−1 .(6)
It is based on the following set of observations:
1. Equation (6) describes the exact Fourier coefficients of the function
f(x)=
N−1

j=0
f
j
δ(x − x
j
),(7)
viewed as a periodic function on [0, 2π]. Here, δ(x) denotes the Dirac delta
function. It is clearly not well-resolved by a uniform mesh in x.
2. Let g
τ
(x) denote the one-dimensional periodic heat kernel on [0, 2π], given by
g
τ
(x)=
∞

l=−∞
e
−(x−2lπ)
2
/4τ
.
If we define f
τ
(x) to be the convolution
f
τ
(x)=f ∗ g
τ
(x)=

2π
0
f(y)g
τ
(x − y) dy,(8)
then f
τ
isa2π-periodic C
∞
function and can be well-resolved by a uniform
mesh in x whose spacing is determined by τ (Figure 1). The Fourier coeffi-
cients of f
τ
, namely,
F
τ
(k)=
1
2π

2π
0
f
τ
(x)e
−ikx
dx ,
can be computed with high accuracy using the standard FFT on an oversam-
pled grid
F
τ
(k) ≈
1
M
r
M
r
−1

m=0
f
τ
(2πm/M
r
)e
−ik2πm/M
r
,(9)
where
f
τ
(2πm/M
r
)=
N−1

j=0
f
j
g
τ
(2πm/M
r
− x
j
).(10)
3. Once the values F
τ
(k) are known, an elementary calculation shows that
F (k)=

π
τ
e
k
2
τ
F
τ
(k).(11)
(This is a direct consequence of the convolution theorem and the fact that
the Fourier transform of g
τ
is G
τ
(k)=
√
2τe
−k
2
τ
.)

446 LESLIE GREENGARD AND JUNE-YUB LEE
j
image
image
x
π2
0
Fig. 1 In the version of the NUFFT described here, each delta function source at a point such as
x
j
in (7) is replaced by a Gaussian. This smears the source strength to nearby regular grid
points. The regular grid must be fine enough to resolve the smeared function f
τ
in (8). Note
that we include 2π-periodic images of the sources in the definition of the heat kernel g
τ
.
They decay sufficiently rapidly that all but the nearest ones can be ignored.
Recall that a type-2 transformation evaluates a regular Fourier series at irregular
target points. Thus, in one dimension, for x
j
∈ [0, 2π], the type-2 NUFFT is defined
by the calculation of
f(x
j
)=
M
2
−1

k=−
M
2
F (k) e
ikx
j
.(12)
It is based on a closely related idea:
1. We first deconvolve the Fourier coefficients, defining F
−τ
(k)by
F
−τ
(k)=

π
τ
e
k
2
τ
F (k),(13)
and evaluate the corresponding function f
−τ
(x) on a uniform mesh with M
r
points on [0, 2π] using the FFT
f
−τ
(x)=
M
r
−1

k=0
F
−τ
(k)e
ikx
.(14)
In the preceding expression, we set F
−τ
(k)=0for
M
2
≤ k<M
r
−
M
2
and
view F
−τ
(k)asanM
r
-periodic function F
−τ
(k)=F
−τ
(k − M
r
).
2. We then compute the desired values f(x
k
) from
f(x
j
)=f
−τ
∗ g
τ
(x
j
)=
1
2π

2π
0
f
−τ
(x)g
τ
(x
j
− x) dx
≈
1
M
r
M
r
−1

m=0
f
−τ
(2πm/M
r
) g
τ
(x
j
− 2πm/M
r
).(15)
This is again a direct consequence of the convolution theorem except that we
deconvolve the effect of Gaussian smoothing in (13) before we actually carry
out the smoothing in (15)!

ACCELERATING THE NONUNIFORM FAST FOURIER TRANSFORM 447
Remark 2. There are a number of details that need to be fixed here, including
the selection of τ, the convolution with a Gaussian in both procedures, and the length
M
r
of the FFTs used. We will not repeat the analysis of [5], since the relevant results
can be summarized very simply: with M
r
=2M and τ =12/M
2
, Gaussian spreading
of each source to the nearest 24 grid points yields about 12 digits of accuracy. With
τ =6/M
2
, Gaussian spreading of each source to the nearest 12 grid points yields
about 6 digits of accuracy.
3. Fast Gaussian Gridding. The dominant task in the NUFFT is the calcula-
tion of f
τ
(2πm/M
r
) in (10) and f(x
j
) in (15). Following standard practice, we will
refer to these processes as gridding and the M
r
-point mesh as the oversampled mesh.
In d dimensions, gridding requires 12
d
N exponential evaluations for single precision
accuracy and about 24
d
N exponential evaluations for double precision accuracy. In
order to avoid that cost using existing schemes, one can precompute all the necessary
quantities, incurring a storage cost of 12
d
N or 24
d
N and a computational cost of
12
d
N or 24
d
N multiplications.
This cost (in either storage or CPU time or both) becomes a significant burden in
two, three, and higher dimensions. It is sometimes called the curse of dimensionality;
in the absence of a separable coordinate system, interpolation-type processes have
costs that grow exponentially with dimension. In the remainder of this paper, we
don’t overcome the curse, but we show that (1 + d)N exponential evaluations or
(1+d)N storage is sufficient, followed by (12d+12
d
)N or (24d+24
d
)N multiplications,
depending on the required accuracy. The net reduction in CPU time is by a factor of
5–10.
By inspection of (10), it is evident that we only need values of the function f
τ
at
equispaced points on the oversampled mesh. For this, we have
f
τ
(2πm/M
r
)=
N−1

j=0
f
j
∞

l=−∞
e
−(x
j
−2πm/M
r
−2lπ)
2
/4τ
.
This expression for f
τ
looks much more expensive than it actually is. Since the
Gaussian sources are sharply peaked (in a manner dependent on τ ), each source of
strength f
j
located at x
j
is nonnegligible only at nearby grid points. As mentioned
above, we only need to compute its effect at the nearest 12 or 24 points (Figure 1)
to achieve either 6- or 12-digit accuracy. Thus, for the purpose of computation, we
change our point of view from the receiving point (2πm/M
r
) to the source point x
j
and consider one Gaussian source at a time. An elementary calculation shows that
e
−(x
j
−2πm/M
r
)
2
/4τ
= e
−x
2
j
/4τ

e
x
j
π/M
r
τ

m
e
−(πm/M
r
)
2
/τ
.(16)
But this means that we can compute and store two exponentials, e
−x
2
j
/4τ
and e
x
j
π/M
r
τ
,
for each source point. The third exponential, e
−(πm/M
r
)
2
/τ
, is independent of x
j
.
The gridding algorithm is straightforward:
• Let ξ =2πm/M
r
denote the nearest regular grid point that is less than or
equal to x
j
on the oversampled grid. Beginning at ξ, let M
sp
denote the
number of grid points to which the spreading will be accounted for in each
direction.
• Compute the two exponentials: E
1
= e
−(x
j
−ξ)
2
/4τ
, E
2
= e
(x
j
−ξ)π/M
r
τ
.
• The spreading contribution to f
τ
(2π(m+m

)/M
r
) for −M
sp
<m

≤ M
sp
is
f
j
E
1
·E
m

2
·E
3
(m

), where M
sp
= 6 for single precision, M
sp
= 12 for double
precision, and E
3
(m

)=e
−(πm

/M
r
)
2
/τ
.

Frequently asked questions

Q1. What are the contributions mentioned in the paper "Accelerating the nonuniform fast fourier transform∗" ?

In this paper, the authors observe that one of the standard interpolation or “ gridding ” schemes, based on Gaussians, can be accelerated by a significant factor without precomputation and storage of the interpolation weights. 

Q2. What is the purpose of the MRI?

In most clinical systems, the device is designed to acquire data on a uniform Cartesian mesh, from which a standard FFT can be used for image reconstruction. 

Q3. What is the simplest way to do this?

2. Precompute E3(l) = e−(πl/Mr)2/τ for 0 ≤ l ≤ Msp and E4(k) = E4(M−k) = eτk2 for |k| ≤ M2 .Step C: Convolution for Each Source Point (xj , yj) 1. 

Q4. How many factors affect the computation time of the fast gridding algorithm?

The actual computation time depends on a number of factors, including compiler options, type of CPU, performanceof the math coprocessor, cache size, etc.