%==============================================================================
%\documentclass[11pt,a4paper,twoside,bibtotoc]{scrartcl}
\documentclass[acmtoms]{acmtrans2m}
%==============================================================================
\acmVolume{V}
\acmNumber{N}
\acmYear{08}
\acmMonth{M}

% LaTeX packages
\usepackage{caption}
\usepackage{a4wide}
\usepackage{amsfonts}
\usepackage{amsmath}
\usepackage{theorem}
\usepackage{color}
\usepackage{graphicx}
\usepackage{subfigure}
\usepackage[plain,Alg.]{algorithm}
\usepackage{algorithmic}
\usepackage{multicol}
\usepackage{multirow}
%\usepackage{showkeys}

% Allow for more figures on each side
\renewcommand{\topfraction}{1}
\renewcommand{\textfraction}{0.2}
\setcounter{totalnumber}{4}

% Shortcuts for maths symbols.
\newcommand{\N}{\ensuremath{\mathbb{N}}}
\newcommand{\T}{\ensuremath{\mathbb{T}}}
\renewcommand{\S}{\ensuremath{\mathbb{S}}}
\newcommand{\NZ}{\ensuremath{\mathbb{N}_{0}}}
\newcommand{\Z}{\ensuremath{\mathbb{Z}}}
\newcommand{\R}{\ensuremath{\mathbb{R}}}
\newcommand{\Rp}{\ensuremath{\mathbb{R}^{+}}}
\newcommand{\Rn}{\ensuremath{\mathbb{R}^n}}
\newcommand{\Rnn}{\ensuremath{\mathbb{R}^{n \times n}}}
\newcommand{\C}{\ensuremath{\mathbb{C}}}
\newcommand{\nfft}{\textsf{NFFT}}
\newcommand{\cO}{\ensuremath{\mathcal{O}}}
\newcommand{\tT}{\ensuremath{\text{\tiny{T}}}}
\newcommand{\ti}{\ensuremath{\text{\scriptsize{i}}}}
\newcommand{\e}{{\ensuremath{\mathrm{e}}}}
\newcommand{\eim}[1]{\ensuremath{\mathrm{e}^{-2\pi{\ti} #1}}}
\newcommand{\eip}[1]{\ensuremath{\mathrm{e}^{ 2\pi{\ti} #1}}}
\renewcommand{\mathbf}[1]{\ensuremath{\boldsymbol{#1}}}
\newcommand{\ds}{\displaystyle}
\newcommand{\sinc}{{\rm sinc}}
\newcommand{\dist}{{\rm dist}}
\newcommand{\adj}{{\vdash \hspace*{-1.72mm} \dashv}}
\DeclareMathOperator{\diag}{diag}
\DeclareMathOperator{\dd}{d}
\newcommand{\bend}
{
  \hspace*{0ex} \hfill
  \hbox{\vrule height 1.5ex \vbox{\hrule width 1.4ex \vskip
  1.4ex\hrule  width 1.4ex} \vrule height 1.5ex}
}

% Environment Numbering
\numberwithin{equation}{section}
\numberwithin{table}{section}
\numberwithin{figure}{section}

\DeclareCaptionLabelSeparator{periodem}{.\hskip 1em}
\captionsetup{font=footnotesize,labelfont=,labelsep=periodem}

%\captionnamefont{\footnotesize}
%\captiontitlefont{\footnotesize}

%===============================================================================
% Title Page
%===============================================================================
\markboth{Keiner, Kunis, Potts}{NFFT, nonequispaced fast Fourier transforms}

\title{Using NFFT~3 -- a software library for various nonequispaced fast
  Fourier transforms}
\author{
  Jens Keiner\\
  University of L\"ubeck\\
  23560 L\"ubeck, Germany\\
  keiner@math.uni-luebeck.de
  \and
  Stefan Kunis\\
  Chemnitz University of Technology\\
  09107 Chemnitz, Germany\\
  kunis@mathematik.tu-chemnitz.de
  \and
  Daniel Potts\\
  Chemnitz University of Technology\\
  09107 Chemnitz, Germany\\
  potts@mathematik.tu-chemnitz.de
  }
\begin{abstract}
  NFFT~3 is a software library that implements the nonequispaced fast Fourier
  transform (NFFT) and a number of related algorithms, e.g. nonequispaced fast
  Fourier transforms on the sphere and iterative schemes for inversion. This is
  to provide a survey on the mathematical concepts behind the NFFT and its variants,
  as well as a general guideline for using the library. Numerical examples for
  a number of applications are given.
\end{abstract}

\category{G.1}{F.2.1}{Numerical analysis, Computation of transforms}
\terms{Algorithms, Documentation, Theory}
\keywords{fast Fourier transforms, approximative algorithms}

\begin{document}
\maketitle
%===============================================================================
\section{Introduction} \label{sect:1}
%===============================================================================
\setcounter{equation}{0}

NFFT~3 \cite{nfft3} is a C software library for computing
the nonequispaced fast Fourier transform (NFFT); see Appendix \ref{A.4} for
alternative denominations, further approaches, and references.
At the heart of the library is a fast algorithm that generalizes the
ubiquitous FFT \cite{CoTu65} from equally spaced to arbitrary sampling
points/spatial nodes, hence the name NFFT; see Appendix \ref{appendix.nfft}.
Our implementation is based on the popular FFTW library \cite{FFTW05}.
Unlike the FFT, the NFFT is not trivially inverted, which compells one to
resort to iterative techniques. Our library also implements these iterative
algorithms for inversion, as well as several other discrete transforms related to
or based on the NFFT. This paper is to present an overview over NFFT~3, the
underlying algorithmic concepts, and intents to serve as a general guideline.

We start in Section \ref{sect:nfft} with the definition of the nonequispaced
discrete Fourier transform (NDFT) as a generalization of the discrete
Fourier transform (DFT) to arbitrary nodes. The NFFT algorithm, its key properties,
and differences to the classical FFT are discussed.
Having settled the relevant definitions, the basic principles
of using the NFFT~3 software library are described in Section \ref{sect:lib}.
%shows how to compile and install the software package on
%your computer system.
%A small example in one dimension demonstrates the basic principles.
An overview over the software interface to the NFFT library routines concludes
this section.
In Section \ref{sect:gen}, several transforms that derive from the original NFFT
are described. This includes, among others, the NFFT on the unit sphere and
iterative schemes for inversion of the NFFT. Section \ref{sect:sum} is to summarize
the theoretical arithmetic costs of our algorithms and presents actual
performance measurements for our implementation.
In the final Section \ref{sect:exap}, a number of elementary numerical examples
as well as more advanced applications are given.

We have restricted ourselves to describe brief recipes which should be
sufficient for most users. For the experienced user, a comprehensive software
interface documentation is provided within the downloadable software package
\cite[doc/refman.pdf]{nfft3}.

%===============================================================================
\section{Nonequispaced discrete Fourier transform} \label{sect:nfft}
%===============================================================================
\setcounter{equation}{0}

Fast Fourier transforms (FFTs) have recently been generalized to
arbitrary sampling situations; see
\cite{duro93,bey95,AnDa96,st97,ware98,postta01} and the references in Appendix
\ref{A.4}.
A convenient approach is to define the discrete Fourier transform (DFT) as the
evaluation of a trigonometric polynomial, given by its Fourier coefficients,
at equispaced spatial nodes. This can then be generalized to arbitrary nodes.

\subsection{Discrete Fourier transform}
Following standard conventions \cite{loan,FFTW05} the (forward) discrete Fourier
transform (DFT) is defined as the calculation of the sums
\begin{equation} \label{eq.dft}
  f_{j} = \sum_{k = 0}^{N-1} \hat{f}_{k} \mathrm{e}^{-2\pi\mathrm{i}jk/N} \qquad
  (j=0,\ldots,N-1)
\end{equation}
for $N \in \mathbb{N}$ and given coefficients $\hat{f}_{k}\in\C$.
The transform is readily inverted, i.e.,
\begin{equation}
  \label{eq.idft}
  \hat{f}_{k} = \frac{1}{N} \sum_{j = 0}^{N-1} f_{j} \mathrm{e}^{2\pi\mathrm{i}kj/N}.
\end{equation}
Clearly, the computation of these transforms requires $\mathcal{O}(N^{2})$
arithmetic operations in general. Fast Fourier transforms (FFTs) are
$\mathcal{O}(N \log N)$ algorithms to compute the same result.

In matrix-vector notation, the DFT can be written as the matrix-vector product
\begin{align*}
  \mathbf{f} = \mathbf{F} \, \mathbf{\hat{f}}, \quad
\mathbf{F} =
(\mathrm{e}^{-2\pi\mathrm{i}jk/N})_{j,k=0,\ldots,N}
\end{align*}
with vectors $\mathbf{f} = (f_j)_{j=0,\ldots,N-1}$,
$\mathbf{\hat{f}} = (\hat{f}_k)_{k=0,\ldots,N-1}$,
and the Fourier matrix $\mathbf{F}$.
%= (\mathrm{e}^{-2\pi\mathrm{i}jk/N})_{j,k=0,\ldots,N}$.
The inverse transform is equivalent to $\mathbf{\hat{f}} = \tfrac{1}{N}
\mathbf{F}^{\adj} \, \mathbf{f}$.
It is straightforward to generalize these concepts to multidimensional
transforms; see Section \ref{subsect.definition} below. But before we
define the NDFT, some notation is at order.

\subsection{Notation}
A $d$-dimensional nonequispaced discrete Fourier transform is defined by a set
of arbitrary spatial nodes $\mathcal{X}$ and a frequency bandwidth vector
$\mathbf{N}$.
Each node $\mathbf x_j$ in the sampling set
${\cal X}:=\{\mathbf{x}_j \in \T^d:\,j=0,\hdots,M-1\}$ is drawn from the
$d$-dimensional torus $\mathbb{T}^d \cong [-\tfrac{1}{2},\tfrac{1}{2})^d$
with the number of nodes equal to $|\mathcal{X}| = M$.

For each dimension $t = 0,\ldots,d-1$, the bandwidth $N_{t} \in 2\N$ is
defined as a fixed even number. These bandwidths are collected in the
vector $\mathbf{N} := \left({N_0},\ldots,N_{d-1}\right)^{\top}$.
Let us define the multi-index set
\begin{equation*}
  I_{\mathbf{N}} :=
  \mathbb{Z}^d \cap \prod_{t=0}^{d-1} \left[-\tfrac{N_t\rule[-0.4ex]{0ex}{1ex}}{2},
\tfrac{\rule[-0.4ex]{0ex}{1ex}N_t}{2}\right)
=\left\{ \mathbf{k}=\left(k_t\right)_{t=0,\hdots,d-1} \in
  \Z^d: - \frac{N_t}{2} \le k_t < \frac{N_t}{2} ,\;t=0,\hdots,d-1\right\}.
\end{equation*}
Then the set $I_{\mathbf{N}}$ is a representation of all possible frequencies
in a transform. For simplicity, we use the convention that a multi-index
$\mathbf{k}$ can also reference elements in vectors or rows and columns of matrices
as if it were linearized to $k:=\sum_{t=0}^{d-1}
(k_t+\frac{N_t}{2}) \prod_{t'=t+1}^{d-1} N_{t'}$.

The inner product between a frequency index $\mathbf{k}$ and a node $\mathbf{x}$
is defined as usual, $\mathbf{k} \mathbf{x}:=k_0 x_0 +\hdots+
k_{d-1} x_{d-1}$.
Furthermore, we define the component-wise vector product $\mathbf{\sigma} \odot
\mathbf{N}:=\left(\sigma_0 N_0, \hdots, \sigma_{d-1} N_{d-1}, \right)^{\top}$,
and the reciprocal of a vector $\mathbf{N}$ with nonzero components, $\mathbf{N}^{-1}:=
\left({N_0}^{-1},\ldots,{N_{d-1}^{-1}}\right)^{\top}$.

The class of functions $f: \T^d \rightarrow \C$ that is naturally associated
with an index set $I_{\mathbf{N}}$ is the space of $d$-variate, one-periodic
trigonometric polynomials of degree at most $N_t$ along each dimension $t$,
$
  T_{\mathbf{N}}:={\rm
span}\left(\eim{\mathbf{k}\,\mathbf{\cdot}\,}:\,\mathbf{k}
  \in I_{\mathbf{N}}\right)
$.
This space has dimension $\dim T_{\mathbf{N}} = |I_{\mathbf{N}}| = N_{0} \cdot
\ldots \cdot N_{d-1}$.

%-------------------------------------------------------------------------------
\subsection{Definition}
\label{subsect.definition}
%-------------------------------------------------------------------------------
We are now ready to define the \emph{nonequispaced discrete Fourier transform (NDFT)}
as a generalization of the forward DFT \eqref{eq.dft} to arbitrary nodes.
Given Fourier coefficients $\hat f_{\mathbf{k}} \in \C$,
$\mathbf{k}\in I_{\mathbf{N}}$, as input, the NDFT is defined as the evaluation
of the corresponding trigonometric polynomial $f \in T_{\mathbf{N}}$ at the
set of $M$ arbitrary nodes $\mathcal{X}$, i.e., the calculation of the sums
\begin{equation}
  \label{eq:nfft}
  f_j = \sum_{\mathbf{k}\in I_{\mathbf{N}}} \hat{f}_{\mathbf{k}}
  \eim{\mathbf{k}\mathbf{x}_j} \qquad (j=0,\hdots,M-1).
\end{equation}
To see how this compares with the definition of the DFT, take
$\mathbf{N}=(N,\hdots,N)^{\top}$, $N\in 2\N$, and $N^d$
equispaced nodes $\mathbf{x}_{\mathbf{j}}={\frac{1}{N}}\mathbf{j}$,
$\mathbf{j} \in I_{\mathbf{N}}$. Then \eqref{eq:nfft}
reduces (up to an index shift) to a usual multidimensional forward
DFT.

The NDFT \eqref{eq:nfft} too can be written as a matrix-vector product,
\begin{equation*}
  \mathbf{f} = \mathbf{A} \mathbf{\hat{f}},
\end{equation*}
with the vectors $\mathbf{f} := \left(f_j\right)_{j=0,\hdots,M-1}$,
$\mathbf{\hat{f}} :=(\hat{f}_{\mathbf{k}})_{\mathbf{k}\in
I_{\mathbf{N}}}$, and the nonequispaced Fourier matrix $\mathbf{A}
:= \left(\eim{\mathbf{k}\mathbf{x}_j}\right)_{j=0,\hdots,M-1;\;\mathbf{k}\in
I_{\mathbf{N}}}$.
Typically, the matrix $\mathbf{A}$ is not square. Even if this should be the case,
it is usually neither orthogonal nor has an inverse. Therefore,
the definition of an inverse NDFT transform is not canonical.
Instead, it is customary to define the \emph{adjoint NDFT} by
the matrix-vector product
\[
  \mathbf{\hat{h}} = \mathbf{A}^{\adj} \mathbf{f}.
\]
This is equivalent to the sums
\begin{equation}
  \label{eq:nfftH}
  \hat{h}_{\mathbf{k}}  = \sum_{j=0}^{M-1} \hat{f}_{ j} \eip{\mathbf{k}
  \mathbf{x}_j} \qquad (\mathbf{k}\in I_{\mathbf{N}}).
\end{equation}

The NFFT algorithm, implemented in our library, is a fast approximate algorithm
to compute the sums in \eqref{eq:nfft}. It contains also a fast algorithm to
compute the adjoint transform \eqref{eq:nfftH}.

%-------------------------------------------------------------------------------
\subsection{Key features}
%-------------------------------------------------------------------------------
A $d$-dimensional DFT of degree $N$ along each dimension has $N^d$
input and output coefficients. This generally requires $\mathcal{O}(N^{d+1})$
to compute the transform. FFTs reduce this to $\mathcal{O}(N^{d} \log N)$.

The cost for nonequispaced fast transforms is similar. However, it is important
to note that there the number of nodes as well as their configuration are arbitrary.
This generally implies $M \neq |I_{\mathbf{N}}|$.
The direct calculation of the sums in \eqref{eq:nfft}, an algorithm which we
also denote NDFT, requires $\cO(M |I_{\mathbf{N}}|)$ floating point operations.
%and uses no additional scratch memory.
Since the number of nodes $M$ is typically of order $|I_{\mathbf{N}}|$, this cost
is considered too expensive for most applications.
The NFFT~3 library implements a fast approximate algorithm, called
NFFT, which computes the same result using only
${\cal O}\left(|I_{\mathbf{N}}| \log |I_{\mathbf{N}}|  + |\log \varepsilon|^d
M\right)$ operations.
Here, $\varepsilon$ is the desired accuracy of the computation. The
constant contained in this notation depends on the amount of extra
memory spent for precomputation; see Appendix \ref{sect:windows}.

Let us briefly compare NDFT/NFFT and classical DFT/FFT to observe some
key features and differences:
While a multidimensional DFT/FFT is uniquely defined by the bandwidth
vector $\mathbf{N}$ alone, i.e., the sampling set $\mathcal{X}$ is fixed,
the NDFT/NFFT also depend on the actual choice of $\mathcal{X}$;
see Subsection \ref{subsect.definition}.
FFT algorithms are usually exact algorithms with small errors in the result
caused only by limitations of floating point arithmetic.
The NFFT algorithm purposely introduces a systematic error in the computation
to achieve its favorable arithmetic cost.
This additional error can be controlled, and, if deemed necessary, can be
reduced to the order of machine precision.
To this end, the NFFT uses an \emph{oversampled} FFT internally (in our case
we use routines from FFTW) together with a pre- and post-processing step.
%While the pre-processing step is exact, the post-processing step contains a
%\emph{truncated} approximation scheme.
This introduces two additional parameters,
an \emph{oversampling factor} and a \emph{truncation parameter}, that
control the accuracy of the NFFT. More details are found in Appendix
\ref{appendix.nfft}.

%\textbf{Dem ersten Referee folgend, koennten hier mehr Details zum Algorithmus
%stehen, insbesondere welche zusaetzlichen Parameter der Benutzer jetzt
%einstellen kann. Gut waeren sicher Verweise auf Stellen weiter hinten, die auf
%einzelne Punkte genauer eingehen. Ich denke, eine gute Erklaerung der
%zusaetzlichen Freiheitsgrade und Parameter ist sehr wichtig fuer das
%Verstaendnis.}
%For now, it suffices to note that the NFFT algorithm has a similar complexity
%than the classical FFT algorithm but that it is an approximate algorithm whose
%accuracy is controlled by the user.

%===============================================================================
\section{Guidelines for using the library} \label{sect:lib}
%===============================================================================
NFFT~3 is a free software library completely written in C.
The downloadable source code comes as a GNU-style package and follows standard
conventions. Brief installation instructions and a source tree can be found in
Appendix \ref{sect:inst}.
To compile and link programmes that use NFFT~3, you also need the FFTW library
\cite{fftw} in version 3.0.0 or better.

%-------------------------------------------------------------------------------
\subsection{Interface} \label{sect:general_proc}
%-------------------------------------------------------------------------------
The interface of NFFT~3 resembles that of FFTW 3. A transform is
specified by creating a \emph{plan} -- a data structure that contains the
information to determine the characteristics of a transform. In contrast to
FFTW 3, it is also necessary to call a precomputation routine, once the
nodes $\cal X$ have been provided. After these initial steps, transforms are
executed on input data by invoking an execution routine.
At the end of its lifetime, a plan is finalized and all associated data
structures are destroyed.

%-------------------------------------------------------------------------------
\subsection{Typical workflow} \label{sect:general_ex}
%-------------------------------------------------------------------------------
This brief example is to illustrate the basic principles. It shows how to
compute a one-dimensional NFFT.

%-------------------------------------------------------------------------------
\subsubsection*{Creating a plan}
%-------------------------------------------------------------------------------
An NFFT plan \verb"p" is declared by
\begin{verbatim}
  nfft_plan p;
\end{verbatim}
The simplest way to initialize the plan for a one-dimensional transform is to
call
\begin{verbatim}
  nfft_init_1d(&p,N0,M);
\end{verbatim}
The arguments to \verb"nfft_init_1d" are the pointer \verb"&p" to the yet
uninitialized plan, the number of Fourier
coefficients \verb"N0", and the number of nodes \verb"M".
This creates and initializes all data structures inside the
plan. Memory for Fourier coefficients, nodes, and
function values is allocated automatically.

%But to completely specify the transform, we also have to provide the
%nodes.

%-------------------------------------------------------------------------------
\subsubsection*{Providing the nodes}
%-------------------------------------------------------------------------------
The nodes have to be stored in the member array \verb"p.x".
In one dimension, the coordinate of the $j$th node is a number in
$[-\tfrac{1}{2},\tfrac{1}{2})$ that is stored at the index $j$, i.e.,\\[2ex]
%\begin{verbatim}
\verb"  p.x[j] = "\textsl{your choice in [-0.5,0.5)}\verb";"
%\end{verbatim}

%-------------------------------------------------------------------------------
\subsubsection*{Precomputation}
%-------------------------------------------------------------------------------
The precomputation procedure for the plan is invoked (conditionally) by
\begin{verbatim}
  if (p.nfft_flags & PRE_ONE_PSI)
    nfft_precompute_one_psi(&p);
\end{verbatim}
NFFT~3 uses different strategies here which are free of choice by the user via
passing the appropriate flag to one of the more advanced plan initialization
routines. A default strategy is chosen by the simple initialization routines.
For convenience, the flag \verb+PRE_ONE_PSI+ indicates that precomputation
is necessary; see Appendix \ref{sect:windows} for details.

%-------------------------------------------------------------------------------
\subsubsection*{Providing the Fourier coefficients}
%-------------------------------------------------------------------------------
Fourier coefficients are stored in the one-dimensional member array
\verb"p.f_hat". In the one-dimensional case this array contains the Fourier
coefficients in the order $\hat{f}_{-N/2},\hdots,\hat{f}_{N/2-1}$.
They are accessed by\\[2ex]
%\begin{verbatim}
\verb"  p.f_hat[k] = "\textsl{your choice in the complex numbers}\verb";"
%\end{verbatim}

%-------------------------------------------------------------------------------
\subsubsection*{Execution}
%-------------------------------------------------------------------------------
To compute \eqref{eq:nfft} with the NFFT algorithm, invoke
\begin{verbatim}
  nfft_trafo(&p);
\end{verbatim}
The $M$ output values $f_j$ are stored in the one-dimensional member array
\verb"p.f". To compute a second NFFT with different input data, simply
update the Fourier coefficients in \verb"p.f_hat" and invoke \verb"nfft_trafo"
a second time.

The adjoint transform \eqref{eq:nfftH} is computed by \verb+nfft_adjoint+. Here,
the input is put in \verb"p.f" and the output is written to \verb"p.f_hat".

%-------------------------------------------------------------------------------
\subsubsection*{Finalization}
%-------------------------------------------------------------------------------
Since the initialization routine for the plan has allocated memory, it is
essential to finalize the plan, once it is no longer
needed. This is done by
\begin{verbatim}
  nfft_finalize(&p);
\end{verbatim}

%-------------------------------------------------------------------------------
\subsubsection*{Example}
%-------------------------------------------------------------------------------
Here is a brief example that computes a one-dimensional NFFT with 14
Fourier coefficients and 19 nodes. The routine is part of the example program
\verb"simple_test.c" which is found under \verb"examples/nfft/" in the
software package. It uses some auxilliary routines from the NFFT~3 header
\verb"util.h" to initialize vectors with random entries and to write them
to the standard output.
\begin{verbatim}
  void simple_test_nfft_1d()
  {
    nfft_plan p; int N=14; int M=19;
    nfft_init_1d(&p,N,M);
    nfft_vrand_shifted_unit_double(p.x,p.M_total);
    if(p.nfft_flags & PRE_ONE_PSI) nfft_precompute_one_psi(&p);
    nfft_vrand_unit_complex(p.f_hat,p.N_total);
    nfft_vpr_complex(p.f_hat,p.N_total,"given Fourier coefficients, f_hat");
    ndft_trafo(&p);
    nfft_vpr_complex(p.f,p.M_total,"ndft, f");
    nfft_trafo(&p);
    nfft_vpr_complex(p.f,p.M_total,"nfft, f");
    nfft_finalize(&p);
  }
\end{verbatim}

%-------------------------------------------------------------------------------
\subsection{Multivariate transforms}
%-------------------------------------------------------------------------------
In $d$ dimensions, the nodes $\mathbf{x}_j$ in the sampling set
$\mathcal{X}$, as well as the frequency indices $\mathbf{k}$ in the set
$I_{\mathbf{N}}$ are $d$-dimensional vectors. This affects the storage layout of
the arrays \verb"x" and \verb"f_hat". Generally, both remain
one-dimensional and all $d$-dimensional quantities are linearized.
For \verb"x" this means that the $t$th coordinate of the $j$th node
($t=0,\ldots,d-1$; $j = 0,\ldots,M-1$) is defined via\\[\parskip]

\noindent\verb"  p.x[d*j+t] = "\textsl{your choice in [-0.5,0.5)}\verb";"\\[\parskip]

\noindent This linear ordering has been chosen because the node configuration
may be completely unstructured.
Analogously, the Fourier coefficient $\hat{f}_{\mathbf{k}}$ is stored in
\verb"f_hat[k]" at the linearized index \verb"k"$ = \sum_{t=0}^{d-1}
(k_t+\frac{N_t}{2}) \prod_{t'=t+1}^{d-1} N_{t'}$.

Due to software engineering considerations, we do not offer alternative data
layouts, even if the nodes belong to a multidimensional structure, e.g. a
tensor product grid. Moreover, NFFT~3 does currently not offer support for strided
data or even more general FFTW-style i/o-tensors.

%-------------------------------------------------------------------------------
\subsection{Programming interface and data structures}
%-------------------------------------------------------------------------------
The programming interface for the NFFT module is summarized in
Table \ref{table:nfft_user_fun}. To allow for a comparison with the NDFT algorithm,
the routine \verb"ndft_trafo" computes the transform via direct
evaluation of \eqref{eq:nfft}. The routines \verb"ndft_adjoint"
and \verb"nfft_adjoint" compute adjoint transforms \eqref{eq:nfftH},
respectively. The Initialization routines \verb"nfft_init_1d",
\verb"nfft_init_2d", \verb"nfft_init_3d" allow for convenient initialization
of plans for one-, two-, and three-dimensional transforms with default
parameters. Plans for general multivariate transforms are created using
\verb"nfft_init". For the experienced user, the function \verb"nfft_init_guru"
allows control over all relevant parameters.

The most important members of the \verb"nfft_plan" data type are listed in Table
\ref{table:nfft_members}. Only \verb"f_hat", \verb"f", and
\verb"x" are to be modified by the user.

\begin{acmtable}{0.72\textwidth}
  \centering%
    \vspace{1.3ex}%
    \begin{tabular}{l}
      \verb"void ndft_trafo(nfft_plan* p)" \\
      \verb"void ndft_adjoint(nfft_plan* p)" \\
      \verb"void nfft_trafo(nfft_plan* p)" \\
      \verb"void nfft_adjoint(nfft_plan* p)" \\
      \verb"void nfft_init_1d(nfft_plan* p, int N0, int M)" \\
      \verb"void nfft_init_2d(nfft_plan* p, int N0, int N1, int M)" \\
      \verb"void nfft_init_3d(nfft_plan* p, int N0, int N1, int N2, int M)" \\
      \verb"void nfft_init(nfft_plan* p, int d, int* N, int M)" \\
      \verb"void nfft_init_guru(nfft_plan* p, int d, int* N, int M, int* n," \\
      \verb"                    int m, unsigned nfft_flags, unsigned fftw_flags)" \\
      \verb"void nfft_precompute_one_psi(nfft_plan* p)" \\
      \verb"void nfft_check(nfft_plan* p)" \\
      \verb"void nfft_finalize(nfft_plan* p)"
    \end{tabular}%
    \vspace{0.5ex}%
  \caption{Interface to NDFT/NFFT transforms in NFFT~3.}%
  \label{table:nfft_user_fun}%
\end{acmtable}

\begin{acmtable}{0.72\textwidth}
  \centering%
    \vspace{1.3ex}%
    \begin{tabular}{r|l|l|l}%
      Type & Name & Size & Description
      \\
      \hline
      \verb"int" & \verb"d" & 1 & Spatial dimension $d$
      \\
      \verb"int*" & \verb"N" & $d$ & Multibandwidth $\mathbf{N}$
      \\
      \verb"int" & \verb"N_total" & 1 & Number of coefficients $|I_{\mathbf{N}}|$
      \\
      \verb"int" & \verb"M_total" & 1 & Number of nodes $M$
      \\
      \verb"double complex*" & \verb"f_hat" & $|I_{\mathbf{N}}|$ & Fourier coefficients $\mathbf{\hat{f}}$ or
      \\
      & & & adjoint coefficients $\mathbf{\hat{h}}$
      \\
      \verb"double complex*"& \verb"f" & $M$ & Samples $\mathbf{f}$
      \\
      \verb"double*" & \verb"x" & $dM$ & Sampling set ${\cal X}$
    \end{tabular}%
    \vspace{0.5ex}%
  \caption{Most important members of the structure \texttt{nfft\_plan}.}
  \label{table:nfft_members}
\end{acmtable}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%
%\pagebreak
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%

%===============================================================================
\section{Generalizations} \label{sect:gen}
%===============================================================================

This section surveys further generalizations of the NFFT algorithm that are
implemented by the NFFT~3 library. The interfaces resemble that of the NFFT
module described above, with the prefix \verb+nfft_+ replaced by the
acronym related to the specific transform in all names. Details can be found in
the API documentation accompanying the library. All algorithms, with exception of
the discrete polynomial transform, are based on the same computational procedure
as the NFFT.

%-------------------------------------------------------------------------------
\subsection{NFST/NFCT - Nonequispaced fast sine/cosine transforms}
%-------------------------------------------------------------------------------
Discrete (co)sine transforms (DCT/DST) are closely related real variants of the
discrete Fourier transform DFT. In the same way, nonequispaced discrete (co)sine
transforms (NDCT/NDST) are real transforms related to the NDFT. They are
parameterised by arbitrary nodes $\mathbf{x}_j \in [0,\frac{1}{2}]^d$.
Frequencies $\mathbf{k}$ are taken from one of the index sets
\begin{align*}
  I_{\mathbf{N}}^C &:= \left\{ \mathbf{k}=\left(k_t\right)_{t=0,\hdots,d-1}
    \in \Z^d: 0 \le k_t < N_t ,\;t=0,\hdots,d-1\right\}, \\
  I_{\mathbf{N}}^S &:= \left\{ \mathbf{k}=\left(k_t\right)_{t=0,\hdots,d-1} \in
    \Z^d: 1 \le k_t < N_t ,\;t=0,\hdots,d-1\right\},
\end{align*}
where $I_{\mathbf{N}}^C$ corresponds to the NDCT and $I_{\mathbf{N}}^S$ to the
NDST, respectively. For given real Fourier coefficients
$\hat f_{\mathbf{k}} \in \R$, NDCT and NDST are defined by
\begin{equation*}
  f_j  = \sum_{\mathbf{k}\in I_{\mathbf{N}}^C}
\hat{f}_{\mathbf{k}}
  \cos(2\pi(\mathbf{k}\odot\mathbf{x}_j))\quad
  (j=0,\hdots,M-1),
  \qquad
  f_j  = \sum_{\mathbf{k}\in I_{\mathbf{N}}^S}
\hat{f}_{\mathbf{k}}
  \sin(2\pi(\mathbf{k}\odot\mathbf{x}_j)),
\end{equation*}
respectively. Here, for notational convenience, we have defined
$\cos(\mathbf{x}) := \cos(x_0) \cdot \hdots \cdot \cos(x_{d-1})$ and
$\sin(\mathbf{x}) := \sin(x_0) \cdot \hdots \cdot \sin(x_{d-1})$.

Straightforward algorithms for these transforms need $\cO(M
|I_{\mathbf{N}}^C|)$ and $\cO(M |I_{\mathbf{N}}^S|)$ arithmetic operations.
The fast algorithms (NFST/NFCT) need $\cO(|I_{\mathbf{N}}^S| \log
|I_{\mathbf{N}}^S| + |\log \varepsilon|^d M)$ and $\cO(|I_{\mathbf{N}}^C| \log
|I_{\mathbf{N}}^C| + |\log \varepsilon|^d M)$ operations, respectively.
Note that for these real transforms, the adjoint variants coincide with the
transposed versions. Details about the implemented algorithms can be found in
\cite{tili00,po01,fepo02}.

%-------------------------------------------------------------------------------
\subsection{NSFFT - Nonequispaced sparse fast Fourier transform}
%-------------------------------------------------------------------------------
In higher dimensions, the so-called ``curse of dimensionality''
\cite{Zeng91,Spr00,BuZe04} often requires the use of Fourier series where
frequency indices are restricted to the hyperbolic cross
\begin{equation*}
  H_N^d := \bigcup_{\mathbf{N}\in\N^d,\, |I_{\mathbf{N}}| = N} I_{\mathbf{N}},
  \qquad N=2^{J+2},\;J\in\NZ,
\end{equation*}
which has only $\cO(N \log^{d-1} N)$ degrees of freedom compared to
$\cO(N^d)$ in the plain case. If compared to ordinary Fourier series,
the approximation error commited, which is typically measured in a norm of
dominated mixed smoothness, can be shown to get worse by a factor of only
$\log^{d-1} N$ (see \cite{Spr00}).
The \emph{nonequispaced sparse discrete Fourier transform (NSDFT)} is defined as
the evaluation of
\begin{equation*}
  f_j  = \sum_{\mathbf{k}\in H_N^d}
  \hat f_{\mathbf{k}} \eim{\mathbf{k} \mathbf{x}_j} \qquad (j=0,\ldots,M-1)
\end{equation*}
with given nodes $\mathbf{x}_j\in \T^d$ and Fourier coefficients
$\hat f_{\mathbf{k}}\in\C$. The algorithms implemented in NFFT~3 reduce the
number of arithmetic operations from $\cO(M N\log^{d-1}N)$ to
$\cO(N\log^2 N+|\log \varepsilon|^2 M)$ in two dimensions ($d = 2$), and to
$\cO(N^{3/2} \log N+ |\log \varepsilon|^3 M)$ in three dimensions ($d = 3$).
Details are found in \cite{FeKuPo} for details.

%-------------------------------------------------------------------------------
\subsection{NNFFT - fast Fourier
transform for nonequispaced data in space and frequency domain}
%-------------------------------------------------------------------------------
The NNFFT is a variant of the NFFT where also frequencies are no
longer equispaced integer multi-indices.
The transform is defined by
\begin{equation*}
  f_j  = \sum_{l=0}^{L-1} \hat f_l \eim{\left(\mathbf{v}_l \odot
  \mathbf{N}\right)\mathbf{x}_j} \qquad (j =0,\hdots,M-1)
\end{equation*}
with Fourier coefficients $\hat{f}_l \in\C$, and arbitrary frequencies
$\mathbf{v}_l \in \T^d$ and nodes $\mathbf{x}_j \in \T^d$.
Here, $\mathbf{N}\in \N^d$ is called the \emph{nonharmonic
bandwidth}. This transform is known as \textit{fast Fourier
transform for nonequispaced data in space and frequency domain (NNFFT)}
\cite{duro93,elst,postta01} or as \emph{type 3 nonuniform FFT} \cite{LeGr05}
and takes $\cO(|I_{\mathbf{N}}|\log|I_{\mathbf{N}}| + |\log
\varepsilon|^d(L+M))$ arithmetic operations to compute.

%------------------------------------------------------------------------------
\subsection{NFSFT - nonequispaced fast spherical Fourier transform}
\label{sect:gen:nfsft}
%------------------------------------------------------------------------------
Denote by $\S^2:= \big\{\mathbf{x} \in \R^{3}:\;\|\mathbf{x}\|_2=
1\big\}$ the two-dimensional unit sphere.
In spherical cordinates, any point on $\S^2$ is identified with a pair of
angles $(\vartheta,\varphi)^{\top}\in  [0,\pi]\times [0,2\pi)$.
Hence, a sampling set on the sphere can be defined by
$
  \mathcal{X}: = \left\{\left(\vartheta_{j},\varphi_{j}\right)
  %\in \mathbb{S}^2
  :\,j=0,\hdots,M-1\right\}.
$
A convenient basis of orthogonal functions on $\S^2$ are spherical harmonics
$Y_k^n : \S^2 \rightarrow \C$, $k \in \mathbb{N}_{0}$, $|n| \le k$,
defined by
\begin{equation*}
  Y_k^n(\vartheta,\varphi) := P_k^{|n|}(\cos\vartheta)
  \, \mathrm{e}^{\mathrm{i} n \varphi}.
\end{equation*}
Here, $P_k^{|n|}$ are the associated Legendre
functions of degree $k$ and order $|n|$ (see \cite[pp. 331]{abst}).
Frequencies for a finite transform up to bandwidth $N$ are collected in the
index set
\[
J_{N} := \{(k,n)\in\Z^2:\,k=0,\ldots,N;\; n=-k,\ldots,k\}.
\]
The \emph{nonequispaced discrete spherical Fourier transform (NDSFT)} is
defined as the calculation of
\begin{equation*}
  f_j  =  \sum_{(k,n) \in \, J_{N}} \hat{f}_k^n \,
  Y_k^n\left(\vartheta_{j},\varphi_{j}\right) \qquad (j=0,\ldots,M-1).
\end{equation*}

The direct method (NDSFT) clearly needs $\cO(M N^2)$ arithmetic operations.
The approximate $\cO(N^2 \log^2 N + |\log{\varepsilon}|^2M)$ algorithm
implemented in the NFFT~3 library (\emph{NFSFT}) is a combination of the
fast polynomial transform (FPT; see next subsection) and the NFFT.
%, and are given in Algorithms \ref{alg:NFSFT} and
%\ref{alg:adjointNFSFT}.
It was introduced in \cite{kupo02} while the adjoint variant was
developed in \cite{KePo06}. Other algorithms, mostly for particular
sampling sets, exist (e.g. \cite{drhe,SuTa01,roty}).

%\begin{algorithm}[tb]
%  \caption[NFSFT]{Nonequispaced fast spherical Fourier transform (NFSFT).}
%  \label{alg:NFSFT}
%  \begin{tabbing}
%    Input: \= \kill
%    Input: \> $M, N \in \mathbb{N}$, $\left(\vartheta_{j},
%      \varphi_{j} \right) \in [0,\pi]\times [0,2\pi)$, $j=0,\ldots,M-1$,
%      $\hat{f}_{k}^n \in \mathbb{C}$, $(k,n) \in J_{N}$.
%  \end{tabbing}
%  \begin{algorithmic}[1]
%    \itemsep=1.1ex
%    \STATE For $n=-N,\hdots,N$ convert $\sum_{k=|n|}^N \hat{f}_k^n
%    P_k^{n}(\cos\vartheta_j)$ into $\sum_{k=-N}^N \hat{g}_k^n \, \mathrm{e}^{k
%    \vartheta_{j}}$ by FPTs.
%    \STATE For $j=0,\hdots,M-1$ compute $f_j=\sum_{n=-N}^N\sum_{k=-N}^N
%    \hat{g}_k^n \, \mathrm{e}^{n\varphi_j+k \vartheta_{j}}$ by an NFFT.\\[2ex]
%  \end{algorithmic}
%  Output: Approximated function values $f_j$, $j=0,\ldots,M-1$.\\
%  Complexity: $\mathcal{O}\left(N^2 \log^2 N + M\right)$.
%\end{algorithm}

%\begin{algorithm}[tb]
%  \caption[Adjoint NFSFT]{Adjoint nonequispaced fast spherical Fourier
%    transform (adjoint NFSFT).}
%  \label{alg:adjointNFSFT}
%  \begin{tabbing}
%    Input: \= \kill
%    Input: \> $M, N \in \mathbb{N}$, $\left(\vartheta_{j},
%      \varphi_{j} \right) \in [0,\pi]\times [0,2\pi)$, $f_j \in\C$,
%      $j=0,\ldots,M-1$.
%  \end{tabbing}
%  \begin{algorithmic}[1]
%    \itemsep=1.1ex
%    \STATE For $k,n = -N,\hdots,N$ compute $\tilde{h}_{k}^n=\sum_{j=0}^{M-1}
%      f_j \mathrm{e}^{n\varphi_j+k \vartheta_{j}}$ by an adjoint NFFT.
%    \STATE
%      For $n=-N,\hdots,N$ convert $\sum_{k=-N}^N \tilde{h}_{k}^n \,
%      \mathrm{e}^{k \vartheta_{j}}$ into $\sum_{k=|n|}^N \hat{h}_k^n
%    P_k^{n}(\cos\vartheta_j)$ by adjoint FPTs.\\[2ex]
%  \end{algorithmic}
%  Output: Approximated coefficients $\hat{h}_{k}^n$, $(k,n) \in J_N$.\\
%  Complexity: $\mathcal{O}\left(N^2 \log^2 N + M\right)$.
%\end{algorithm}

%-------------------------------------------------------------------------------
\subsection{FPT - fast polynomial transform} \label{sect:fpt}
%-------------------------------------------------------------------------------
A \emph{discrete polynomial transform} (DPT) is a generalization of the DFT
from complex exponentials $\mathrm{e}^{\mathrm{i} k x}$ to
systems of algebraic polynomials that satisfy a three-term recurrence;
see \cite[pp. 773]{abst}.
More precisely, let $p_{0},p_{1},\hdots:[-1,1]\rightarrow\R$ be a sequence of
polynomials satisfying the recurrence
\begin{align}
  \label{eq.threeterm}
  p_{k+1}(x) & = (\alpha_{k} x  + \beta_{k})p_{k}(x) + \gamma_{k} p_{k-1}(x) &
  (k \ge 0),
\end{align}
where $p_{-1}(x) := 0$, $p_{0}(x) := 1$, and $\alpha_{k} \neq 0$, $k\ge 0$.
Then every $p_{k}$ is a degree-$k$ polynomial. Typical examples are the
Legendre polynomials $P_k(x)$ or the Chebyshev polynomials of first
kind $T_{k}(x) := \cos(k \arccos x)$. Let a finite linear expansion of the form
\begin{equation*}
  f(x) = \sum_{k=0}^{N} a_{k} p_{k}(x)
\end{equation*}
be given. The DPT in our sense is defined as the
computation of the coefficients $b_{k}$ in the expansion
\[
  f(x) = \sum_{k=0}^N b_{k} T_{k}(x).
\]
Implemented directly, this computation takes $\cO(N^2)$ arithmetic operations.
The \emph{fast polynomial transform (FPT)} is an $\cO(N \log^2 N)$
algorithm to compute the same result. It is implemented in the
NFFT~3 library following the approach in \cite{postta98}. This is based on using
the three-term-recurrence relation \eqref{eq.threeterm} repeatedly, together
with a fast method for polynomial multiplication in the Chebyshev
basis, and a cascade-like summation process. More details can be found in
\cite{drhe,drhero,postta98,po01} and references therein.

%------------------------------------------------------------------------------
\subsection{Inversion and solver module} \label{sect:general_proc_infft}
%------------------------------------------------------------------------------
The NDFT differs from the DFT also in that an inversion formula similar to
\eqref{eq.idft} generally does not exist; see Section \ref{subsect.definition}.
A more general view is at order to define what should be an inverse NDFT/NFFT
transform. This section is devoted to explain our approach for the NFFT. But
everything also applies to the mentioned NFFT variants.

Inversion of the NDFT means the computation of Fourier coefficients
$\hat{f}_{\mathbf{k}}$, $\mathbf{k} \in I_{\mathbf{N}}$, from given function
samples $y_{j}$ at nodes $\mathbf{x}_j$, $j = 0,\ldots,M-1$.
In matrix-vector notation, this is equivalent to solving the linear system
\begin{equation}\label{eq:infft}
  \mathbf{A} \mathbf{\hat{f}} \approx \mathbf{y}.
\end{equation}
%for the vector $\mathbf{\hat{f}}\in\C^{|I_{\mathbf{N}}|}$.
This can be either and overdetermined system if $|I_{\mathbf{N}}| \le
M$ (this includes the quadratic case) or an underdetermined system
if $|I_{\mathbf{N}}| > M$.
Generally, this compels us to look for a pseudo-inverse solution
$\mathbf{\hat{f}}^{\dagger}$ (see e.g. \cite[p. 15]{Bj96}).
For this, we also require that the nonequispaced Fourier matrix $\mathbf{A}$
have full rank. Eigenvalue estimates in \cite{FeGrSt95,BaGr03,KuPo04} indeed assert that this
condition is satisfied if the sampling set $\cal X$ is uniformly dense or uniformly
separated with respect to the inverse bandwidth $\mathbf{N}^{-1}$.

For the overdetermined case, we now consider a weighted least squares problem,
while for the consistent underdetermined case, we look for a solution of an
interpolation problem. Both can be solved by iterative algorithms using
NFFT and adjoint NFFT to realize fast matrix-vector multiplications
involving $\mathbf{A}$ and $\mathbf{A}^{\adj}$, respectively.

%-------------------------------------------------------------------------------
\subsubsection*{Weighted least squares problem}
%-------------------------------------------------------------------------------
If $|I_{\mathbf{N}}|\le M$, the linear system \eqref{eq:infft} is overdetermined
which typically implies that the given data $y_j\in \C,\;j=0,\hdots,M-1,$ can
only be approximated up to a residual $\mathbf{r}:=\mathbf{y}-\mathbf{A}
\mathbf{\hat{f}}$. Therefore, we consider the \emph{weighted least squares problem}
\begin{equation*}
\sum_{j=0}^{M-1} w_j |y_j-f(x_j)|^2
  \stackrel{\mathbf{\hat{f}}}{\rightarrow} \min
\end{equation*}
with weights $w_j> 0$, $\mathbf{W}:=\diag(w_j)_{j=0,\hdots,M-1}$. This weights
might be used to compensate for clusters in the sampling set ${\cal X}$.
The weighted least squares problem is equivalent to solving
the \emph{weighted normal equation of first kind},
$
\mathbf{A}^{\adj} \mathbf{W} \mathbf{A} \mathbf{\hat{f}} = \mathbf{A}^{\adj}
\mathbf{W} \mathbf{y}.
$
This can be done, e.g., by using the Landweber (or Richardson) iteration,
the steepest descent method, or the conjugate gradient method for least
squares problems. The latter method is given in pseudo-code in Algorithm
\ref{algo:wCGNR(E)} (left). The NFFT~3 library nevertheless implements
all three algorithms.

\setlength{\columnseprule}{0.5pt}

\begin{algorithm}[tp]
  \hrule
  \footnotesize
  \rule[-1.5ex]{0ex}{5ex}Input: $\mathbf{y} \in \mathbb{C}^M,\, {\mathbf{\hat{f}}}_0 \in
    \mathbb{C}^{|I_{\mathbf{N}}|}$.
    \begin{multicols}{2}
      \begin{algorithmic}%[1]
      \topsep=1.1ex
      \itemsep=1.1ex
      \STATE $\mathbf{r}_0=\mathbf{y} - \mathbf{A} \mathbf{\hat{f}}_0$
      \STATE $\mathbf{\hat{z}}_0=\mathbf{A}^{\adj} \mathbf{W} \mathbf{r}_0$
      \STATE $\mathbf{\hat{p}}_0=\mathbf{\hat{z}}_0$
      \FOR{$l=0,\hdots$}
        \topsep=1.1ex
        \itemsep=1.1ex
        \STATE $\mathbf{v}_l=\mathbf{A} \mathbf{\hat{W}} \mathbf{\hat{p}}_l$
        \STATE $\alpha_l=\frac{\mathbf{\hat{z}}_l^{\adj} \mathbf{\hat{W}}
          \mathbf{\hat{z}}_l} {\mathbf{v}_l^{\adj} \mathbf{W} \mathbf{v}_l}$
        \STATE $\mathbf{\hat{f}}_{l+1}=\mathbf{\hat{f}}_{l} +
          \alpha_l \mathbf{\hat{W}} \mathbf{\hat{p}}_l$
        \STATE $\mathbf{r}_{l+1}=\mathbf{r}_{l}- \alpha_l \mathbf{v}_l$
        \STATE $\mathbf{\hat{z}}_{l+1}=\mathbf{A}^{\adj} \mathbf{W}
          \mathbf{r}_{l+1}$
        \STATE $\beta_l=\frac{\mathbf{\hat{z}}_{l+1}^{\adj}
          \mathbf{\hat{W}}\mathbf{\hat{z}}_{l+1}}
          {\mathbf{\hat{z}}_l^{\adj} \mathbf{\hat{W}}\mathbf{\hat{z}}_l}$
        \STATE $\mathbf{\hat{p}}_{l+1}=\beta_l \mathbf{\hat{p}}_l +
          \mathbf{\hat{z}}_{l+1}$
      \ENDFOR
    \end{algorithmic}
\columnbreak
      \begin{algorithmic}
      \topsep=1.1ex
      \itemsep=1.1ex
    \STATE $\mathbf{r}_0=\mathbf{y} - \mathbf{A} \mathbf{\hat{f}}_0$
    \STATE $\mathbf{\hat{p}}_0=\mathbf{A}^{\adj} \mathbf{W} \mathbf{r}_0$
    \FOR{$l=0,\hdots$}
    \topsep=1.1ex
    \itemsep=1.1ex
    \STATE $\alpha_l=\frac{\mathbf{r}_l^{\adj} \mathbf{W} \mathbf{r}_l}
           {\mathbf{\hat{p}}_l^{\adj} \mathbf{\hat{W}} \mathbf{\hat{p}}_l}$
           \STATE $\mathbf{\hat{f}}_{l+1}=\mathbf{\hat{f}}_{l}
           +\alpha_l \mathbf{\hat{W}} \mathbf{\hat{p}}_l$

           \STATE $\mathbf{r}_{l+1}=\mathbf{r}_{l}- \alpha_l \mathbf{A}
\mathbf{\hat{W}} \mathbf{\hat{p}}_l$

           \STATE $\beta_l=\frac{\mathbf{r}_{l+1}^{\adj} \mathbf{W}
\mathbf{r}_{l+1}}
                  {\mathbf{r}_l^{\adj} \mathbf{W} \mathbf{r}_l}$

                  \STATE $\mathbf{\hat{p}}_{l+1}=\beta_l \mathbf{\hat{p}}_l +
\mathbf{A}^{\adj} \mathbf{W} \mathbf{r}_{l+1}$
                  \ENDFOR
    \end{algorithmic}
\end{multicols}
    Output: The $l$th approximation to the solution vector
    $\mathbf{\hat{f}}_{l}\in\mathbb{C}^{|I_{\mathbf{N}}|}$.\rule[-1.5ex]{0ex}{3ex}
  \hrule
  \caption{Variants of the conjugate gradient method for normal equations of first
  (CGNR = residual minimization; left) and second kind (CGNE = error minimization;
  right).}
  \label{algo:wCGNR(E)}
\end{algorithm}

%-------------------------------------------------------------------------------
\subsubsection*{Interpolation problem}
%-------------------------------------------------------------------------------
If $|I_{\mathbf{N}}|>M$, and if the linear system \eqref{eq:infft} is
consistent, the data $y_j\in \C,\;j=0,\hdots,M-1,$ can be interpolated exactly.
But since there exist multiple solutions, we consider the
\emph{constrained minimization problem}
\begin{equation*}
  \label{eq:ip}
  \sum_{\mathbf{k}\in I_{\mathbf{N}}} \frac{|\hat f_{\mathbf{k}}|^2}{\hat
  w_{\mathbf{k}}} \stackrel{\mathbf{\hat{f}}}{\rightarrow} \min \quad
  \text{subject to} \quad \mathbf{A} \mathbf{\hat{f}} = \mathbf{y},
\end{equation*}
which incorporates ``damping factors'' $\hat w_{\mathbf{k}}> 0$,
$\mathbf{\hat{W}}:= \diag(\hat w_{\mathbf{k}})_{\mathbf{k}\in
I_{\mathbf{N}}}$. A smooth solution, i.e., a solution with rapid decay of
Fourier coefficients $\hat{f}_{\mathbf{k}}$, is favored if the damping
factors $\hat{w}_{\mathbf{k}}$ are decreasing themselves.
The interpolation problem is equivalent to the
\emph{damped normal equation of second kind}
$
  \mathbf{A} \mathbf{\hat{W}} \mathbf{A}^{\adj} \mathbf{\tilde{f}}
  = \mathbf{y}
$,
$\mathbf{\hat{f}} =
  \mathbf{\hat{W}} \mathbf{A}^{\adj} \mathbf{\tilde{f}}
$.
The conjugate gradient method applied to this linear system is given in
Algorithm \ref{algo:wCGNR(E)} (right). Of course, the NFFT~3 library also
implements this scheme.

%==============================================================================
\section{Summary and performance analysis}\label{sect:sum}
%==============================================================================

This section is to provide a performance analysis of the NFFT algorithm.
The asymptotic bounds on the number of floating point operations for the
nonequispaced fast Fourier transform (NFFT) and its variant on the sphere
(NFSFT) are compared to those for the discrete/fast Fourier transform in
Table \ref{table:bigO}.
\begin{acmtable}{0.85\textwidth}
  \centering%
  \vspace{1.3ex}%
  \begin{tabular}{r|cc|cc}
    & \multicolumn{2}{c|}{\rule[-1.5ex]{0ex}{4.25ex}$M={\mathcal{O}}\left(N^d\right)$, $\cal X$ grid}
    & \multicolumn{2}{c}{$M$, $\cal X$ arbitrary}
    \\
    \hline
    & DFT & FFT & NDFT & NFFT
    \\
    \hline
    \rule[-1.5ex]{0ex}{4.25ex}%
    Torus $\T^d$ &
    ${\mathcal{O}}\left(N^{d+1}\right)$ &
    ${\mathcal{O}}\left(N^d \log N\right)$ &
    ${\mathcal{O}}\left(M N^d\right)$ &
    ${\mathcal{O}}\left(N^d \log N +|\log \varepsilon|^d M\right)$\\
    \rule[-1.5ex]{0ex}{4.25ex}%
    Sphere $\S^2$ &
    ${\mathcal{O}}\left(N^3\right)$ &
    ${\mathcal{O}}\left(N^2 \log^2 N\right)$ &
    ${\mathcal{O}}\left(M N^2\right)$ &
    ${\mathcal{O}}\left(N^2 \log^2 N +|\log \varepsilon|^2 M\right)$
  \end{tabular}
  \vspace{0.5ex}%
  \caption{Asymptotic number of floating point operations with respect to the
  polynomial degree $N$ (so that $\mathbf{N}=(N,\hdots,N)^{\top}$), the number of nodes
  $M$, spatial dimension $d$, and target accuracy $\varepsilon$ for different transforms.}
  \label{table:bigO}
\end{acmtable}
The rest of the section is to provide results of actual performance
measurements for our implementation of the NFFT in spatial
dimensions $d=1,2,3$. Computation time measurements and a numerical
confirmation of the error estimates for the window functions
(\ref{gaussian}-\ref{kaiser}) with respect to the so-called cut-off
parameter are given.

For our tests, we used a computer system with an AMD
$\text{Athlon}^{\text{\texttrademark}}$ XP 2700+ processor, 1GB memory,
SuSE Linux with kernel 2.4.20-4GB-athlon, and the GCC 3.3 compiler.
All computations were carried out in double precision using the compiler
option \verb"-O3". For all tests with random input, the nodes
$\mathbf{x}_j$ and the Fourier coefficients $\hat f_{\mathbf{k}}$ were
drawn from a uniform pseudo-random distribution in the respective domains,
i.e., $\mathbf{x}_j \in \mathbb{T}^d$ and
$\hat f_{\mathbf{k}} \in [0,1] \times [0,1]{\rm i}$. The code for these
and more examples and applications can be found in the directories
\verb"examples/" and \verb"applications/" of the NFFT~3 package.

%------------------------------------------------------------------------------
\subsection{Computation time with respect to problem size}
%------------------------------------------------------------------------------
We compared computation times for ordinary FFTs (using FFTW
3.0.0 \cite{FFTW05,fftw} in double precision with the flag
\texttt{FFTW\_MEASURE}),
direct nonequispaced discrete Fourier transforms (NDFTs), and nonequispaced fast
Fourier transforms (NFFTs) for different problem sizes
$|I_{\mathbf{N}}|$ in dimensions $d=1,2,3$.
The multibandwidth $\mathbf{N}$ was chosen to
$\mathbf{N}=(N,\hdots,N)^{\top},\;N\in 2\N$.
The nodes for FFTs are fixed and restricted to the lattice $\mathbf{N}^{-1}\odot
I_{\mathbf{N}}$.
For NDFTs and NFFTs, we chose $M = N^d$ random nodes and used an oversampling
factor $\sigma=2$, cut-off parameters $m = 2$ and $m = 4$, and the
Kaiser-Bessel window function (precomputation flags \verb"PRE_PSI" and
\verb"PRE_PHI_HUT").
With these options, the accuracy measured by
\begin{equation}\label{Einf}
E_{\infty}:={\max_{0\le
    j<M}|f_j-s_j|} / {\sum_{\mathbf{k} \in I_{\mathbf{N}}} |\hat
  f_{\mathbf{k}}|}
\end{equation}
is approximately $E_{\infty}\approx 10^{-4}$ ($m=2$) and
$E_{\infty}\approx 10^{-8}$ ($m=4$) for $d=1,2,3$.
Here, $s_j$ denotes the NFFT-approximation to $f_j$; see the Appendix for
details. The corresponding example program can be found under
\verb"examples/nfft/nfft_times".
As a byproduct of the measurements shown in Table \ref{tab:TimeSpace_1}, we
also obtain a performance index measured in million floating point operations
per second (MFLOPS). We assume that the number of floating point operations are
\begin{equation*}
  5 |I_{\mathbf{N}}|\log_2|I_{\mathbf{N}}|, \quad \text{and} \quad
  |I_{\mathbf{N}}| + 5\cdot 2^d |I_{\mathbf{N}}|\log(2^d|I_{\mathbf{N}}|) +
  4(2m+1)^d |I_{\mathbf{N}}|
\end{equation*}
for the FFT and the NFFT ($M = |I_{\mathbf{N}}|$ nodes, oversampling factor
$\sigma=2$, and precomputation flags \verb"PRE_PSI" and \verb"PRE_PHI_HUT"),
respectively. More details are given in the Appendix.

\begin{acmtable}{0.85\textwidth}
  \centering%
  \vspace{1.3ex}%
    \begin{tabular}{r|r|r|r|r r|r r}
    & & & & \multicolumn{2}{|c|}{$m=4$, $E_{\infty} \approx 10^{-8}$} & \multicolumn{2}{|c}{$m=2$, $E_{\infty}\approx 10^{-4}$}
    \\
    \hline
    & \rule[-1.8ex]{0ex}{4.5ex}
    $l_N$ & FFT & NDFT & NFFT & \parbox[c]{1cm}{\scriptsize\centering penalty\\factor}
    & NFFT & \parbox[c]{1cm}{\scriptsize\centering penalty\\factor}
    \\
    \hline
    \multirow{20}{*}{$d=1$}
& \verb"3"&       \verb"1.3e-07" &        \verb"8.6e-06" &        \verb"1.3e-06" & \verb"(10.0)"  &        \verb"1.1e-06" & \verb"(8.6)"\\
& \verb"4"&       \verb"2.1e-07" &        \verb"3.5e-05" &        \verb"2.7e-06" & \verb"(12.9)"  &        \verb"2.2e-06" & \verb"(10.5)"\\
& \verb"5"&       \verb"4.2e-07" &        \verb"1.4e-04" &        \verb"5.2e-06" & \verb"(12.3)"  &        \verb"4.2e-06" & \verb"(11.4)"\\
& \verb"6"&       \verb"9.2e-07" &        \verb"5.7e-04" &        \verb"1.1e-05" & \verb"(11.5)"  &        \verb"9.1e-06" & \verb"(9.8)"\\
& \verb"7"&       \verb"2.2e-06" &        \verb"2.3e-03" &        \verb"2.1e-05" & \verb"(9.5)"   &        \verb"1.9e-05" & \verb"(9.0)"\\
& \verb"8"&       \verb"5.3e-06" &        \verb"9.2e-03" &        \verb"4.4e-05" & \verb"(8.3)"   &        \verb"3.9e-05" & \verb"(7.7)"\\
& \verb"9"&       \verb"1.1e-05" &        \verb"3.7e-02" &        \verb"9.8e-05" & \verb"(8.8)"   &        \verb"8.3e-05" & \verb"(7.4)"\\
& \verb"10"&      \verb"2.5e-05" &        \verb"1.5e-01" &        \verb"2.0e-04" & \verb"(8.1)"   &        \verb"1.8e-04" & \verb"(7.1)"\\
& \verb"11"&      \verb"6.0e-05" &        \verb"6.0e-01" &        \verb"6.5e-04" & \verb"(10.8)"  &        \verb"5.4e-04" & \verb"(8.9)"\\
& \verb"12"&      \verb"1.5e-04" &        \verb"2.4e+00" &        \verb"1.5e-03" & \verb"(9.8)"   &        \verb"1.4e-03" & \verb"(8.9)"\\
& \verb"13"&      \verb"5.5e-04" &        \verb"9.6e+00" &        \verb"3.5e-03" & \verb"(6.4)"   &        \verb"3.3e-03" & \verb"(5.6)"\\
& \verb"14"&      \verb"1.7e-03" &        \verb"4.0e+01" &        \verb"7.8e-03" & \verb"(4.6)"   &        \verb"7.1e-03" & \verb"(4.3)"\\
& \verb"15"&      \verb"4.0e-03" &        \verb"1.6e+02" &        \verb"1.6e-02" & \verb"(4.1)"   &        \verb"1.5e-02" & \verb"(3.7)"\\
& \verb"16"&      \verb"8.4e-03" &        \verb"*"        &       \verb"3.4e-02" & \verb"(4.1)"   &       \verb"3.2e-02" & \verb"(3.8)"\\
& \verb"17"&      \verb"2.0e-02" &        \verb"*"        &       \verb"7.5e-02" & \verb"(3.8)"   &       \verb"7.2e-02" & \verb"(3.6)"\\
& \verb"18"&      \verb"4.6e-02" &        \verb"*"        &       \verb"1.6e-01" & \verb"(3.4)"   &       \verb"1.5e-01" & \verb"(3.2)"\\
& \verb"19"&      \verb"9.5e-02" &        \verb"*"        &       \verb"3.2e-01" & \verb"(3.4)"   &       \verb"3.0e-01" & \verb"(3.1)"\\
& \verb"20"&      \verb"2.1e-01" &        \verb"*"        &       \verb"6.9e-01" & \verb"(3.3)"   &       \verb"6.4e-01" & \verb"(3.3)"\\
& \verb"21"&      \verb"4.3e-01" &        \verb"*"        &       \verb"1.5e+00" & \verb"(3.5)"   &       \verb"1.4e+00" & \verb"(3.2)"\\
& \verb"22"&      \verb"1.0e+00" &        \verb"*"        &       \verb"3.2e+00" & \verb"(3.1)"   &       \verb"3.0e+00" & \verb"(2.9)"\\
      \hline
%      \hline \multicolumn{3}{|l|}{$d=2$} & \multicolumn{2}{|c|}{$m=4$, $E_{\infty}\approx
%  10^{-8}$} & \multicolumn{2}{|c|}{$m=2$, $E_{\infty}\approx 10^{-4}$} \\
% \hline $l_N$ & FFT & NDFT & NFFT & (/FFT) & NFFT & (/FFT)\\
    \multirow{9}{*}{$d=2$}
& \verb"6"&       \verb"9.0e-07" &        \verb"6.0e-04" &        \verb"6.3e-05"  & \verb"(70.4)" &        \verb"3.2e-05" & \verb"(33.2)"\\
& \verb"8"&       \verb"4.4e-06" &        \verb"9.5e-03" &        \verb"2.5e-04"  & \verb"(58.0)" &        \verb"1.3e-04" & \verb"(29.0)"\\
& \verb"10"&      \verb"2.2e-05" &        \verb"1.5e-01" &        \verb"1.2e-03"  & \verb"(55.0)" &        \verb"6.0e-04" & \verb"(27.4)"\\
& \verb"12"&      \verb"1.2e-04" &        \verb"2.4e+00" &        \verb"6.4e-03"  & \verb"(52.5)" &        \verb"4.0e-03" & \verb"(32.2)"\\
& \verb"14"&      \verb"1.7e-03" &        \verb"4.0e+01" &        \verb"4.0e-02"  & \verb"(23.4)" &        \verb"3.1e-02" & \verb"(18.1)"\\
& \verb"16"&      \verb"2.2e-02" &        \verb"*"        &       \verb"1.7e-01"  & \verb"(7.8)"  &       \verb"1.3e-01" & \verb"(6.1)"\\
& \verb"18"&      \verb"8.7e-02" &        \verb"*"        &       \verb"6.7e-01"  & \verb"(7.7)"  &       \verb"5.0e-01" & \verb"(5.8)"\\
& \verb"20"&      \verb"3.3e-01" &        \verb"*"        &       \verb"3.0e+00"  & \verb"(9.1)"  &       \verb"2.1e+00" & \verb"(6.4)"\\
& \verb"22"&      \verb"1.4e+00" &        \verb"*"        &       \verb"1.4e+01"  & \verb"(10.1)" &       \verb"1.1e+01" & \verb"(7.9)"\\
      \hline
 %     \hline \multicolumn{3}{|l|}{$d=3$} & \multicolumn{2}{|c|}{$m=4$, $E_{\infty}\approx
 % 10^{-8}$} & \multicolumn{2}{|c|}{$m=2$, $E_{\infty}\approx 10^{-4}$} \\
 % \hline $l_N$ & FFT & NDFT & NFFT & (/FFT) & NFFT & (/FFT)\\
    \multirow{5}{*}{$d=3$}
& \verb"9"&       \verb"1.0e-05" &        \verb"3.8e-02" &        \verb"4.4e-03" & \verb"(423.4)"  &        \verb"1.1e-03" & \verb"(109.5)"\\
& \verb"12"&      \verb"1.1e-04" &        \verb"2.4e+00" &        \verb"4.2e-02" & \verb"(369.1)"  &        \verb"1.4e-02" & \verb"(115.8)"\\
& \verb"15"&      \verb"3.5e-03" &        \verb"1.6e+02" &        \verb"3.9e-01" & \verb"(110.8)"  &        \verb"1.5e-01" & \verb"(41.8)"\\
& \verb"18"&      \verb"3.9e-02" &        \verb"*"        &       \verb"3.9e+00" & \verb"(99.3)"   &       \verb"1.7e+00" & \verb"(42.4)"\\
& \verb"21"&      \verb"9.4e-01" &        \verb"*"        &       \verb"8.4e+01" & \verb"(89.0)"   &       \verb"1.8e+01" & \verb"(19.6)"\\
    \end{tabular}
  \vspace{0.5ex}%
  \caption{Computation times in seconds with respect to the problem size $l_N=\log_2|I_{\mathbf{N}}|$.
    Numbers in parentheses give the penalty factor, i.e., the quotient of computation times
    for NFFT and FFT.
    Times for small $l_N$ have been averaged over a series of transforms. An
    asterisk (\texttt{*}) indicates values not calculated due to excessive computation
    time.}
  \label{tab:TimeSpace_1}
\end{acmtable}

Some conclusions can be drawn from Table \ref{tab:TimeSpace_1}:
Both, FFT and NFFT, require a computation time of order
$\cO(|I_{\mathbf{N}}|\log|I_{\mathbf{N}}|)$ as expected.
Doubling $|I_{\mathbf{N}}|$ indeed results in only slightly more
than twice the computation time. The NDFT also scales as the expected
$\cO(|I_{\mathbf{N}}|^2)$.
For FFT and NDFT, the constant in the $\cO$-notation is independent of the
space dimension $d$, whereas the computation time of the NFFT increases
considerably with $d$ fixing $|I_{\mathbf{N}}|$.
Moreover, Figure \ref{fig:flops} shows a relative performance gain of the FFT
for small and moderate sized problems.
In contrast, the overhead contained in the approximation scheme of the NFFT algorithm
seems to attenuate this effect for the NFFT.

\goodbreak

%-------------------------------------------------------------------------------
\subsection{Accuracy}
%-------------------------------------------------------------------------------
The actual choice of the window function and the cut-off parameter $m$ (cf.
Appendix) have direct influence on the accuracy achieved by the NFFT algorithm.
To confirm this, we compare values $f_j$ in \eqref{eq:nfft} computed by the
NDFT with approximations $s_{j}$ computed using an NFFT with different
choices of the window function and the cut-off parameter $m$.
The error is measure is again \eqref{Einf}.
Figure \ref{fig:accuracy} shows results for $1$, $2$, and $3$-dimensional
transforms with equal total numbers of Fourier coefficients $|I_{\mathbf{N}}|$.

\begin{figure}[tp]
  \centering
  \includegraphics[width=0.9\textwidth]{images/nfft_times_flops}
  \caption{Performance measured in MFLOPS achieved by the NFFT~3 library
    with respect to the logarithm of the total problem size $l_N=\log_2|I_{\mathbf{N}}|$.
    Shown are results for $d=1$ (solid), $d=2$ (dotted), and $d=3$ (dashed).
    The tested algorithms are the FFT ($\circ$), the NFFT with $m=4$ ($\times$),
    and the NFFT with $m=2$ ($+$).}
  \label{fig:flops}
\end{figure}

\begin{figure}[tp]
  \subfigure[$d = 1$, $N =
2^{12}$.]{\includegraphics[width=0.3\textwidth]{images/accuracy4}}
  \hfill
  \subfigure[$d = 2$, $N =
2^{6}$.]{\includegraphics[width=0.3\textwidth]{images/accuracy5}}
  \hfill
  \subfigure[$d = 3$, $N =
2^{4}$.]{\includegraphics[width=0.3\textwidth]{images/accuracy6}}
  \caption{Error $E_{\infty}$ for increasing cut-off parameter $m=0,\hdots,14$ and
    $d=1,2,3$. In each case, the degree $N$ was chosen to be equal
    along each dimension such that $|I_{\mathbf{N}}| = 2^{12}$. We fixed
    $\sigma = 2$ (cf. Appendix) and $M = 10000$. Shown are results for
    the Kaiser Bessel ($\circ$), the Sinc ($\times$), the B-spline ($+$), and
    the Gaussian window function ($\triangle$). \label{fig:accuracy}}
 \end{figure}

%==============================================================================
\section{Applications}\label{sect:exap}
%==============================================================================

The rest of this tutorial is devoted to illustrating the usage of the NFFT~3
library for more advanced applications.
We start with two small problems from scattered data approximation which
use the solver module and the NFFT and NFSFT routines.
Then we briefly sketch some more important applications - fast
summation schemes and medical imaging.
These and more examples and applications can be found in the directories
\verb"examples/" and \verb"applications/" of the package.


%-------------------------------------------------------------------------------
\subsection{Scattered data approximation - computing an inverse
transform}\label{sect:rec_example}
%-------------------------------------------------------------------------------

Our first example demonstrates the reconstruction of the surface of a glacier
from data available on level sets. %\cite{Franke.Daten}
This is done by solving a damped optimal interpolation problem.
Figure \ref{fig:glacier} shows the result using
the $M = 8345$ available data points and bandwidths $N_0 = N_1 = 256$ for
the two dimensions after $40$ iterations of the CGNE method.
The corresponding example program \verb+glacier+ can be found in
\verb+examples/solver/+. There is also a simple test program called
\verb"simple_test" to illustrate basic principles.
\begin{figure}[tp]
  \centering
  \subfigure[Surface plot.]{\includegraphics[width=0.45\textwidth]{images/glacier1}}
  \hfill
  \subfigure[Contour plot.]{\includegraphics[width=0.45\textwidth]{images/glacier2}}
 \caption{Reconstruction of the glacier from data on level sets at $M=8345$%, cf.
%\cite[\texttt{vol87.dat}]{Franke.Daten},
with $N_0=N_1=256$ after $40$ iterations. \label{fig:glacier}}
\end{figure}

The second example uses a map from the NASA AMSU mission \cite{NASA_AIRS}
containing global atmospheric temperature data of the earth from 5 November 2006.
The map contains empty strips not covered by the trajectory of the satellite
used to take the measurements. To fill these gaps, we solve a weighted linear least
squares problem to obtain a spherical harmonics
approximation of degree $N = 128$ to the given data (see \cite{KeKuPo06b}).
We then use the values of the approximation to fill the gaps. Figure
\ref{fig:temperature} shows the original data and the map with the gaps filled.

\begin{figure}[tp]
  \centering
  \subfigure[Given data.]
    {\includegraphics[width=0.45\textwidth]{images/temperaturemap_gaps} } \hfill
  \subfigure[Approximation.]
    {\includegraphics[width=0.45\textwidth]{images/temperaturemap_filled}}
  \caption{Original atmospheric temperature data in Kelvin (a) and the same
    data with values filled in from the least squares approximation at degree
    $N = 128$.\label{fig:temperature}}
\end{figure}

The following subsections are devoted to the demonstration of further important
applications of NFFT~3 routines. The source code is available in
the directory \verb"applications/" in the package.

%-------------------------------------------------------------------------------
\subsection{Fast summation}
%-------------------------------------------------------------------------------
An important task in many applications is to efficiently evaluate a linear
combinations of radial functions, i.e.,
\begin{align*}
  g\left(\mathbf{y}_j\right) & := \sum_{l=0}^{L-1} \alpha_l
K\left(\left\|\mathbf{y}_j-
  \mathbf{x}_l\right\|_2\right) & (j=0,\ldots,M-1)
\end{align*}
with source and target nodes $\mathbf{x}_l, \mathbf{y}_j \in \R^d$,
and real coefficients $\alpha_l \in \R$.

For kernels $K$ with a certain smoothness, e.g. the Gaussian
$K(x)=\e^{-x^2/c^2}$, the multiquadric $K(x)=\sqrt{x^2+c^2}$, or the inverse
multiquadric $K(x)=1/\sqrt{x^2+c^2}$, all with a parameter $c>0$,
our algorithm requires $\cO(L+M)$ arithmetic operations instead of
$\mathcal{O}(LM)$ via direct calculation.
For singular kernels $K$, e.g.
\begin{equation*}
  \frac{1}{x^2},\; \frac{1}{|x|}, \; \log |x|, \;
  x^2 \log |x|, \; \frac{1}{x}, \;
  \frac{\sin(c x)}{x},\; \frac{\cos(c x)}{x}, \; \cot(cx)
\end{equation*}
it is imperative to apply an additional regularization procedure. The cost then
becomes $\cO(L \log L +M)$ or $\cO(M\log M + L)$ provided that
either the source nodes $\mathbf{x}_l$ or the target nodes $\mathbf{y}_j$ are
uniformly distributed to a reasonable degree.
This fast method, which was proposed in \cite{post02,postni04,FeSt03},
generalizes the diagonalization of convolution matrices by Fourier matrices to
settings with arbitrary nodes.
Also, our method yields nearly the same cost as the FMM \cite{beatgreen97}
while allowing for an easy exchange of the kernel $K$.
A recent application in particle simulation is given in \cite{popori05}.
The directory \texttt{applications/fastsum} contains C and {\sc MatLab} programs
showing how to use the fast summation method.

%-------------------------------------------------------------------------------
%\subsection{Fast Gauss transform}
%-------------------------------------------------------------------------------
As a special case, let us consider the fast Gauss transform. This is the calculation of sums of the form
\begin{align*}
  g\left(y_j\right) & :=\sum_{l=0}^{L-1} \alpha_l {\rm e}^{-\sigma|y_j-x_l|^2} &
  (j=0,\ldots,M-1)
\end{align*}
for source and target nodes $x_l,y_{j}\in[-\frac{1}{4},\frac{1}{4}]$,
complex coefficients $\alpha_l \in \C$, and a complex parameter
$\sigma =a +{\rm i} b$, $a>0$, $b\in\R$ with $\mathcal{O}(L + M)$ operations.
For details see \cite{AnBe04,KuPoSt06}.
Numerical examples found in \cite{KuPoSt06} can be reproduced with the programs
found in the directory \texttt{applications/fastgauss}.

%-------------------------------------------------------------------------------
%\subsection{Summation of zonal functions on the sphere}
%-------------------------------------------------------------------------------
On the sphere $\mathbb{S}^2$, the equivalent of radial functions are \emph{zonal
functions}. Fast summation of zonal functions on the sphere is the computation
of
sums of the form
\begin{equation*}
  g\left(\mathbf{\xi}_j\right) := \sum_{l = 0}^{L-1} \alpha_l K\left(\mathbf{\eta}_l
  \cdot \mathbf{\xi}_j\right) \qquad (j=0,\ldots,M-1)
\end{equation*}
for source and target nodes $\mathbf{\eta}_l,\mathbf{\xi}_j \in \S^2$, and real
coefficients $\alpha_l \in \R$ with $\cO(L+M)$ operations.
Our algorithm is based on the nonequispaced fast spherical Fourier
transform. It and can be easily adapted to different kernels $K$, e.g. the
Poisson kernel, smooth locally supported kernels, or the spherical Gaussian
kernel. For details see \cite{KeKuPo06}. Numerical examples
are found under \verb+applications/fastsumS2/+.


%-------------------------------------------------------------------------------
\subsection{Applications in medical imaging}
%-------------------------------------------------------------------------------
%-------------------------------------------------------------------------------
%\subsection{Iterative reconstruction in magnetic resonance imaging}
%-------------------------------------------------------------------------------
In magnetic resonance imaging (MRI), raw data is measured in \emph{k-space}, the
domain of spatial frequencies. Methods that use non-Cartesian sampling grids
in k-space, e.g. a spiral, have received increasing attention recently.
Reconstruction is usually done by resampling the data onto a Cartesian
grid to use the standard FFT. This is often called \emph{gridding}.
A so-called inverse model is based on an implicit discretization and both
methods can be implemented efficiently using the NFFT and the inverse NFFT,
respectively (see \cite{KnKuPo} and references therein).
Furthermore, approaches to field inhomogeneity correction have been proposed
in \cite{SuNoFe03,EgKnPo06}.
Numerical tests are found in \verb+applications/mri/+.

%-------------------------------------------------------------------------------
%\subsection{Radon transform, computer tomography, and ridgelet transform}
%-------------------------------------------------------------------------------
In computerized tomography (CT) an $N \times N$ (medical) image is to be
reconstructed from its Radon transform.
The standard reconstruction algorithm, the filtered backprojection, yields
reasonably good images at the expense of $\cO( N^3 )$ arithmetic operations.
Fourier reconstruction methods reduce the number of arithmetic operations to
${\cal O}( N^2 \log N )$.
Unfortunately, straightforward Fourier reconstruction algorithms suffer
from unacceptable artifacts making them useless in practice. Better quality of
the reconstructed images can be achieved by our algorithm which is based on
NFFTs. For details see \cite{PoSt01IMA,PoSt00SPIE,PoSt02JFAA} and references
therein, for examples we refer to \verb+applications/radon/+.

Another application of the discrete Radon transform is the \emph{discrete
ridgelet transform}, see e.g. \cite{CaDeDoYi} and the references therein.
A simple test program for denoising images by applying hard thresholding to
ridgelet coefficients \cite{FeMa05} can be found in \texttt{applications/radon}.
It uses the NFFT-based discrete Radon transform and the translation-invariant
discrete Wavelet transform from the \textsc{MatLab} toolbox WaveLab850
\cite{WaveLab}.

%-------------------------------------------------------------------------------
%\subsection{Computation of the polar FFT}
%-------------------------------------------------------------------------------
Our last example to be mentioned here is the \emph{polar FFT}, which is a special case
of the NFFT for a particular set of nodes.
Of course, the polar as well as a so-called pseudo-polar FFT can be calculated
efficiently accurately with the NFFT.
Furthermore, the reconstruction of a $2d$ signal from its Fourier transform
samples on a (pseudo-)polar grid by means of the inverse nonequispaced FFT is
possible under certain density assumptions.
For details see \cite{AvCoDoElIs,FeKuPo06,BeKuMo07}. Numerical tests can be
found in the directory \verb+applications/polarFFT/+.

\appendix
\section*{APPENDIX}

%-------------------------------------------------------------------------------
\section{Installation} \label{sect:inst}
%-------------------------------------------------------------------------------

The installation of the NFFT~3 software library follows the usual GNU
principles.
Download the most recent version from
\verb"http://www.tu-chemnitz.de/~potts/nfft".
We assume that you use a \verb"bash" compatible shell. Uncompress the archive
and
change to the newly created directory:
\begin{verbatim}
  tar xfvz nfft3.x.x.tar.gz
  cd nfft3.x.x
\end{verbatim}
%Before you compile and install the library, you might want to adjust the
%environment variables \verb"CPPFLAGS" and \verb"LDFLAGS" so that NFFT~3 can find
%the FFTW library if it has been installed in a non-standard location and/or if
%the NFFT~3 library shall be compiled with specific compiler flags. Leave this
%step
%out, if you are unsure what to do. Return to this step if the configure script,
%the compiler or the linker complain about missing files or libraries. Here is
%an example that adds the compiler flag \verb"-O3" for optimisations and
%adds a path to include files and libraries in a non-standard location:
%\begin{itemize}
%  \itemsep-1ex
%  \item
%    \verb+export CPPFLAGS="$(CPPFLAGS) -O3 -I/path/to/include/files"+
%  \item \verb+export LDFLAGS="$(LDFLAGS) -L/path/to/libraries"+
%\end{itemize}
Figure \ref{nfft:fig:svn} gives an overview over the directory structure of
the NFFT~3 package.
\begin{figure}[tp]%tb!]
  \centering
  \input{tree}
  \caption{Directory structure of the NFFT~3 package.
    \label{nfft:fig:svn}}
\end{figure}
The package is now ready to be configured. To this end, run the configure script
\verb"configure" which will adapt the build process to your computer system.
After that, compile the library by invoking the \verb"make" utility:
\begin{verbatim}
  ./configure
  make
\end{verbatim}
When the build process has finished, you can run \verb"make install" to install the
library permanently on your system.
%If for example \verb"/usr" is your default installation directory,
%\verb"make install"
%will copy the NFFT~3 library to \verb+/usr/lib+, the C header files to
%\verb+/usr/include/nfft+ and the documentation together with all examples and
%applications to \verb+/usr/share/nfft+.

For more information about options offered by the \verb"configure" script,
run the script with the option \verb+--help+. For example, to change the
installation location you can use the \verb"--prefix".
Here is a short list of some important options:
By default, NFFT~3 uses the Kaiser-Bessel window function
(see Appendix \ref{appendix.nfft}). This can be changed
by passing the option \verb"--with-window=ARG",
where \verb"ARG" is replaced by one of \verb"kaiserbessel" (Kaiser-Bessel),
\verb"gaussian" (Gaussian), \verb"bspline" (B-spline), or
\verb"sinc" (sinc power). This makes all routines use the
respective window function.
NFFT~3 routines can also be configured to measure the computation time for
several steps during execution. You can enable this behavior with
the options
\verb"--enable-measure-time" and/or \verb"--enable-measure-time-fftw".

%-------------------------------------------------------------------------------
\section{NFFT - nonequispaced fast Fourier transform}
\label{appendix.nfft}
%-------------------------------------------------------------------------------
In this section, we describe the approximate NFFT algorithm in more detail.
To illustrate the principles, we restrict ourselves to the one-dimensional
case, hence the computation of
\begin{align}
  \label{trigPoly1d}
  f\left(x\right) & = \sum_{k\in I_{N}} \hat{f}_k \eim{k x} &
  (j=0,\hdots,M-1)
\end{align}
for nonequispaced nodes $x_j \in \T$.
In this case, the NFFT has cost
${\cal O}\left(N \log N  + |\log \varepsilon|
M\right)$, where $\varepsilon$ is the desired accuracy.
The chief idea of the NFFT algorithm is to use standard FFTs in combination with
an approximation scheme that is based on a window function $\varphi$. This function needs to
be mutually well localized in time/spatial and frequency domain.
%$\mathbb{R}$ and at the same time also in the frequency domain $\mathbb{R}$.
Several such window functions have been proposed in
\cite{duro93,bey95,st97,Fou02,fesu02}.

\subsection*{The ansatz}
The first step is to approximate the trigonometric polynomial $f$ of degree
$N$ in \eqref{trigPoly1d} by a linear combination $s_{1}$ of shifted
one-periodic window functions $\tilde \varphi$,
\begin{equation}\label{s1}
  s_1\left(x\right) := \sum_{l \in I_n} g_{l} \; \tilde \varphi\left(x -
  \frac{l}{n}\right) \,,
\end{equation}
where $n:=\sigma N$ for some oversampling factor $\sigma >1$. The value of $n$ determines the
length of the ordinary FFT used below.

\subsection*{The window function}
The window function $\varphi:\R\rightarrow\R$ is chosen so that
the one-periodic version
$
  \tilde \varphi\left(x\right):=\sum_{r \in \Z} \varphi\left(x+r\right)
$
has a uniformly convergent Fourier series
$
  \tilde \varphi\left(x\right)=
  \sum_{k\in \Z} c_k\left(\tilde\varphi\right) \eim{k x},
$
and at the same time is well localized in the time/spatial domain $\T$
and in the frequency domain $\mathbb{Z}$.
The Fourier coefficients $c_k\left( \tilde \varphi \right)$ of this one-periodic version are samples
of the Fourier transform $\hat \varphi$ of the original window function at integers
points, i.e.,
\begin{align*}
  c_k\left( \tilde \varphi \right)
  & :=\int\limits_{\T} \tilde \varphi \left(x\right) \eip{k x} \, {\rm d} x
  = \int\limits_{\mathbb{R}} \varphi \left(x\right) \eip{k x} \, {\rm d} x
  =: \hat \varphi \left(k\right) & (k\in\Z).
\end{align*}
An example is shown in Figure \ref{fig:window_fun}.
\begin{figure}[tp]
  \centering
   \subfigure[Window function.]{
     \includegraphics[width=4.8cm,height=3.8cm]{images/window_fct1}}
   \qquad
     %\includegraphics[width=4.8cm,height=3.8cm]{images/window_fct2} &
   \subfigure[Fourier transform.]{
     \includegraphics[width=4.8cm,height=3.8cm]{images/window_fct3}}
   \caption{Example for a suitable window function, in this case the
     Gaussian window function $\varphi$ (left)
     %its one-periodic version $\tilde \varphi$ (center),
     and the sampled integral Fourier-transform $\hat \varphi$ (right) with
     pass, transition, and stop band for $N=24,\,\sigma=\frac{4}{3},\,n=32$.}
   \label{fig:window_fun}
 \end{figure}

\subsection*{The first approximation - cut-off in frequency domain}
If we take the definition in \eqref{s1} and switch to the frequency domain, we obtain
\begin{equation*}
  s_1\left(x\right)
  =\sum_{k \in I_n} \hat g_k \, c_k\left(\tilde \varphi\right) \, \eim{k x}
  +\sum_{r \in \Z \backslash \left\{0\right\}} \sum_{k \in I_n} \hat g_k \,
  c_{k + n r} \left(\tilde \varphi\right) \, \eim{(k + n r)x }
\end{equation*}
with discrete Fourier coefficients
\begin{equation} \label{coeff}
  \hat g_k = \sum_{l \in I_n} g_l \, \eip{\frac{k l}{n}}.
\end{equation}
A comparison of \eqref{trigPoly1d} and \eqref{s1} suggests to
define
\begin{equation} \label{coeff2}
  \hat g_k :=
    \begin{cases}
      \frac{\hat f_k}{c_k \left(\tilde \varphi\right)} & \text{for $k \in I_N$},
      \\[1ex]
      0 & \text{for $k \in I_n \backslash I_N$},
    \end{cases}
\end{equation}
assuming that $c_k\left(\tilde\varphi\right)$ is small for $|k| \ge n-\frac{N}{2}$.
The values $g_l$ can now be computed by applying the Fourier inversion theorem
to \eqref{coeff}. This step can be realized by an FFT of length $n$.
This first approximation causes an aliasing error.

\subsection*{The second approximation - cut-off in time/spatial domain}
If the window function $\varphi$ is well localized in time/spatial domain then it can be
approximated by a truncated version
\begin{equation*}
  \psi\left(x\right):=\varphi\left(x\right) \chi_{[-{\frac{m}{n}},
  {\frac{m}{n}}]}\left(x\right),
\end{equation*}
with ${\rm supp} \, \psi = \left[-{\frac{m}{n}},{\frac{m}{n}}\right], \;m \ll
n,\; m\in \N$.
Here too, we can define a one periodic version $\tilde \psi$ with compact support
in $\T$ by
$
  \tilde \psi\left(x\right)=\sum_{r \in \Z} \psi\left(x+r\right).
$
With the help of the index set
$
  I_{n,m} \left(x_j\right) := \left\{ l \in I_n : n x_j - m \le l \le n x_j +
    m \right\}
$
a further approximation to the first approximation $s_1$ at the nodes $x_{j}$ can be defined by
\begin{equation} \label{approx}
  s\left(x_j\right) :=
  \sum_{l \in I_{n,m}\left(x_j\right) } g_l \,
  \tilde\varphi\left(x_j - {\frac{l}{n}}\right) =
  \sum_{l \in I_n } g_l \,
  \tilde\psi\left(x_j - {\frac{l}{n}}\right) \, .
\end{equation}
Note, that for fixed $x_j\in \T$, the above sum contains at most
$(2m+1)$ nonzero summands.
This second approximation causes a truncation error.

\subsection*{The case $d>1$}
The described approximation method can be generalized to arbitrary dimensions.
But first, a few generalizations have to be introduced.
A multivariate window function is now defined by
$\varphi\left(\mathbf{x}\right):=\prod_{t=0}^{d-1} \varphi_t\left(x_t\right)$
with a univariate window function $\varphi_t$.
A simple consequence is that the Fourier coefficients have the form $c_{\mathbf{k}}\left(\tilde \varphi\right)  =
\prod_{t=0}^{d-1} c_{k_t}\left(\tilde \varphi_t\right)$.
The ansatz function $s_{1}$ reads now
\begin{equation*}
  s_1\left(\mathbf{x}\right) := \sum_{\mathbf{l} \in I_{\mathbf{n}}}
 g_{\mathbf{l}} \; \tilde
  \varphi\left(\mathbf{x} - \mathbf{n}^{-1}\odot\mathbf{l}\right).
\end{equation*}
The size of the FFT size is determined to $\mathbf{n}:=\mathbf{\sigma} \odot
\mathbf{N}$, and the oversampling factor is
$\mathbf{\sigma}=\left(\sigma_0,\hdots,\sigma_{d-1}\right)^{\top}$.
Following \eqref{coeff2} we define
\begin{equation*}
  \hat g_{\mathbf{k}} := \left\{
    \begin{array}{ll}
      \frac{\hat f_{\mathbf{k}}}{c_{\mathbf{k}} \left(\tilde \varphi\right)}
 & \text{for
      } \mathbf{k} \in I_{\mathbf{N}} , \\[1ex]
      0 & \text{for } \mathbf{k} \in I_{\mathbf{n}} \backslash I_{\mathbf{N}}.
    \\
    \end{array}
  \right.
\end{equation*}
The values $g_{\mathbf{l}}$ are obtained by a (multivariate) FFT of size
$n_0\times n_1 \times \hdots \times n_{d-1}$.
Using the compactly supported function $\psi\left(\mathbf{x}\right)=\varphi
\left(\mathbf{x}\right)\chi_{[-{\frac{m}{n}},{\frac{m}{n}}]^d} \left(\mathbf{x}
\right)$ gives the final form,
\begin{equation*}
  s\left(\mathbf{x}_j\right) := \sum_{\mathbf{l} \in I_{\mathbf{n},m}
 \left(\mathbf{x}_j\right) }
  g_{\mathbf{l}} \, \tilde\varphi\left(\mathbf{x}_j - \mathbf{n}^{-1}\odot
 \mathbf{l}\right) =
\sum_{\mathbf{l} \in I_{\mathbf{n}} }
g_{\mathbf{l}} \, \tilde\psi\left(\mathbf{x}_j - \mathbf{n}^{-1}\odot
    \mathbf{l}\right)\,.
\end{equation*}
Again, $\tilde \psi$ denotes the one-periodic version of $\psi$.
Note that the multi-index set
$
  I_{\mathbf{n},m} \left(\mathbf{x}_j\right) := \left\{ \mathbf{l} \in
 I_{\mathbf{n}} : \mathbf{n}
  \odot \mathbf{x}_j - m \mathbf{1} \le \mathbf{l} \le \mathbf{n} \odot
 \mathbf{x}_j + m
  \mathbf{1}\right\}
$
has at most $(2m+1)^d$ entries for every $\mathbf{x}_j$.

\subsection*{The final algorithm}
The NFFT algorithm is given in pseudo-code in Algorithm~\ref{algo:nfft}. In
addition to the evaluation of the window functions, it requires roughly
$|I_{\mathbf{N}}|+|I_{\mathbf{n}}| \log |I_{\mathbf{n}}| + 2(2m+1)^d M$
floating point operations.
In matrix-vector notation, Algorithm~\ref{algo:nfft} can be written as
$
  \mathbf{A} \mathbf{\hat{f}} \approx \mathbf{D} \mathbf{F} \mathbf{B}
  \mathbf{\hat{f}},
$
where $\mathbf{B}$ denotes the sparse real $M \times
\left|I_{\mathbf{n}}\right|$
matrix
\begin{equation*} \label{spB}
  \mathbf{B} := \left( \tilde\psi\left(\mathbf{x}_j - \mathbf{n}^{-1}\odot
  \mathbf{l}\right)
  \right)_{j=0,\hdots,M-1;\, \mathbf{l}\in I_{\mathbf{n}}}.
\end{equation*}
The matrix
$\mathbf{F}:=(\eim{\mathbf{k}(\mathbf{N^{-1}}\odot\mathbf{j})})_{\mathbf{j},
  \mathbf{k}\in I_{\mathbf{n}}}$ is the ordinary Fourier matrix of size
$\left|I_{\mathbf{n}}
\right| \times
\left|I_{\mathbf{n}}\right|$ and $\mathbf{D}$ is a real $\left|I_{\mathbf{n}}
\right|
\times |I_{\mathbf{N}}|$ ``diagonal'' matrix defined by
\begin{equation*}
  \mathbf{D} := \bigotimes_{t=0}^{d-1} \left(\mathbf{O}_t \, | \, \mbox{diag}
  \left(1/\,
      c_{k_t}\left(\tilde \varphi_t\right)\right)_{k_t\in I_{N_t}} \, | \,
      \mathbf{O}_t \right)^{\top}
\end{equation*}
with zero matrices $\mathbf{O}_t$ of size $N_t \times \frac{n_t-N_t}{2}$.
The corresponding adjoint matrix-vector product is given by
$
  \mathbf{A}^{\adj} \mathbf{\hat{f}} \approx \mathbf{B}^{\top}
\mathbf{F}^{\adj}
    \mathbf{D}^{\top} \mathbf{\hat{f}}
$.
With the help of the transposed index set
$
  I_{\mathbf{n},m}^{\top}(\mathbf{l}) := \{j=0,\hdots,M-1 :
  \mathbf{l}-m
  \mathbf{1} \le \mathbf{n} \odot \mathbf{x}_j \le \mathbf{l}+m \mathbf{1}
  \}
$,
one obtains Algorithm~\ref{algo:nfftH} for the adjoint NFFT.
\begin{algorithm}[tp]
  \hrule
  \footnotesize
  \rule[-1.5ex]{0ex}{5ex}\raggedright Input: $d$, $M \in \N$, $\mathbf{N} \in 2\N^d$,
      ${\mathbf{x}}_j \in [-\frac{1}{2},\frac{1}{2})^d$, $j=0,\ldots,M-1$, and
      $\hat  f_{\mathbf{k}}\in \C$, $\mathbf{k} \in I_{\mathbf{N}}$.
  \begin{algorithmic}[1]
    \itemsep=1.1ex
    \STATE For $\mathbf{k} \in I_{\mathbf{N}}$ compute
      $\hat g_{\mathbf{k}} :=
      |I_{\mathbf{n}}|^{-1} \cdot \hat f_{\mathbf{k}}/
       c_{\mathbf{k}}(\tilde \varphi)$.

    \STATE For $\mathbf{l} \in I_{\mathbf{n}}$ compute
    $g_{\mathbf{l}} := \sum_{\mathbf{k} \in I_{\mathbf{N}}} \hat
    g_{\mathbf{k}} \, \eim{\mathbf{k} (\mathbf{n}^{-1} \odot \mathbf{l})}$
    by a $d$-variate FFT.

    \STATE For $j=0,\hdots,M-1$ compute
    $f_j := \sum_{\mathbf{l} \in I_{\mathbf{n},m}
      ({\mathbf{x}}_j)} g_{\mathbf{l}} \,
      \tilde\varphi({\mathbf{x}}_j - \mathbf{n}^{-1}\odot
      \mathbf{l})$.\\[2ex]
  \end{algorithmic}

  Output: Approximate function values $f_j,\;j=0,\hdots,M-1$.\\
  Arithmetic cost: $|I_{\mathbf{N}}|+|I_{\mathbf{n}}| \log |I_{\mathbf{n}}| +
  2(2m+1)^d M$ + evaluations of the window function.\rule[-1.5ex]{0ex}{3ex}
  \hrule
  \caption{Nonequispaced fast Fourier transform (NFFT).\label{algo:nfft}}
\end{algorithm}
\begin{algorithm}[tp]
  \hrule
  \footnotesize
  \rule[-1.5ex]{0ex}{5ex}\raggedright Input: $d$, $M \in \N$, $\mathbf{N}\in 2\N^d$,
      ${\mathbf{x}}_j\in [-\frac{1}{2},\frac{1}{2})^d$,
      and $f_j\in \C$, $j=0,\hdots,M-1$.\\[\baselineskip]
  \begin{algorithmic}[1]
    \itemsep=1.1ex
    \STATE For $\mathbf{l} \in I_{\mathbf{n}}$ compute
      $g_{\mathbf{l}} := \sum\limits_{j \in
      I_{\mathbf{n},m}^{\top} (\mathbf{l}) } f_j \,
      \tilde\varphi({\mathbf{x}}_j - \mathbf{n}^{-1}\odot \mathbf{l})$.
    \STATE For $\mathbf{k} \in I_{\mathbf{N}}$ compute
    $\hat g_{\mathbf{k}} := \sum\limits_{\mathbf{l} \in
      I_{\mathbf{n}}} g_{\mathbf{l}} \, {\rm e}^{+2 \pi {\rm i} \mathbf{k}
      (\mathbf{n}^{-1} \odot \mathbf{l})}$
    by $d$-variate (backward) FFT.
    \STATE For $\mathbf{k} \in I_{\mathbf{N}}$ compute
      $\hat h_{\mathbf{k}} := |I_{\mathbf{n}}|^{-1} \cdot
      \hat{g}_{\mathbf{k}} / c_{\mathbf{k}}(\tilde \varphi)$.\\[2ex]
  \end{algorithmic}

  Output: Approximate coefficients $\hat h_{\mathbf{k}},\;\mathbf{k} \in
    I_{\mathbf{N}}$.\\
 Arithmetic cost: $|I_{\mathbf{N}}|+|I_{\mathbf{n}}| \log |I_{\mathbf{n}}| +
  2(2m+1)^d M$ + evaluation of the window function.\rule[-1.5ex]{0ex}{3ex}
  \hrule
  \caption{Adjoint nonequispaced fast Fourier transform (adjoint NFFT).}
  \label{algo:nfftH}
\end{algorithm}

%Owing to the characterization of the nonzero elements in the matrix $\mathbf{B}$
%by
%$
%  \bigcup_{j=0}^{M-1} \{j\} \times I_{\mathbf{n},m}({\mathbf{x}}_j)
%  =
%  \bigcup_{\mathbf{l} \in I_{\mathbf{n}}} I_{\mathbf{n},m}^{\top}(\mathbf{l})
%  \times \{\mathbf{l}\},
%$
%the multiplication with $\mathbf{B}^{\top}$ is implemented in
%a ``transposed'' way in the library with summation as outer loop and only using
%the multi-index sets $I_{\mathbf{n},m} ({\mathbf{x}}_j)$.

%-------------------------------------------------------------------------------
\section{Available window functions and evaluation techniques}
\label{sect:windows}
%-------------------------------------------------------------------------------
For clarity of the presentation, we restrict ourselves to the case $d=1$, again.
To keep the aliasing error and the truncation error small, several functions
$\varphi$ with good localization in time and frequency domain have been
proposed,
e.g. the (dilated) \emph{Gaussian} \cite{duro93,st97,DuSc} (with shape
parameter $b := \frac{2\sigma}{2 \sigma -1} \frac{m}{\pi}$),
\begin{align}
\label{gaussian}
\varphi\left(x\right)
& =
\left(\pi b\right)^{-1/2} \,
{\rm e}^{-\frac{\left(nx\right)^2}{b}}
, &
\hat \varphi \left(k\right)
& =
\frac{1}{n} \, {\rm e}^{-b\left(\frac{\pi k}{n}\right)^2},
\intertext{
(dilated) \emph{cardinal central $B$--splines} \cite{bey95,st97},}
\label{bspline}
\varphi\left(x\right)
& =
M_{2m}\left(n x\right), &
\hat \varphi \left(k\right)
& =
\frac{1}{n} \sinc^{2m} \left( k \pi/n\right),
\intertext{where $M_{2m}$ denotes the centered cardinal $B$--Spline of order $2m$,
(dilated) \emph{Sinc functions} \cite{pottshabil},}
\label{sinc}
\varphi\left(x\right)
& =
\frac{N\left(2\sigma-1\right)}{2m} \sinc^{2m} \left(\frac{\left(\pi N x
\left(2\sigma-1\right)\right)}{2m} \right), &
\hat \varphi \left(k\right)
& =
M_{2m}\left(\frac{2mk}{\left(2\sigma-1\right)N}\right),
\end{align}
and
(dilated) \emph{Kaiser--Bessel functions} \cite{Ja,Fou02} (with shape
parameter $b :=  \pi \left(2 - \frac{1}{\sigma} \right)$),
\begin{equation}
  \label{kaiser}
  \begin{aligned}
    \varphi\left(x\right)
    & = \frac{1}{\pi}
    \begin{cases}
      \ds \frac{\sinh\left(b\sqrt{m^2 - n^2 x^2}\right)}{\sqrt{m^2 - n^2 x^2}},
      & \text{for $|x| \le \frac{m}{n}$},
      \\[4ex]
      \ds \frac{\sin\left(b\sqrt{n^2 x^2 - m^2}\right)}{\sqrt{n^2 x^2-m^2}},
      & \text{otherwise},
    \end{cases}
    \\[1ex]
    \hat \varphi \left(k\right)
    & = \frac{1}{n}
    \begin{cases}
      \ds I_0 \left(m \sqrt{b^2 - \left(2 \pi k/n \right)^2} \right),
      & \text{for $k= -n \left(1-\frac{1}{2\sigma} \right), \ldots,n \left(1-\frac{1}{2\sigma}\right)$},
      \\
      0, & \text{otherwise},
    \end{cases}
  \end{aligned}
\end{equation}
where $I_0$ denotes the \emph{modified order-zero Bessel function}.
For these functions $\varphi$ it has been proved that
\[
|f\left(x_j\right) - s\left(x_j\right)| \le
C\left(\sigma,m\right) \sum_{k\in I_N} |\hat{f}_k|,
\]
where
\[
C\left(\sigma,m\right) :=
\begin{cases}
  4 \, \e^{-m \pi \left(1-1/\left(2\sigma-1\right)\right)}, & \text{for \eqref{gaussian}},
  \\
  4 \, \left(\frac{1}{2\sigma -1} \right)^{2m}, & \text{for \eqref{bspline}},
  \\[1ex]
  \frac{1}{m-1}\left(\frac{2}{\sigma^{2m}}+\left(\frac{\sigma}{2\sigma-1}\right)^{2m}\right), &
    \text{for \eqref{sinc}},
  \\[1ex]
  4\pi\left(\sqrt{m}+m\right)\sqrt[4]{1-\frac{1}{\sigma}} \, \e^{-2 \pi m \sqrt{1-1/\sigma}}, &
    \text{for \eqref{kaiser}}.
\end{cases}
\]
Thus, for fixed $\sigma >1$ the approximation error introduced by the NFFT
decays exponentially with the number of summands $m$ in \eqref{approx}.
These estimates can of course be generalized to higher dimensions; see \cite{elst}.

It must be noted that the cost of the NFFT also increases with $m$.
In the following, we will suggest different methods for efficient storage and
application of the matrix $\mathbf{B}$. These are all available in the
NFFT~3 library by choosing particular flags during the initialization phase.
%These methods do not yield a different asymptotic performance but rather yield
%a lower constant in the amount of computation.

\subsection*{Fully precomputed window function}
An obvious possibility is to precompute all values $\varphi(\mathbf{x}_j -
\mathbf{n}^{-1}\odot\mathbf{l})$ for $j=0,\hdots,M-1$ and $\mathbf{l}\in
I_{\mathbf{n},m} (\mathbf{x}_j)$ which are needed for $\mathbf{B}$.
This needs to store the possibly large amount of $(2m+1)^d M$ real numbers.
On the other hand, no additional computation is required during the matrix-vector
multiplication apart from the inevitable $2(2m+1)^d M$ floating point operations.
%Furthermore, we store for this method explicitly the row and column index for
%each nonzero entry of the matrix $\mathbf{B}$.
This method, indicated by the flag \texttt{PRE\_FULL\_PSI}, is the fastest procedure
but can only be used if enough memory for storage is available.

\subsection*{Tensor product based precomputation}
Using the fact that the window functions are represented by tensor products one can
resort to storing only the values $\varphi_t((\mathbf{x}_j)_t - \frac{l_t}{n_t})$ for
$j=0,\hdots,M-1$, $t=0,\hdots,d-1$, and $l_t\in I_{n_t,m} ((\mathbf{x}_j)_t)$,
where $(\mathbf{x}_j)_t$ denotes the $t$th component of the $j$th node.
This method needs to store $d(2m+1)M$ real numbers which is less than the first method needs.
However, for each node one possibly needs up to $2(2m+1)^d$ extra multiplications to
obtain the values of the multivariate window function
$\varphi(\mathbf{x}_j - \mathbf{n}^{-1}\odot\mathbf{l})$ for $\mathbf{l}\in
I_{\mathbf{n},m}(\mathbf{x}_j)$.
This technique is available for every window function discussed
and can be used by setting the flag \texttt{PRE\_PSI}. This is also the
default method used in our library.

\subsection*{Linear interpolation from a lookup table}
For a large number of nodes $M$, the amount of memory can be further reduced
by using lookup table techniques.
For a recent example within the framework of gridding see \cite{BeNiPa05}.
We suggest to precompute from the even window function the equidistant samples
$\varphi_t(\frac{rm}{Kn_t})$ for $t=0,\hdots,d-1$ and $r=0,\hdots,K$,
$K\in\N$. Then for the actual node $\mathbf{x}_j$ during the NFFT
the values $\varphi_t((\mathbf{x}_j)_t - \frac{l_t}{n_t})$ for
$t=0,\hdots,d-1$ and $l_t\in I_{n_t,m} ((\mathbf{x}_j)_t)$ are computed
by linear interpolation from its two neighboring values.
This method needs to store only $dK$ real numbers in total, where $K$
depends only on the desired accuracy but neither on the number of
nodes $M$ nor on the multibandwidth $\mathbf{N}$.
If we choose $K$ to be a multiple of $m$, we can further reduce the computational
costs during the interpolation step since the distance from $(\mathbf{x}_j)_t-
\frac{l_t}{n_t}$ to the two neighboring values is the same for all
$l_t\in I_{n_t,m}((\mathbf{x}_j)_t)$.
During the actual multiplication with the matrix $\mathbf{B}$, this method requires
$2(2m+1)^d$ extra multiplications per node. It can be used by setting the flag
\texttt{PRE\_LIN\_PSI}. Error estimates for this additional
approximation are given in \cite{KuPo06}.

\subsection*{Fast Gaussian gridding}
Two useful properties of the Gaussian window function \eqref{gaussian} that can be exploited
within the presented framework have recently been reviewed in \cite{GrLe04}.
Beside the tensor product structure for $d>1$, which is also valid for all other
window functions, it is remarkable that the number of evaluations of the
exponential function \texttt{exp()} can be reduced by a substantially.
More precisely, for $d=1$ and a fixed node $x_j$ the evaluation of
$\varphi(x_j-\frac{l'}{n})$, $l'\in I_{n,m}(x_j)$ can be computed more
efficiently by writing this as
\begin{equation*}
  \sqrt{\pi b}\varphi\left(x_j-\frac{l'}{n}\right)
  =\e^{-\frac{(nx_j-l')^2}{b}}
  =\e^{-\frac{(nx_j-u)^2}{b}}
  \left(\e^{-\frac{2(nx_j-u)}{b}}\right)^l
  \e^{-\frac{l^2}{b}}\,,
\end{equation*}
where $u=\min I_{n,m}(x_j)$ and $l=0,\hdots,2m$.
The first factor and the exponential within the brackets are a
constant for each fixed node $x_j$.
Once we have evaluated the second exponential, its $l$th power can be computed
by repeated multiplications only.
Furthermore, the last exponential is independent of $x_j$ such that these $2m+1$
different values need to be precomputed only once -- usually a negligible amount.
Thus, it is sufficient to store or evaluate $2M$ exponentials for $d=1$.
The case $d>1$ needs to store or evaluate $2dM$ numbers. This follows from the
tensor product structure.
The described method is employed when the flags \texttt{FG\_PSI} and possibly
\texttt{PRE\_FG\_PSI} are used.

\subsection*{No precomputation of the window function}
The last method considered does not precompute any values of the window function
at all. Instead, all values are calculated on-line. This needs to be done at most
$\frac{3}{2}(2m+1)^d M$ times. The computation time strongly depends on how fast
one can evaluate the particular window function. Since this method does not need
any additional storage, it is attractive whenever large problems are considered
that would push the other methods beyond the memory limitations.

%-------------------------------------------------------------------------------
\section{Further NFFT approaches}\label{A.4}
%-------------------------------------------------------------------------------
Besides \cite{nfft3}, the {\sc MatLab} toolbox \cite{fematlab} also implements
the NFFT.
In several papers, approximations techniques similar to the NFFT have been
described.
Common names for these schemes are
\emph{non-uniform fast Fourier transform} \cite{fesu02},
\emph{generalized fast Fourier transform} \cite{duro93},
\emph{unequally-spaced fast Fourier transform} \cite{bey95},
\emph{fast approximate Fourier transforms for irregularly spaced data}
\cite{ware98},
\emph{non-equispaced fast Fourier transform} \cite{Fou02} or
\emph{gridding} \cite{scsc,Ja,Pe}.

Also, different window functions have been considered in various papers, e.g. a Gaussian
pulse tapered with a Hanning window \cite{DuSc}, Gaussian kernels combined
with Sinc kernels \cite{Pe}, and particular optimized windows \cite{Ja,DuSc}.
Special approaches based on scaling vectors \cite{NgLi99}, minimizing
the Frobenius norm of certain error matrices \cite{st01}, or min-max
interpolation \cite{fesu02} have been proposed.
While these approaches lead to better accuracy for the Gaussian or B-Spline
windows, improvements for the Kaiser-Bessel window are marginal.
Recently, the time consuming multiplication with the matrix $\mathbf{B}$ has
been accelerated by using commodity graphics hardware \cite{SoScNoHa08}.

The following two approaches have also been considered for the
univariate case \hbox{$d=1$:}
The authors of \cite{AnDa96} use for each node $x_j$ an $m$th order Taylor
expansion of the trigonometric polynomial \eqref{trigPoly1d} around the nearest
neighboring point on an oversampled equispaced lattice $I_n/n$. Here again we have
$n=\sigma N,\, \sigma\gg 1$.
This algorithm utilizes $m$ FFTs of size $n$ (compared to only one in our approach)
and uses a fair amount of extra memory; see also \cite{KuPo06,Ku08GAMM}.
Another approach is considered in \cite{DuRo95}. It is based on a Lagrange
interpolation technique.
After computing an FFT of the vector $\mathbf{\hat{f}}\in\C^N$ in \eqref{eq:nfft}
an exact polynomial interpolation scheme is used to obtain the values of the
trigonometric polynomial $f$ at the nonequispaced nodes $x_j$.
This is the most time consuming part of this method. In an approximate way,
it can however be realized with the help of the fast multipole method.
This approach is appealing since it allows also for computing an
inverse transform.
Nevertheless, numerical experiments in \cite{DuRo95} indicate that this
approach is far more time consuming than Algorithm \ref{algo:nfft}. Also, the
inversion can only be computed in a stable way for almost equispaced nodes
\cite{DuRo95}.

For a comparison of different approaches to nonequispaced Fourier transforms, we also
refer to \cite{ware98,st01,fesu02,KuPo06}.

\begin{acks}
We thank G.~Steidl for numerous fruitful and enlightning discussions.
Contributions to the source code made by M. B\"ohme, M. Fenn,
T. Knopp, and S. Klatt are gratefully acknowledged. For details, see the
\verb"AUTHORS" file in the package directory. We also thank the anonymous
referees and the editor for man helpful comments and valuable suggestions.
\end{acks}

\bibliographystyle{acmtrans}
%\bibliographystyle{abbrv}
\bibliography{../references/references}
\begin{received}
Received 09/2007; revised 02/2008; accepted Month Year
\end{received}
\end{document}
