%==============================================================================
\documentclass[11pt,a4paper,twoside,bibtotoc]{scrartcl}
%==============================================================================

% LaTeX packages
\usepackage{a4wide}
\usepackage{amsfonts}
\usepackage{amsmath}
\usepackage{theorem}
\usepackage{color}
\usepackage{graphicx}
\usepackage{subfigure}
\usepackage[boxed,Algorithm]{algorithm}
\usepackage{algorithmic}

% Allow for more figures on each side
\renewcommand{\topfraction}{1}
\renewcommand{\textfraction}{0}
\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{\nfftversion}{3.0}
\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}
}

\def\invisible#1{\textcolor{white}{#1}}

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

%==============================================================================
% Title Page
%==============================================================================
\begin{document}

\title{NFFT \nfftversion\ - Tutorial}
\author{Jens Keiner\thanks{keiner@math.uni-luebeck.de,University of
    L\"ubeck, Institute of  Mathematics, 23560 L\"ubeck} \and
  Stefan Kunis\thanks{kunis@mathematik.tu-chemnitz.de, Chemnitz University of
    Technology, Department of Mathematics, 09107 Chemnitz, Germany} \and
  Daniel Potts\thanks{potts@mathematik.tu-chemnitz.de, Chemnitz University of
    Technology, Department of Mathematics, 09107 Chemnitz, Germany}}
\date{}
\maketitle

\begin{center}
  \vspace{-8ex}
  {\large {\tt http://www.tu-chemnitz.de/$\sim$potts/nfft}}
  \vspace{8ex}
\end{center}

%\begin{abstract}
%\noindent {\it Key words and phrases}:  \\
%\noindent {\it Last revision}: \today \\
%\noindent {\it 2000 AMS Mathematics Subject Classification}: 65T50, 65T40
%\end{abstract}


\begin{center}
  \begin{tabular}{cc}
    \includegraphics[width=7.5cm]{images/front_pic1} &
    \includegraphics[width=7.5cm]{images/front_pic2}
  \end{tabular}
\end{center}

\clearpage
\tableofcontents


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\clearpage
\section{Introduction}\label{sect:1}
\setcounter{equation}{0}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

This tutorial surveys the fast Fourier transform at nonequispaced nodes (NFFT),
its generalisations, its inversion, and the related C library NFFT \nfftversion.
The goal of this manuscript is twofold -- to shed some light on the
mathematical background, as well as to provide an overview over the C library
to allow the user to use it effectively. Following this, this document splits
into two parts:

In Sections \ref{sect:nfft}, \ref{sect:gen} and \ref{sect:inv}, we
introduce the NFFT which, as a starting point, leads to several concepts for
further generalisation and inversion. We focus on a mathematical
description of the related algorithms. However, we complement the necessarily
incomplete quick glance at the ideas by references to related mathematical
literature.
In Section \ref{sect:nfft}, we introduce the problem of computing the discrete
Fourier transform at nonequispaced nodes (NDFT), along with the notation
used and related works in the literature. The algorithms (NFFT) are developed
in pseudo code.
Furthermore, we turn for recent generalisations, including in particular NFFTs
on the sphere and iterative schemes for inversion of the nonequispaced FFTs,
in Section \ref{sect:gen} and \ref{sect:inv}, respectively.

Having laid the theoretical foundations, Section \ref{sect:lib} provides an
overview over the NFFT \nfftversion\ C library and the general principles for using the
available algorithms in your own code. Rather than describing the library
interface in every detail, we restrict to simple recipes in order to
familiarise with the very general concepts sufficient for most everyday tasks.
For the experienced user, we provide full interface documentation
(see \verb+doc/html/index.html+ in the package directory) which gives
detailed information on more advanced options and parameter settings for
each routine. Figure \ref{nfft:fig:svn} gives an overview over the directory
structure of the NFFT \nfftversion\ package.
Finally, Section \ref{sect:ex} gives some simple examples for using the
library and some more advanced applications are given in Section
\ref{sec:appl}
along with numerical results.

\begin{figure}[p]
  \centering
  \input{tree}
  \caption{Directory structure of the NFFT \nfftversion\ package
    \label{nfft:fig:svn}}
\end{figure}

%The computational complexity of the proposed algorithms is $\cO(|I_{\mathbf{N}}|
%\log |I_{\mathbf{N}}| + M)$ for $|I_{\mathbf{N}}|$ equispaced frequencies and $M$
%nonequispaced samples.
%Besides the core module for the nonequispaced FFT, the following tree
%structure shows various generalisations, a module for inverse transforms, some
%examples, and recent applications.

%The remainder of this tutorial is organised as follows:
%The implementation and usage of the library is described in Section
%\ref{sect:library}.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\clearpage
\section{Notation, the NDFT, and the NFFT}\label{sect:nfft}
\setcounter{equation}{0}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

This section summarises the mathematical theory and ideas behind the NFFT.
Let the torus
\begin{equation*}
  \T^d
  := \left\{ \mathbf{x}=\left(x_t\right)_{t=0,\hdots,d-1} \in \R^d: - \frac{1}{2}
  \le x_t < \frac{1}{2},\; t=0,\hdots,d-1 \right\}
\end{equation*}
of dimension $d\in\N$ be given.
It will serve as domain from which the nonequispaced nodes $\mathbf{x}$ are taken.
Thus, the sampling set is given by ${\cal X}:=\{\mathbf{x}_j \in
\T^d:\,j=0,\hdots,M-1\}$.

Possible frequencies $\mathbf{k}\in\Z^d$ are collected in the multi-index set
\begin{equation*}
  I_{\mathbf{N}} := \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*}
where $\mathbf{N}=\left(N_t\right)_{t=0,\hdots,d-1}$ is the {\bf EVEN}
multibandlimit, i.e., $N_t\in 2\N$.
To keep notation simple, the multi-index $\mathbf{k}$ addresses elements of
vectors and matrices as well, i.e., the plain index $\tilde
k:=\sum_{t=0}^{d-1} (k_t+\frac{N_t}{2}) \prod_{t'=t+1}^{d-1} N_{t'}$
is not used here.
The inner product between the frequency index $\mathbf{k}$ and the time/spatial
node $\mathbf{x}$ is defined in the usual way by $\mathbf{k} \mathbf{x}:=k_0 x_0 + k_1 x_1
+\hdots+ k_{d-1} x_{d-1}$.
Furthermore, two vectors may be combined by the component-wise product $\mathbf{\sigma}
\odot \mathbf{N}:=\left(\sigma_0 N_0, \sigma_1 N_1, \hdots, \sigma_{d-1} N_{d-1},
\right)^{\top}$ with its inverse $\mathbf{N}^{-1}:=\left(\frac{1}{N_0},
  \frac{1}{N_1}, \hdots,  \frac{1}{N_{d-1}} \right)^{\top}$.

The space of all $d$-variate, one-periodic functions $f: \T^d \rightarrow \C$
is restricted to the space of $d$-variate trigonometric polynomials
\begin{equation*}
  T_{\mathbf{N}}:={\rm span}\left(\eim{\mathbf{k}\mathbf{\cdot}}:\,\mathbf{k} \in I_{\mathbf{N}}\right)
\end{equation*}
with degree $N_t\; (t=0,\hdots,d-1)$ in the $t$-th dimension.
The dimension $\dim T_{N}$ of the space of $d$-variate trigonometric
polynomials $T_{N}$ is given by $\dim T_{N} = |I_{\mathbf{N}}| =
\prod\limits_{t=0}^{d-1} N_t$.

%------------------------------------------------------------------------------
\subsection{NDFT - nonequispaced discrete Fourier transform}
%------------------------------------------------------------------------------
The first problem to be addressed can be regarded as a matrix vector
multiplication. For a finite number of given Fourier coefficients
$\hat f_{\mathbf{k}} \in \C$, $\mathbf{k}\in I_{\mathbf{N}}$, we consider the evaluation of the
trigonometric polynomial
\begin{equation} \label{trigPoly}
  f\left(\mathbf{x}\right) := \sum_{\mathbf{k}\in I_{\mathbf{N}}} \hat{f}_{\mathbf{k}}
  \eim{\mathbf{k}\mathbf{x}}
\end{equation}
at given nonequispaced nodes $\mathbf{x}_j \in \T^d$.
Thus, our concern is the evaluation of
\begin{equation}\label{eq:nfft}
  f_j = f\left(\mathbf{x}_j\right) := \sum_{\mathbf{k}\in I_{\mathbf{N}}} \hat{f}_{\mathbf{k}}
  \eim{\mathbf{k}\mathbf{x}_j},
\end{equation}
$j=0,\hdots,M-1$.
In matrix vector notation this reads
\begin{equation} \label{ndft}
  \mathbf{f}=\mathbf{A} \mathbf{\hat{f}}
\end{equation}
where
\begin{equation*}
  \mathbf{f}:=\left(f_j\right)_{j=0,\hdots,M-1},\quad
  \mathbf{A}:=\left(\eim{\mathbf{k}\mathbf{x}_j}\right)_{j=0,\hdots,M-1;\;\mathbf{k}\in I_{\mathbf{N}}},
  \quad
  \mathbf{\hat{f}}:=\left(\hat{f}_{\mathbf{k}}\right)_{\mathbf{k}\in I_{\mathbf{N}}}.
\end{equation*}
The straightforward algorithm for computing this matrix vector product,
which is called NDFT, takes $\cO(M |I_{\mathbf{N}}|)$ arithmetical operations.

A closely related matrix vector product is the adjoint NDFT
\begin{equation*}
  \mathbf{\hat{h}}=\mathbf{A}^{\adj} \mathbf{f}, \qquad \hat h_{\mathbf{k}}=\sum_{j=0}^{M-1}
  \hat{f}_{ j} \eip{\mathbf{k}\mathbf{x}_j},
\end{equation*}
where $\mathbf{A}^{\adj}=\overline {\mathbf{A}}^{\top}$ denotes the conjugate transpose
of the nonequispaced Fourier matrix $\mathbf{A}$.

\subsubsection*{Equispaced nodes}
For $d\in\N,\;N_t=N\in 2\N$, $t=0,\hdots,d-1$ and $M=N^d$ equispaced nodes
$\mathbf{x}_{\mathbf{j}}={1 \over N}\mathbf{j}$, $\mathbf{j} \in I_{\mathbf{N}}$ the computation of
\eqref{ndft} is known as multivariate discrete Fourier transform (DFT).
In this special case, the input data $\hat{f}_{\mathbf{k}}$ are called discrete
Fourier coefficients and the samples $f_{\mathbf{j}}$ can be computed by the well
known fast Fourier transform (FFT) with only $\cO(|I_{\mathbf{N}}| \log |
  I_{\mathbf{N}}|)$ arithmetic operations.
Furthermore, one has the inversion formula
\begin{equation*}
  \mathbf{A} \mathbf{A}^{\adj}=\mathbf{A}^{\adj} \mathbf{A}= |I_{\mathbf{N}}| \mathbf{I},
\end{equation*}
which does {\bf NOT} hold true for the nonequispaced case in general.

%------------------------------------------------------------------------------
\subsection{NFFT - nonequispaced fast Fourier transform}
\label{sec.nfft}
%------------------------------------------------------------------------------
For the sake of simplicity, we explain the ideas behind the NFFT for the
one-dimensional case $d=1$ and the algorithm NFFT.
The generalisation of the FFT is an approximative algorithm and has
computational complexity ${\cal O}\left(N \log N  + \log
  \left(1/\varepsilon\right) M\right)$, where $\varepsilon$ denotes the
desired accuracy.
The main idea is to use standard FFTs and a window function $\varphi$ which is
well localised in the time/spatial domain $\mathbb{R}$ and in the frequency
domain $\mathbb{R}$.
Several window functions were proposed in \cite{duro93,bey95,st97,four,fesu02}.

The considered problem is the fast evaluation of
\begin{equation} \label{trigPoly1d}
  f\left(x\right) = \sum_{k\in I_N} \hat{f}_k \eim{k x}
\end{equation}
at arbitrary nodes $x_j \in \T,\;j=0,\hdots,M-1$.

\subsubsection*{The ansatz}
One wants to approximate the trigonometric polynomial $f$ in
(\ref{trigPoly1d}) by a linear combination of shifted 1-periodic window
functions $\tilde \varphi$ as
\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}
With the help of an oversampling factor $\sigma >1$, the FFT length is given by
$n:=\sigma N$.

\subsubsection*{The window function}
Starting with a reasonable window function $\varphi:\R\rightarrow\R$, one
assumes that its 1-periodic version $\tilde\varphi$, i.e.
\begin{equation*}
  \tilde \varphi\left(x\right):=\sum_{r \in \Z} \varphi\left(x+r\right)
\end{equation*}
has an uniformly convergent Fourier series and is well localised in the
time/spatial domain $\T$ and in the frequency domain $\mathbb{Z}$.
The periodic window function $\tilde \varphi$ may be represented by its
Fourier series
\begin{equation*}
  \tilde \varphi\left(x\right)=
  \sum_{k\in \Z} c_k\left(\tilde\varphi\right) \eim{k x}
\end{equation*}
with the Fourier coefficients
\begin{equation*}
  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),\quad k\in\Z.
\end{equation*}
\begin{figure}[ht!]
   \begin{center}
    \begin{tabular}{ccc}
     \includegraphics[width=4.8cm,height=3.8cm]{images/window_fct1} &
     \includegraphics[width=4.8cm,height=3.8cm]{images/window_fct2} &
     \includegraphics[width=4.8cm,height=3.8cm]{images/window_fct3}
    \end{tabular}
   \end{center}
   \caption{From left to right: Gaussian window function $\varphi$, its
     1-periodic version $\tilde \varphi$, and the integral Fourier-transform
     $\hat \varphi$ (with pass, transition, and stop band) for
     $N=24,\,\sigma=\frac{4}{3},\,n=32$.}
   \label{fig:window_fun}
 \end{figure}

\subsubsection*{The first approximation - cut-off in frequency domain}
Switching from the definition (\ref{s1}) to the frequency domain, one obtains
\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 the discrete Fourier coefficients
\begin{equation} \label{coeff}
  \hat g_k := \sum_{l \in I_n} g_l \, \eip{\frac{k l}{n}}.
\end{equation}
Comparing (\ref{trigPoly1d}) to (\ref{s1}) and assuming
$c_k\left(\tilde\varphi\right)$ small for $|k|\ge n-\frac{N}{2}$ suggests to
set
\begin{equation} \label{coeff2}
  \hat g_k := \left\{
    \begin{array}{ll}
      \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{array}
  \right.
\end{equation}
Then the values $g_l$ can be obtained from (\ref{coeff}) by
\begin{equation*}
  g_l = \frac{1}{n}\sum\limits_{k \in I_N} \hat g_k \, \eim{\frac{k l}{n}}
  \qquad (l \in I_n),
\end{equation*}
a FFT of size $n$.

This approximation causes an aliasing error.

\subsubsection*{The second approximation - cut-off in time/spatial domain}
If $\varphi$ is well localised in time/space domain $\mathbb{R}$ it can be
approximated by a function
\begin{equation*}
  \psi\left(x\right)=\varphi\left(x\right) \chi_{[-{m \over n},{m \over
      n}]}\left(x\right)
\end{equation*}
with ${\rm supp} \, \psi \, \left[-{m \over n},{m \over n}\right], \;m \ll
n,\; m\in \N$.
Again, one defines its one periodic version $\tilde \psi$ with compact support
in $\T$ as
\begin{equation*}
  \tilde \psi\left(x\right)=\sum_{r \in \Z} \psi\left(x+r\right).
\end{equation*}
With the help of the index set
\begin{equation*}
  I_{n,m} \left(x_j\right) := \left\{ l \in I_n : n x_j - m \le l \le n x_j +
    m \right\}
\end{equation*}
an approximation to $s_1$ is defined by
\begin{equation} \label{approx}
  s\left(x_j\right) := \sum_{l \in I_{n,m}\left(x_j\right) } g_l \,
  \tilde\psi\left(x_j - {l \over n}\right) \, .
\end{equation}
Note, that for fixed $x_j\in \T$, the above sum contains at most
$(2m+1)$ nonzero summands.

This approximation causes a truncation error.

\subsubsection*{The case $d>1$}
Starting with the original problem of evaluating the multivariate
trigonometric polynomial in (\ref{trigPoly}) one has to do a few
generalisations.
The window function is given by
\begin{equation*}
  \varphi\left(\mathbf{x}\right):=\varphi_0\left(x_0\right)
  \varphi_1\left(x_1\right) \hdots \varphi_{d-1}\left(x_{d-1}\right)
\end{equation*}
where $\varphi_t$ is an univariate window function.
Thus, a simple consequence is
\begin{equation*}
  c_{\mathbf{k}}\left(\tilde \varphi\right)
  = c_{k_0}\left(\tilde \varphi_0\right) c_{k_1}\left(\tilde \varphi_1\right)
  \hdots c_{k_{d-1}}\left(\tilde \varphi_{d-1}\right).
\end{equation*}
The ansatz is generalised to
\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*}
where the FFT size is given by $\mathbf{n}:=\mathbf{\sigma} \odot \mathbf{N}$ and the
oversampling factors by
$\mathbf{\sigma}=\left(\sigma_0,\hdots,\sigma_{d-1}\right)^{\top}$.
Along the lines of (\ref{coeff2}) one defines
\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}}$ can be obtained by a (multivariate) FFT of size $n_0
\times n_1 \times \hdots \times n_{d-1}$ as
\begin{equation*}
  g_{\mathbf{l}} = \frac{1}{\left|I_{\mathbf{n}}\right|}\sum\limits_{\mathbf{k} \in
  I_{\mathbf{N}}} \hat g_{\mathbf{k}} \, \eim{\mathbf{k} \left(\mathbf{n}^{-1}\odot \mathbf{l}\right)},
  \quad \mathbf{l} \in I_{\mathbf{n}}.
\end{equation*}
Using the compactly supported function $\psi\left(\mathbf{x}\right)=\varphi\left(\mathbf{x}\right)\chi_{[-{m \over n},{m \over
    n}]^d}\left(\mathbf{x}\right)$,
one obtains
\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\psi\left(\mathbf{x}_j - \mathbf{n}^{-1}\odot\mathbf{l}\right) \, ,
\end{equation*}
where $\tilde \psi$ again denotes the one periodic version of $\psi$ and the
multi-index set is given by
\begin{equation*}
  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\}\,.
\end{equation*}

\subsubsection*{The algorithm}
In summary, the following Algorithm~\ref{algo:nfft} is obtained for the fast
computation of \eqref{ndft} with ${\cal O} \left(\left|I_{\mathbf{n}}\right| \log
  \left|I_{\mathbf{n}}\right| + m M\right)$ arithmetic operations.
\begin{algorithm}[ht]
  \caption{\tt NFFT\label{algo:nfft}}
  Input: $d,M \in \N$, $\mathbf{N}\in 2\N^d$\\
  {\invisible{Input:}} ${\mathbf{x}}_j\in [-\frac{1}{2},\frac{1}{2}]^d,\;
  j=0,\hdots,M-1$, and $\hat f_{\mathbf{k}}\in \C,\; \mathbf{k} \in I_{\mathbf{N}}$,\\[1ex]
                                %{\invisible{Input:}} ($\mathbf{\sigma} \in \mathbb{R}^d ,\;\sigma_t > 1,\;t=0,\hdots,d-1$, $m \in \N$, flags)\\[1ex]
                                %Precomputation: $\mathbf{n} := \sigma \odot \mathbf{N}$,\\
                                %{\invisible{Precomputation:}} $c_{\mathbf{k}}\left(\tilde \varphi\right), \; \mathbf{k} \in I_{\mathbf{N}}$
                                %if the flag (see Section \ref{sect:parameters}) {\tt PRE\_PHI\_HUT} is set,\\
                                %{\invisible{Precomputation:}} $\tilde\psi\left({\mathbf{x}}_j - \mathbf{n}^{-1}\odot \mathbf{l}\right),
                                %\, j=0,\hdots,M-1 , \, \mathbf{l} \in I_{\mathbf{n},\mathbf{m}} \left({\mathbf{x}}_j\right)$ if {\tt PRE\_PSI} is set.
  \begin{algorithmic}[1]
    \topsep=1.5ex
    \itemsep=1.5ex

    \STATE For $\mathbf{k} \in I_{\mathbf{N}}$ compute
    \begin{equation*}
      \hat g_{\mathbf{k}} :=\frac{\hat f_{\mathbf{k}}}{ \left|I_{\mathbf{n}}\right| c_{\mathbf{k}}
        \left(\tilde \varphi\right)}.
    \end{equation*}

    \STATE For $\mathbf{l} \in I_{\mathbf{n}}$ compute by $d$-variate {\rm FFT}
    \begin{equation*}
      g_{\mathbf{l}} := \sum\limits_{\mathbf{k} \in I_{\mathbf{N}} } \hat g_{\mathbf{k}} \,
        \eim{\mathbf{k} \left(\mathbf{n}^{-1}\odot \mathbf{l}\right)}.
    \end{equation*}

    \STATE For $j=0,\hdots,M-1$ compute
    \begin{equation*}
      f_j := \sum\limits_{\mathbf{l} \in I_{\mathbf{n},m} \left({\mathbf{x}}_j\right)
      } g_{\mathbf{l}} \, \tilde\psi\left({\mathbf{x}}_j - \mathbf{n}^{-1}\odot \mathbf{l}\right).
    \end{equation*}
  \end{algorithmic}
  Output: approximate values $f_j,\;j=0,\hdots,M-1$.\\
  Complexity: ${\cal O} (|\mathbf{N}| \log |\mathbf{N}|+M)$.
\end{algorithm}

Algorithm~\ref{algo:nfft} reads in matrix vector notation as
\begin{equation*}
  \mathbf{A} \mathbf{\hat{f}} \approx \mathbf{B} \mathbf{B} \mathbf{B} \mathbf{\hat{f}},
\end{equation*}
where $\mathbf{B}$ denotes the real $M \times \left|I_{\mathbf{n}}\right|$ sparse 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*}
where $\mathbf{B}$ is the Fourier matrix of size $\left|I_{\mathbf{n}}\right| \times
\left|I_{\mathbf{n}}\right|$, and where $\mathbf{B}$ is the real $\left|I_{\mathbf{n}}\right|
\times |I_{\mathbf{N}}|$ 'diagonal' matrix
\begin{equation*}
  \mathbf{B} := \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)^{\rm {T}}
\end{equation*}
with zero matrices $\mathbf{O}_t$ of size $N_t \times \frac{n_t-N_t}{2}$.

The corresponding computation of the adjoint matrix vector product reads as
\begin{equation*}
  \mathbf{A}^{\adj} \mathbf{\hat{f}} \approx \mathbf{B}^{\top}  \mathbf{B}^{\adj} \mathbf{B}^{\top}
  \mathbf{\hat{f}} \, .
\end{equation*}
With the help of the transposed index set
\begin{equation*}
  I_{\mathbf{n},m}^{\top}\left(\mathbf{l}\right) := \left\{j=0,\hdots,M-1 : \mathbf{l}-m
  \mathbf{1} \le \mathbf{n} \odot \mathbf{x}_j \le \mathbf{l}+m \mathbf{1}\right\},
\end{equation*}
one obtains Algorithm~\ref{algo:nfftH} for the adjoint NFFT.
\begin{algorithm}[htbp]
  \caption{\tt NFFT$^{\adj}$\label{algo:nfftH}}
  Input: $d,M \in \N$, $\mathbf{N}\in 2\N^d$\\
  {\invisible{Input:}} ${\mathbf{x}}_j\in [-\frac{1}{2},\frac{1}{2}]^d,\; j=0,\hdots,M-1$, and $f_j\in \C,\; j=0,\hdots,M-1$,\\[1ex]
                                %Input: $\mathbf{N}\in \N^d$, $\mathbf{\sigma} \in \mathbb{R}^d ,\;\sigma_t > 1,\;t=0,\hdots,d-1$, $M \in \N$, $\mathbf{m} \in \N^d$\\
                                %{\invisible{Input:}} ${\mathbf{x}}_j\in {\T}^d,\; j=0,\hdots,M-1$, and $f_j \in \C,\; j=0,\hdots,M-1$\\[1ex]
                                %Precomputation: $\mathbf{n} := \sigma \odot \mathbf{N}$,\\
                                %{\invisible{Precomputation:}} $c_{\mathbf{k}}\left(\tilde \varphi\right), \; \mathbf{k} \in I_{\mathbf{N}}$
                                %if the flag (see Section \ref{sect:parameters}) {\tt PRE\_PHI\_HUT} is set,\\
                                %{\invisible{Precomputation:}} $\tilde\psi\left({\mathbf{x}}_j - \mathbf{n}^{-1}\odot \mathbf{l}\right),
                                %\, j=0,\hdots,M-1 , \, \mathbf{l} \in I_{\mathbf{n},\mathbf{m}} \left({\mathbf{x}}_j\right)$ if {\tt PRE\_PSI} is set.
  \begin{algorithmic}[1]
    \topsep=1.5ex
    \itemsep=1.5ex

    \STATE For $\mathbf{l} \in I_{\mathbf{n}}$ compute
    \begin{equation*}
      g_{\mathbf{l}} := \sum\limits_{j \in I_{\mathbf{n},m}^{\top} (\mathbf{l}) } f_j \,
      \tilde\psi\left({\mathbf{x}}_j - \mathbf{n}^{-1}\odot \mathbf{l}\right)\, .
    \end{equation*}

    \STATE For $\mathbf{k} \in I_{\mathbf{N}}$ compute by $d$-variate (backward) {\rm FFT}
    \begin{equation*}
      \hat g_{\mathbf{k}} := \sum\limits_{\mathbf{l} \in I_{\mathbf{n}}}
      g_{\mathbf{l}} \, {\rm e}^{+2 \pi {\rm i} \mathbf{k} \left(\mathbf{n}^{-1} \odot
      \mathbf{l}\right)}.
    \end{equation*}

    \STATE For $\mathbf{k} \in I_{\mathbf{N}}$ compute
    \begin{equation*}
      \hat h_{\mathbf{k}} := \frac{\hat g_{\mathbf{k}}}{\left|I_{\mathbf{n}}\right|
      c_{\mathbf{k}}(\tilde \varphi)}.
    \end{equation*}
  \end{algorithmic}
  Output: approximate values $\hat h_{\mathbf{k}},\;\mathbf{k} \in I_{\mathbf{N}}$.\\
 Complexity: ${\cal O} (|\mathbf{N}| \log |\mathbf{N}|+M)$.
\end{algorithm}
Due to the characterisation of the nonzero elements of the matrix $\mathbf{B}$,
i.e.,
\begin{equation*}
  \bigcup\limits_{j=0}^{M-1} j \times I_{\mathbf{n},m} \left({\mathbf{x}}_j\right)=
  \bigcup\limits_{\mathbf{l} \in I_{\mathbf{n}}} I_{\mathbf{n},m}^{\top}\left(\mathbf{l}\right)
  \times \mathbf{l}.
\end{equation*}
the multiplication with the sparse matrix $\mathbf{B}^{\top}$ is implemented in a
'transposed' way in the library, summation as outer loop and only using the
multi-index sets $I_{\mathbf{n},m} \left({\mathbf{x}}_j\right)$.

%------------------------------------------------------------------------------
\subsection{Available window functions and evaluation techniques}
\label{sect:windows}
%------------------------------------------------------------------------------
Again, only the case $d=1$ is presented.
To keep the aliasing error and the truncation error small, several functions
$\varphi$ with good localisation in time and frequency domain were proposed,
e.g. the (dilated) {\it Gaussian} \cite{duro93, st97, DuSc}
\begin{eqnarray} \label{gaussian}
\varphi\left(x\right)
&=&
\left(\pi b\right)^{-1/2} \,
{\rm e}^{-\frac{\left(nx\right)^2}{b}}
\qquad
\left(b := \frac{2\sigma}{2 \sigma -1} \frac{m}{\pi} \right),\\
\hat \varphi \left(k\right)
&=&
\frac{1}{n} \, {\rm e}^{-b\left(\frac{\pi k}{n}\right)^2} , \nonumber
\end{eqnarray}
(dilated) {\it cardinal central $B$--splines} \cite{bey95, st97}
\begin{eqnarray} \label{bspline}
\varphi\left(x\right)
&=&
M_{2m}\left(n x\right),\\ \nonumber
\hat \varphi \left(k\right)
&=&
\frac{1}{n} \sinc^{2m} \left( k \pi/n\right),
\end{eqnarray}
where $M_{2m}$ denotes the centred cardinal $B$--Spline of order $2m$,\\
(dilated) {\it Sinc functions} \cite{pottshabil}
\begin{eqnarray} \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)
\nonumber
\end{eqnarray}
and
(dilated) {\it Kaiser--Bessel functions} \cite{Ja, four}
\begin{eqnarray} \label{kaiser}
\varphi\left(x\right)
&=& \frac{1}{\pi} \left\{
\ds  \begin{array}{ll}
\ds \frac{\sinh\left(b
\sqrt{m^2 - n^2 x^2}\right)}{\sqrt{m^2 - n^2 x^2}}
& \mbox{for } |x| \le \frac{m}{n}
\qquad
\left(b :=  \pi \left(2 - \frac{1}{\sigma} \right) \right) ,\\[4ex]
\ds \frac{\sin\left(b
\sqrt{n^2 x^2 - m^2}\right)}{\sqrt{n^2 x^2-m^2}}
& \mbox{otherwise},
 \end{array}
\right.
\\[1ex]
\hat \varphi \left(k\right)
&=& \frac{1}{n} \left\{
\ds  \begin{array}{cl}
 \, I_0 \left(m \sqrt{b^2 - \left(2 \pi k/n \right)^2} \right)
&
{\rm for} \; \; k= -n \left(1-\frac{1}{2\sigma} \right), \ldots,
n \left(1-\frac{1}{2\sigma}\right),\\
0&{\rm otherwise},
 \end{array}
\right.
\nonumber
\end{eqnarray}
where $I_0$ denotes the {\it modified zero--order Bessel function}.
For these functions $\varphi$ it has been proven that
$$
|f\left(\mathbf{x}_j\right) - s\left(\mathbf{x}_j\right)| \le C\left(\sigma,m\right) \|\mathbf{\hat{f}}\|_1
$$
where
$$
C\left(\sigma,m\right) :=
\left\{
\begin{array}{ll}
4 \, \e^{-m \pi \left(1-1/\left(2\sigma-1\right)\right)}& {\rm for} \; \left(\ref{gaussian}\right) ,\\
4 \, \left(\frac{1}{2\sigma -1} \right)^{2m}& {\rm for} \; \left(\ref{bspline}\right), \\[1ex]
\frac{1}{m-1}\left(\frac{2}{\sigma^{2m}}+\left(\frac{\sigma}{2\sigma-1}\right)^{2m}\right)&
{\rm for} \; \left(\ref{sinc}\right), \\[1ex]
4\pi\left(\sqrt{m}+m\right)\sqrt[4]{1-\frac{1}{\sigma}} \, \e^{-2 \pi m \sqrt{1-1/\sigma}} &  {\rm for} \; \left(\ref{kaiser}\right).
\end{array}
\right.
$$
Thus, for fixed $\sigma >1$, the approximation error introduced by the NFFT
decays exponentially with the number $m$ of summands in (\ref{approx}).
Using the tensor product approach the above error estimates can be generalised
for the multivariate setting \cite{elst}.
On the other hand, the complexity of the NFFT increases with $m$.

In the following, we suggest different methods for the compressed storage and
application of the matrix $\mathbf{B}$ which are all available within our NFFT
library by choosing particular flags in a simple way during the initialisation
phase.
These methods do not yield a different asymptotic performance but rather yield
a lower constant in the amount of computation.

\subsubsection*{Fully precomputed window function}
One 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)$
explicitly.
Thus, one has to store the large amount of $(2m+1)^d M$ real numbers but uses
no extra floating point operations during the matrix vector multiplication
beside the necessary $(2m+1)^d M$ flops.
Furthermore, we store for this method explicitly the row and column for each
nonzero entry of the matrix $\mathbf{B}$.
This method, included by the flag {\tt PRE\_FULL\_PSI}, is the fastest procedure
but can only be used if enough main memory is available.

\subsubsection*{Tensor product based precomputation}
Using the fact that the window functions are built as tensor products one can
store $\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 uses a medium amount of memory to store $d(2m+1)M$ real numbers
in total.
However, one has to carry out for each node at most $2(2m+1)^d$ extra
multiplications to compute from the factors 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)$.
%where a naive version would take $d(2m+1)^d$ multiplications.
Note, that this technique is available for every window function discussed
here and can be used by means of the flag {\tt PRE\_PSI} which is also the
default method within our software library.

\subsubsection*{Linear interpolation from a lookup table}
For a large number of nodes $M$, the amount of memory can by further reduced
by the use of 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$
and then compute 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)$ by means of the linear interpolation from its two
neighbouring precomputed samples.

This method needs only a storage of $dK$ real numbers in total where $K$
depends solely on the target accuracy but neither on the number of nodes $M$
nor on the multidegree $\mathbf{N}$.
Choosing $K$ to be a multiple of $m$, we further reduce the computational
costs during the interpolation since the distance from $(\mathbf{x}_j)_t-
\frac{l_t}{n_t}$ to the two neighbouring interpolation nodes and hence the
interpolation weights remain the same for all $l_t\in I_{n_t,m}((
\mathbf{x}_j)_t)$.
This method requires $2(2m+1)^d$ extra multiplications per node and is
used within the NFFT by the flag {\tt PRE\_LIN\_PSI}. Error estimates for this approximation are given in \cite{KuPo06}.

\subsubsection*{Fast Gaussian gridding}
Two useful properties of the Gaussian window function \eqref{gaussian} within
the present framework were recently reviewed in \cite{GrLe04}.
Beside its tensor product structure for $d>1$, which also holds for all other
window functions, it is remarkable that the number of evaluations of the form
{\tt exp()} can be greatly decreased.
More precisely, for $d=1$ and a fixed node $x_j$ the evaluations of
$\varphi(x_j-\frac{l'}{n})$, $l'\in I_{n,m}(x_j)$, can be reduced by the
splitting
\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$.
Note, that the first factor and the exponential within the brackets are
constant for each fixed node $x_j$.
Once, we evaluate the second exponential, its $l$-th power can be computed
consecutively by multiplications only.
Furthermore, the last exponential is independent of $x_j$ and these $2m+1$
values are computed only once within the NFFT and their amount is negligible.
Thus, it is sufficient to store or evaluate $2M$ exponentials for $d=1$.
The case $d>1$ uses $2dM$ storages or evaluations by using the general tensor
product structure.
This method is employed by the flags {\tt FG\_PSI} and {\tt PRE\_FG\_PSI} for
the evaluation or storage of $2d$ exponentials per node, respectively.

\subsubsection*{No precomputation of the window function}
The last considered method uses no precomputation at all, but rather evaluates
the univariate window function $(2m+1)^d M$ times.
Thus, the computational time depends on how fast we can evaluate the
particular window function.
However, no additional storage is necessary which suits this approach whenever
the problem size reaches the memory limits of the used computer.

%------------------------------------------------------------------------------
\subsection{Further NFFT approaches}
%------------------------------------------------------------------------------
Several papers have described fast approximations for the NFFT.
Common names for NFFT are
{\it non-uniform fast Fourier transform} \cite{fesu02},
{\it generalised fast Fourier transform} \cite{duro93},
{\it unequally-spaced fast Fourier transform} \cite{bey95},
{\it fast approximate Fourier transforms for irregularly spaced data} \cite{ware98},
{\it non-equispaced fast Fourier transform} \cite{Fou02} or
{\it gridding} \cite{scsc, Ja, Pe}.

In various papers, different window functions were considered, e.g. Gaussian
pulse tapered with a Hanning window in \cite{DuSc}, Gaussian kernels combined
with Sinc kernels in \cite{Pe}, and special optimised windows in
\cite{Ja,DuSc}.

A simple but nevertheless fast scheme for the computation of \eqref{ndft} in
the univariate case $d=1$ is presented in \cite{AnDa96}.
This approach uses for each node $x_j \in [-\frac{1}{2},\frac{1}{2})$ a $m$-th
order Taylor expansion of the trigonometric polynomial in \eqref{trigPoly}
about the nearest neighbouring point on the oversampled equispaced lattice
$\{n^{-1} k - \frac{1}{2}\}_{k=0,\hdots,n-1}$ where again $n=\sigma N,\,
\sigma\gg 1$.
Besides its simple structure and only $\cO(N\log N+M)$ arithmetic operations,
this algorithm utilises $m$ FFTs of size $n$ compared to only one in the NFFT
approach, uses a medium amount of extra memory, and is not suited for highly
accurate computations, see \cite{KuPo06}.
Furthermore, its extension to higher dimensions has not been considered so
far.

Another approach for the univariate case $d=1$ is considered in \cite{DuRo95}
and based on a Lagrange interpolation technique.
After taking a $N$-point FFT of the vector $\mathbf{\hat{f}}$ in \eqref{ndft} one
uses an exact polynomial interpolation scheme to obtain the values of the
trigonometric polynomial $f$ at the nonequispaced nodes $x_j$.
Here, the time consuming part is the exact polynomial interpolation scheme
which can however be realised fast in an approximate way by means of the fast
multipole method.
This approach is appealing since it allows also for the inverse transform.
Nevertheless, numerical experiments in \cite{DuRo95} showed that this
approach is far more time consuming than Algorithm \ref{algo:nfft} and the
inversion can only be computed in a stable way for almost equispaced nodes
\cite{DuRo95}.

Furthermore, special approaches based on scaling vectors \cite{NgLi99}, based
on minimising the Frobenius norm of certain error matrices \cite{st01} or
based on min-max interpolation \cite{fesu02} are proposed.
While these approaches gain some accuracy for the Gaussian or B-Spline
windows, no reasonable improvement is obtained for the Kaiser-Bessel window
function.

For comparison of different approaches, we refer to
\cite{ware98,st01,fesu02,KuPo06}.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\clearpage
\section{Generalisations and inversion} \label{sect:gen}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

We consider generalisations of the NFFT for
%------------------------------------------------------------------------------
\subsection{NNFFT - nonequispaced in time and frequency fast Fourier
  transform}
%------------------------------------------------------------------------------

Now we are interested in the computation of
\begin{equation*}
  f_j= \sum_{k=0}^{N-1} \hat f_k \eim{\left(\mathbf{v}_k \odot \mathbf{N}\right))\mathbf{x}_j}
\end{equation*}
for $j =0,\hdots,M-1$, $\mathbf{v}_k, \mathbf{x}_j \in \T^d$, and $\hat f_k \in \C$, where $\mathbf{N}\in \N^d$
denotes the "nonharmonic" bandwidth.
The corresponding fast algorithm is known as {\it fast Fourier transform for
nonequispaced data in space {\bf and} frequency domain} (NNFFT) \cite{duro93,elst,postta01} or as type 3 nonuniform FFT \cite{LeGr05}.

A straightforward evaluation of this sum with the standard NFFT is not
possible, since the samples in neither domain are equispaced.
However, an so-called NNFFT was first suggested in \cite{duro93} and later
studied in more depth in \cite{elst}, which permits the fast calculation of
the Fourier transform of a vector of nonequispaced samples at a vector of
nonequispaced positions.
It constitutes a combination of the standard NFFT and its adjoint see also
\cite{postta01}.

%------------------------------------------------------------------------------
\subsection{NFCT/NFST - nonequispaced fast (co)sine transform}
%------------------------------------------------------------------------------

Let nonequispaced nodes $\mathbf{x}_j \in [0,\frac{1}{2}]^d$ and frequencies
$\mathbf{k}$ in 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*}
be given.
For notational convenience let furthermore $\cos(\mathbf{k} \odot \mathbf{x}) :=
\cos(k_0 x_0) \cdot \hdots \cdot \cos(k_{d-1} x_{d-1})$ and $\sin(\mathbf{k} \odot
\mathbf{x}) := \sin(k_0 x_0) \cdot \hdots \cdot \sin(k_{d-1} x_{d-1})$.

The {\it nonequispaced discrete cosine} and {\it sine transforms} are given by
\begin{align*}
  f\left(\mathbf{x}_j\right) &= \sum_{\mathbf{k}\in I_{\mathbf{N}}^C} \hat{f}_{\mathbf{k}}
  \cos(2\pi(\mathbf{k}\odot\mathbf{x}_j)),\\
  f\left(\mathbf{x}_j\right) &= \sum_{\mathbf{k}\in I_{\mathbf{N}}^S} \hat{f}_{\mathbf{k}}
  \sin(2\pi(\mathbf{k}\odot\mathbf{x}_j)),
\end{align*}
for $j=0,\hdots,M-1$ and real coefficients $\hat f_{\mathbf{k}} \in \R$,
respectively.

The straight forward algorithm of this matrix vector product, which is called
ndct and ndst, takes $\cO(M |I_{\mathbf{N}}^C|)$ and $\cO(M |I_{\mathbf{N}}^C|)$
arithmetical operations.
For these real transforms the adjoint transforms coincide with the ordinary
transposed matrix vector products.
Our fast approach is based on the NFFT and seems to be easier than the
Chebyshev transform based derivation in \cite{po01} and faster than the
algorithms in \cite{tili00} which still use FFTs.
Instead of FFTs we use fast algorithms for the discrete cosine transform
(DCT--I) and for the discrete sine transform (DST--I).
For details we refer to \cite{fepo02}.

%------------------------------------------------------------------------------
\subsection{NSFFT - nonequispaced sparse fast Fourier transform}
%------------------------------------------------------------------------------

We consider the fast evaluation of trigonometric polynomials from so-called
hyperbolic crosses.
In multivariate approximation one has to deal with the so called `curse of
dimensionality', i.e., the number of degrees of freedom for representing an
approximation of a function with a prescribed accuracy depends exponentially
on the dimensionality of the considered problem.
This obstacle can be circumvented to some extend by the interpolation on
sparse grids and the related approximation on hyperbolic cross points in the
Fourier domain, see, e.g., \cite{Zeng91,Spr00,BuZe04}.

Instead of approximating a Fourier series on the standard tensor product grid
$I_{\mathbf{(N,\hdots,N)}}$ with $\cO(N^d)$ degrees of freedom, it can be
approximated with only $\cO(N \log^{d-1} N)$ degrees of freedom from the
hyperbolic cross
\begin{equation*}
  H_N^d := \bigcup_{\mathbf{N}\in\N^d,\, |I_{\mathbf{N}}|=N} I_{\mathbf{N}},
\end{equation*}
where $N=2^{J+2},\;J\in\NZ$ and we allow $N_t=1$ in the $t$-th coordinate in
the definition of $I_{\mathbf{N}}$.
The approximation error in a suitable norm (dominated mixed smoothness) can be shown to deteriorate only by a factor of $\log^{d-1} N$,
cf. \cite{Spr00}.
\begin{figure}[ht]
  \centering
  \begin{tabular}{cccc}
    \includegraphics[scale=0.6]{images/sparse2D_0}&
    \includegraphics[scale=0.6]{images/sparse2D_1}&
    \includegraphics[scale=0.6]{images/sparse2D_2}&
    \includegraphics[scale=0.6]{images/sparse2D_3}
  \end{tabular}
  \caption{Hyperbolic cross points for $d=2$ and $J=0,\ldots,3$.
    \label{sparse_grid}}
\end{figure}

The {\em nonequispaced sparse discrete Fourier transform} (NSDFT) is the
evaluation of
\begin{equation*}
  f\left(\mathbf{x}_j\right) = \sum_{\mathbf{k}\in H_N^d}
  \hat f_{\mathbf{k}} \eim{\mathbf{k} \mathbf{x}_j}
\end{equation*}
for given Fourier coefficients $\hat f_{\mathbf{k}}\in\C$ and nodes $\mathbf{x}_j\in
\T^d$.
The number of used arithmetical operations is $\cO(M N\log^{d-1}N)$.
This is reduced by our fast schemes to $\cO(N\log^2 N+M)$ for $d=2$ and
$\cO(N^{3/2} \log N+M)$ for $d=3$, see \cite{FeKuPo} for details.

%------------------------------------------------------------------------------
\subsection{FPT - fast polynomial transform} \label{sect:fpt}
%------------------------------------------------------------------------------
A \emph{discrete polynomial transform} (DPT) is a generalisation of the DFT
from the basis of complex exponentials $\mathrm{e}^{\mathrm{i} k x}$ to an
arbitrary systems of algebraic polynomials satisfying a three-term recurrence
relation.
More precisely, let $P_{0},P_{1},\hdots:[-1,1]\rightarrow\R$ be a sequence of
polynomials that satisfies a three-term recurrence relation
\begin{equation*}
  P_{k+1}(x) = (\alpha_{k} x  + \beta_{k})P_{k}(x) + \gamma_{k} P_{k-1}(x),
\end{equation*}
with $\alpha_{k} \neq 0$, $k\ge 0$, $P_{0}(x) := 1$, and $P_{-1}(x) := 0$.
Clearly, every $P_{k}$ is a polynomial of exact degree $k$ and typical
examples are the classical orthogonal Jacobi polynomials
$P_{k}^{(\alpha,\beta)}$.

Now, let $f:[-1,1]\rightarrow\R$ be a polynomial of degree $N\in\N$ given by
the finite linear combination
\begin{equation*}
  f(x) = \sum_{k=0}^{N} a_{k} P_{k}(x).
\end{equation*}
The discrete polynomial transform (DPT) and its fast version, the \emph{fast
  polynomial transform} (FPT), are the transformation of the coefficients
$a_{k}$ into coefficients $b_{k}$ from the orthogonal expansion of $f$ into
the basis of Chebyshev polynomials of the first kind $T_{k}(x) := \cos(k
  \arccos x)$, i.e.,
\begin{equation*}
  f(x) = \sum_{k=0}^N b_{k} T_{k}(x).
\end{equation*}
A direct algorithm for computing this transformation needs $\cO(N^2)$
arithmetic operations.
The FPT algorithm implemented in the NFFT library follows the approach in
\cite{postta98} and is based on the idea of using the three-term-recurrence
relation repeatedly.
Together with a method for fast polynomial multiplication in the Chebyshev
basis and a cascade-like summation process, this yields a method for computing
the polynomial transform with $\cO(N \log^2 N)$ arithmetic operations.
For more detailed information, we refer the reader to
\cite{drhe, drhero, postta98, po01, KePo06} and the references therein.

%------------------------------------------------------------------------------
\subsection{NFSFT - nonequispaced fast spherical Fourier transform}
\label{sect:gen:nfsft}
%------------------------------------------------------------------------------
Fourier analysis on the sphere has, despite other fields, practical
relevance in tomography, geophysics, seismology, meteorology and
crystallography. In analogy to the complex exponentials
$\mathrm{e}^{\mathrm{i} k x}$ on the torus, the spherical harmonics
form the orthogonal Fourier basis with respect to the usual inner product on  the sphere.

\subsubsection*{Spherical coordinates}
Every point in $\R^3$ can be described in spherical coordinates by a vector
$(r,\vartheta,\varphi)^{\top}$ with the radius $r\ge 0$ and two angles
$\vartheta \in [0,\pi]$, $\varphi \in [0,2\pi)$.
We denote by $\S^2$ the two-dimensional unit sphere embedded into $\R^3$,
i.e.
\begin{equation*}
  \S^2 := \left\{\mathbf{x} \in \R^{3}:\; \|\mathbf{x}\|_2=1\right\}
\end{equation*}
and identify a point from $\S^2$ with the corresponding vector
$(\vartheta,\varphi)^{\top}$.
The spherical coordinate system is illustrated in Figure \ref{sphere}.
\begin{figure}[tb]
  \centering
  \includegraphics[width=0.60\textwidth]{images/sphere}
  \caption{The spherical coordinate system in $\R^3$: Every point $\boldsymbol{\xi}$ on a
    sphere with radius $r$ centred at the origin can be described by angles
  $\vartheta \in [0,\pi]$, $\varphi \in [0,2\pi)$ and the radius $r \in \mathbb{R}^+$.
  For $\vartheta = 0$ or $\vartheta = \pi$ the point $\boldsymbol{\xi}$ coincides with the
  North or the South pole, respectively.} \label{sphere}
\end{figure}

\subsubsection*{Legendre polynomials and associated Legendre functions}
The Legendre polynomials $P_k : [-1,1]\rightarrow \R$, $k\ge 0$, as
classical orthogonal polynomials are given by their corresponding Rodrigues
formula
\begin{equation*}
  P_k(x) := \frac{1}{2^k k!} \frac{\text{d}^k}{\text{d} x^k}
  \left(x^2-1\right)^k.
\end{equation*}

The associated Legendre functions $P_k^n : [-1,1] \rightarrow \R$, $k \ge n
\ge 0$ are defined by
\begin{equation*}
  P_k^n(x) := \left(\frac{(k-n)!}{(k+n)!}\right)^{1/2}
  \left(1-x^2\right)^{n/2} \frac{\text{d}^n}{\text{d} x^n} P_k(x).
\end{equation*}
For $n = 0$, they coincide with the Legendre polynomials $P_k = P_k^0$.
The associated Legendre functions $P_{k}^n$ obey the three-term recurrence
relation
\begin{equation*}
  P_{k+1}^n(x) = \frac{2k+1}{((k-n+1)(k+n+1))^{1/2}} x P_k^n(x) -
  \frac{((k-n)(k+n))^{1/2}}{((k-n+1)(k+n+1))^{1/2}} P_{k-1}^n(x)
\end{equation*}
for $k \ge n \ge 0$, $P_{n-1}^n(x) = 0$, $P_{n}^n(x) = \frac{\sqrt{(2n)!}}{2^n
  n!}\left(1-x^2\right)^{n/2}$.
For fixed $n$, the set $\{P_k^n:\: k \ge n\}$ forms a set of
orthogonal functions, i.e.,
\begin{equation*}
  \left\langle P_k^n,P_l^n \right\rangle =
  \int_{-1}^{1}  P_k^n(x) P_l^n(x) \text{d} x =
  \frac{2}{2k+1} \delta_{k,l}.
\end{equation*}
Again, we denote by $\bar{P}_{k}^n = \sqrt{\frac{2k+1}{2}} P_k^n$ the
orthonormal associated Legendre functions.
In the following, we allow also for $n < 0$ and set $P_{k}^n$ := $P_{k}^{-n}$
in this case.

\subsubsection*{Spherical harmonics}
The spherical harmonics $Y_k^n : \S^2 \rightarrow \C$, $k \ge |n|$, $n\in\Z$,
are given by
\begin{equation*}
  Y_k^n(\vartheta,\varphi) := P_k^{n}(\cos\vartheta)
  \, \mathrm{e}^{\mathrm{i} n \varphi}.
\end{equation*}
They form an orthogonal basis for the space of square integrable functions on
the sphere, i.e.,
\begin{equation*}
  \left\langle Y_k^n,Y_l^m \right\rangle  =
  \int_{0}^{2\pi} \int_{0}^{\pi} Y_k^n(\vartheta,\varphi)
  \overline{Y_l^m(\vartheta,\varphi)} \sin \vartheta \; \mathrm{d} \vartheta
  \; \mathrm{d} \varphi =
  \frac{4\pi}{2k+1} \delta_{k,l}\delta_{n,m}.
\end{equation*}
The orthonormal spherical harmonics are denoted by $\bar{Y}_k^n =
\sqrt{\frac{2k+1}{4\pi}} Y_k^n$.

Hence, any square integrable function $f:\S^2\rightarrow\C$ has the expansion
\begin{equation*}
  f = \sum_{k=0}^{\infty} \sum_{n=-k}^{k} \hat{f}(k,n) \bar{Y}_k^n,
\end{equation*}
with the spherical Fourier coefficients $\hat{f}(k,n) = \left\langle f,
\bar{Y}_k^n \right\rangle$.
The function $f$ is called \emph{bandlimited}, if $\hat{f}(k,n)=0$ for $k>N$
and some $N\in\N$.
In analogy to the NFFT on the Torus $\mathbb{T}$, we define the sampling set
\begin{equation*}
  \mathcal{X}: = \left\{\left(\vartheta_{j},\varphi_{j}\right) \in
  \mathbb{S}^2:\,j=0,\hdots,M-1\right\}
\end{equation*}
of nodes on the sphere $\S^2$ and the set of possible ''frequencies''
\begin{equation*}
  \mathcal{I}_{N} := \left\{(k,n):\,k=0,1,\ldots,N;\;
  n=-k,\ldots,k\right\}.
\end{equation*}
The \emph{nonequispaced discrete spherical Fourier transform} (NDSFT) is
defined as the evaluation of the finite spherical Fourier sum
\begin{equation}
  \label{eq.ndsft}
   f_j=
  f\left(\vartheta_{j},\varphi_{j}\right) =
    \sum_{(k,n) \in \, \mathcal{I}_{N}} \hat{f}_k^n \,
      Y_k^n\left(\vartheta_{j},\varphi_{j}\right) =
    \sum_{k=0}^N \sum_{n=-k}^k \hat{f}_k^n \,
      Y_k^n\left(\vartheta_{j},\varphi_{j}\right)
\end{equation}
for $j=0,\hdots,M-1$.
In matrix vector notation, this reads $\mathbf{f} = \mathbf{Y} \; \mathbf{\hat
  f}$ with
\begin{equation*}
  \begin{split}
    \mathbf{f} & := \left(f_j\right)_{j=0}^{M-1} \in
      \mathbb{C}^{M},\ f_{j} := f\left(\vartheta_{j},\varphi_{j}\right),\\
    \mathbf{Y} & :=
      \left(Y_k^n\left(\vartheta_j,\varphi_j\right)\right)_{j=0,\ldots,M-1;\;
      (k,n) \in \mathcal{I}_N} \in \mathbb{C}^{M \times (N+1)^2},\\
    \mathbf{\hat f} & := \left(\hat f_k^n\right)_{(k,n) \in \mathcal{I}_N} \in
      \mathbb{C}^{(N+1)^2}.
  \end{split}
\end{equation*}
The corresponding \emph{adjoint nonequispaced
  discrete fast spherical Fourier transform} (adjoint NDSFT) is defined as
the evaluation of
\begin{equation*}
  \label{eq.adjointndsft}
  \hat h_{k}^n =
  \sum_{j = 0}^{M-1} f\left(\vartheta_{j},\varphi_{j}\right)
  \overline{Y_k^n\left(\vartheta_{j},\varphi_{j}\right)}
\end{equation*}
for all $(k,n) \in \mathcal{I}_{N}$.
Again, in matrix vector notation, this reads $\mathbf{\hat h} =
\mathbf{Y}^{\adj} \, \mathbf{f}$.

The coefficients $\hat h_{k}^n$ are, in general, not
identical to the Fourier coefficients $\hat{f}(k,n)$ of the function $f$.
However, provided that for the sampling set $\mathcal{X}$ a quadrature rule
with weights $w_{j}$, $j=0,\hdots,M-1$, and sufficient degree of exactness
is available, one might recover the Fourier coefficients $\hat{f}(k,n)$ by
evaluating
\begin{equation*}
  \hat{f}_{k}^n = \int_{0}^{2\pi} \int_{0}^{\pi} f(\vartheta,\varphi)
  \overline{Y_{k}^n(\vartheta,\varphi)} \sin\vartheta \; \mathrm{d} \vartheta
  \; \mathrm{d} \varphi =
  \sum_{j = 0}^{M-1} w_{j} f\left(\vartheta_{j},\varphi_{j}\right)
  \overline{Y_k^n\left(\vartheta_{j},\varphi_{j}\right)}
\end{equation*}
for all $(k,n) \in \mathcal{I}_{N}$.

Direct algorithms for computing the NDSFT and adjoint NDSFT transformations
need $\cO(M N^2)$ arithmetic operations.
A combination of the fast polynomial transform and the NFFT leads to
approximate algorithms with $\cO(N^2 \log^2 N + M)$ arithmetic operations.
These are denoted \emph{NFSFT} and \emph{adjoint NFSFT}, respectively.
The NFSFT algorithm using  the FPT and the NFFT was introduced in
\cite{kupo02} while the adjoint NFSFT variant was developed in \cite{KePo06}.

\begin{algorithm}[tb]
  \caption[NFSFT]{Nonequispaced fast spherical Fourier transform (NFSFT)}
  \label{NFSFT:Algorithm:NFSFT}
  \begin{algorithmic}
    \STATE  Input: $N \in \mathbb{N}_0$, $M \in \mathbb{N}$, spherical Fourier
      coefficients $\mathbf{\hat{f}} = (\hat{f}_{k}^n)_{(k,n) \in
      \mathcal{I}_{N}} \in \mathbb{C}^{(N+1)^2}$,\\
    \STATE \invisible{Input:} a sampling set $\mathcal{X} = \left(\vartheta_{j},
      \varphi_{j}\right)_{j=0}^{M-1} \in ([0,\pi]\times [0,2\pi))^{M}$.
    \STATE
    \FOR {$n=-N,\ldots,N$}
      \STATE Compute the Chebyshev coefficients $\left(b_{k}^n\right)_{k=0}^{N}$
        of $g^{n}$ by a fast polynomial transform.
      \STATE Compute the coefficients $\left(c_{k}^n\right)_{k=-N}^N$
        %as in \eqref{NFSFT:afterChebToExp}
        from the coefficients $\left(b_{k}^n\right)_{k=0}^{M}$.
    \ENDFOR
    \STATE Compute the function values
      $\left(f\left(\vartheta_{j},\varphi_{j}\right)\right)_{j=0}^{M-1}$
      by evaluating the Fourier sum %in \eqref{NFSFT:transformed}
      using a fast two-dimensional NFFT.
    \STATE
    \STATE Output: The function values ${\mathbf{f}} =
      \left(f\left(\vartheta_{j},
      \varphi_{j}\right)\right)_{j=0}^{M-1} \in \mathbb{C}^{M}$.
    \STATE Complexity: $\mathcal{O}\left(N^2 \log^2 N +
      M\right)$.
\end{algorithmic}
\end{algorithm}

\begin{algorithm}[tb]
  \caption[Adjoint NFSFT]{Adjoint nonequispaced fast spherical Fourier
    transform (adjoint NFSFT)}
  \label{NFSFT:Algorithm:adjointNFSFT}
  \begin{algorithmic}
    \STATE  Input: $N \in \mathbb{N}_0$, $M \in \mathbb{N}$, a sampling set
      $\mathcal{X} = \left(\vartheta_{j}, \varphi_{j}\right)_{j=0}^{M-1}
      \in ([0,\pi]\times [0,2\pi))^{M}$,\\
    \STATE \invisible{Input:} values $\mathbf{\tilde{f}} =
      (\tilde{f}_{j})_{j=0}^{M-1} \in \mathbb{C}^M$.
    \STATE
    \STATE Compute the coefficients $\left(\tilde{c}_{k}^n\right)_{k,n=-N}^N$ from
      the values $\left(f_{j}\right)_{j=0}^{N-1}$.
      %according to the matrix $\mathbf{F}_{\mathcal{X}}^{\h}$ in
      %\eqref{NFSFT:adjointNFSFTmv}
    \FOR {$n=-N,\ldots,N$}
      \STATE Compute the coefficients $\left(\tilde{b}_{k}^n\right)_{k=0}^N$ from the
        coefficients $\left(\tilde{c}_{k}^n\right)_{k=-N}^{N}$.
        % according to the matrix $\mathbf{B}^{\transp}$ in \eqref{NFSFT:adjointNFSFTmv}
      \STATE Compute the coefficients $\left(\tilde{a}_{k}^n\right)_{k=|n|}^{N}$ from
        the coefficients $\left(\tilde{b}_{k}^n\right)_{k=0}^N$
        %according to the matrix $\mathbf{T}^{\transp}$ in \eqref{NFSFT:adjointNFSFTmv}
        by a fast transposed polynomial transform.
    \ENDFOR
    \STATE
    \STATE Output: Coefficients $\mathbf{h} =
      \left(\tilde{h}_{k}^n\right)_{(k,n)
      \in \mathcal{I}_N} \in \mathbb{C}^{(N+1)^2}$.
    \STATE Complexity: $\mathcal{O}\left(N^2 \log^2 N +
      M\right)$.
\end{algorithmic}
\end{algorithm}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Solver - inverse transforms} \label{sect:inv}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


In the following, we describe the inversion of the NFFT, i.e., the computation
of Fourier coefficients from given samples $(\mathbf{x}_j,y_j) \in \T^d\times \C$,
$j=0,\ldots,M-1$.
In matrix vector notation, we aim to solve the linear system of equations
\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}|}}$.
Note however that the nonequispaced Fourier matrix $\mathbf{A}$ can be replaced by
any other Fourier matrix from Section \ref{sect:gen}.

Typically, the number of samples $M$ and the dimension of the space of the
polynomials $|I_{\mathbf{N}}|$ do not coincide, i.e., the matrix $\mathbf{A}$ is
rectangular.
There are no simple inverses to nonequispaced Fourier matrices in general and
we search for some pseudoinverse solution $\mathbf{\hat{f}}^{\dagger}$, see
e.g. \cite[p. 15]{Bj96}.
However, we conclude from eigenvalue estimates in
\cite{FeGrSt95,BaGr03,KuPo04} that the matrix $\mathbf{A}$ has full rank if
\begin{equation*}
  \max_{0\le t<d} N_t < c_d \delta^{-1}\,, \qquad \delta:=2\max_{
  \mathbf{x}\in\T^d}\min_{j=0,\hdots,M-1}\dist_{\infty}\left(\mathbf{x}_j,\mathbf{x}_l\right),
\end{equation*}
or
\begin{equation*}
  \min_{0\le t<d} N_t>C_d q^{-1}\,, \qquad q:=\min_{0\le j<l<M}
  \dist_{\infty}\left(\mathbf{x}_j,\mathbf{x}_l\right).
\end{equation*}

In what follows, we consider a weighted least squares solution for the
over-determined case and an optimal interpolation for the consistent
underdetermined case.
Both cases are solved by means of iterative algorithms where Algorithm
\ref{algo:nfft} and \ref{algo:nfftH} are used for the matrix vector
multiplication with $\mathbf{A}$ and $\mathbf{A}^{\adj}$, respectively.

%------------------------------------------------------------------------------
\subsubsection*{Least squares}
%------------------------------------------------------------------------------

For a comparative low polynomial degree $|I_{\mathbf{N}}|<M$ the system
\eqref{eq:infft} is over-determined, so that in general the given data $y_j\in
\C,\;j=0,\hdots,M-1,$ will be only approximated up to the {\em residual} $
\mathbf{r}:=\mathbf{f}-\mathbf{A} \mathbf{\hat{f}}$.
One considers the {\em weighted approximation problem}
\begin{equation*}
  \|\mathbf{y} - \mathbf{A} \mathbf{\hat{f}}\|_{\mathbf{W}}
  = \left(\sum_{j=0}^{M-1} w_j |y_j-f(x_j)|^2\right)^{1/2}
  \stackrel{\mathbf{\hat{f}}}{\rightarrow} \min,
\end{equation*}
which may incorporate weights $w_j> 0$, $\mathbf{W}:=\diag(w_j)_{j=0,\hdots,M-1}$,
to compensate for clusters in the sampling set ${\cal X}$.
This least squares problem is equivalent to the
{\em weighted normal equation of first kind}
\begin{equation}\label{eq:no1}
  \mathbf{A}^{\adj} \mathbf{W} \mathbf{A} \mathbf{\hat{f}} = \mathbf{A}^{\adj} \mathbf{W} \mathbf{y}.
\end{equation}

Applying the Landweber (also known as Richardson, ...), the steepest descent,
or the conjugate gradient scheme to \eqref{eq:no1} yields the following
Algorithms \ref{algo:landweber}-\ref{algo:wCGNR(E)}.
\begin{algorithm}[ht!]
  \caption{Landweber}
  \label{algo:landweber}
  Input: $\mathbf{y} \in \mathbb{C}^M,\, {\mathbf{\hat{f}}}_0 \in \mathbb{C}^{|I_{
  \mathbf{N}}|},\, \alpha>0$
  \begin{algorithmic}[1]
    \topsep=1.5ex
    \itemsep=1.5ex
    \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$

    \FOR{$l=0,\hdots$}
    \topsep=1.5ex
    \itemsep=1.5ex
    \STATE $\mathbf{\hat{f}}_{l+1}=\mathbf{\hat{f}}_{l}
    +\alpha \mathbf{\hat{W}} \mathbf{\hat{z}}_l$
    \STATE $\mathbf{r}_{l+1}=\mathbf{y} - \mathbf{A} \mathbf{\hat{f}}_{l+1}$
    \STATE $\mathbf{\hat{z}}_{l+1}=\mathbf{A}^{\adj} \mathbf{W} \mathbf{r}_{l+1}$
    \ENDFOR
  \end{algorithmic}
  Output: $\mathbf{\hat{f}}_{l}$
\end{algorithm}

\begin{algorithm}[ht!]
  \caption{Steepest descent}
  \label{algo:gradient}
  Input: $\mathbf{y} \in \mathbb{C}^M,\, {\mathbf{\hat{f}}}_0 \in \mathbb{C}^{|I_{\mathbf{N}}|}$
  \begin{algorithmic}[1]
    \topsep=1.5ex
    \itemsep=1.5ex
    \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$

    \FOR{$l=0,\hdots$}
    \topsep=1.5ex
    \itemsep=1.5ex
    \STATE $\mathbf{v}_l = \mathbf{A} \mathbf{\hat{W}} \mathbf{\hat{z}}_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{z}}_l$
              \STATE $\mathbf{r}_{l+1}=\mathbf{r}_{l+1} - \alpha_l \mathbf{v}_l$
              \STATE $\mathbf{\hat{z}}_{l+1}=\mathbf{A}^{\adj} \mathbf{W} \mathbf{r}_{l+1}$
              \ENDFOR
  \end{algorithmic}
  Output: $\mathbf{\hat{f}}_{l}$
   \end{algorithm}

\begin{algorithm}[ht!]
    \caption{Conjugate gradients for the normal equations, Residual
      minimisation (CGNR)}
    \label{algo:wCGNR(E)}
    Input: $\mathbf{y} \in \mathbb{C}^M,\, {\mathbf{\hat{f}}}_0 \in \mathbb{C}^{|I_{\mathbf{N}}|}$
    \begin{algorithmic}[1]
      \topsep=1.5ex
      \itemsep=1.5ex
      \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.5ex
      \itemsep=1.5ex
      \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}
    Output: $\mathbf{\hat{f}}_{l}$
\end{algorithm}

%------------------------------------------------------------------------------
\subsubsection*{Interpolation}
%------------------------------------------------------------------------------

For a comparative high polynomial degree $|I_{\mathbf{N}}|>M$ one expects to
interpolate the given data $y_j\in \C,\;j=0,\hdots,M-1,$ exactly.
The (consistent) linear system \eqref{eq:infft} is under-determined.
One considers the {\em damped minimisation problem}
\begin{equation*}
  \label{eq:ip}
  \|\mathbf{\hat{f}}\|_{\mathbf{\hat{W}}}^{-1}
  =\left(\sum_{\mathbf{k}\in I_{\mathbf{N}}} \frac{|\hat f_{\mathbf{k}}|^2}{\hat w_{\mathbf{k}}}
  \right)^{1/2} \stackrel{\mathbf{\hat{f}}}{\rightarrow} \min \quad \text{subject
  to} \quad \mathbf{A} \mathbf{\hat{f}} = \mathbf{y},
\end{equation*}
which may incorporate 'damping factors' $\hat w_{\mathbf{k}}> 0$, $\mathbf{\hat{W}}:=
\diag(\hat w_{\mathbf{k}})_{\mathbf{k}\in I_{\mathbf{N}}}$.
A smooth solution is favoured, i.e., a decay of the Fourier coefficients $\hat
f_{\mathbf{k}},\;\mathbf{k}\in I_{\mathbf{N}}$, for decaying damping factors $\hat w_{
\mathbf{k}}$.
This interpolation problem is equivalent to the
{\em damped normal equation of second kind}
\begin{equation} \label{eq:no2}
  \mathbf{A} \mathbf{\hat{W}} \mathbf{A}^{\adj} \mathbf{\tilde{f}}
  = \mathbf{y}, \qquad \mathbf{\hat{f}} =
  \mathbf{\hat{W}} \mathbf{A}^{\adj} \mathbf{\tilde{f}}.
\end{equation}

Applying the conjugate gradient scheme to \eqref{eq:no2} yields the following
Algorithm \ref{algo:wCGNE(R)}.
\begin{algorithm}[ht!]
  \caption{Conjugate gradients for the normal equations, Error minimisation
    (CGNE)}
  \label{algo:wCGNE(R)}
  Input: $\mathbf{y} \in \mathbb{C}^M,\, {\mathbf{\hat{f}}}_0 \in \mathbb{C}^{|I_{\mathbf{N}}|}$
  \begin{algorithmic}[1]
    \topsep=1.5ex \itemsep=1.5ex
    \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.5ex \itemsep=1.5ex
    \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}
  Output: $\mathbf{\hat{f}}_{l}$
\end{algorithm}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\clearpage
\section{Library}\label{sect:lib}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 The library is completely written in C and uses the FFTW library \cite{fftw}
 which has to be installed on your system.
 The library has several options (determined at compile time) and parameters
 (determined at run time).

%------------------------------------------------------------------------------
\subsection{Installation}
%------------------------------------------------------------------------------
Download from {\tt www.tu-chemnitz.de/$\sim$potts/nfft} the most recent
version {\tt nfft3.x.x.tgz}. In the following, we assume that you use a \verb+bash+
compatible shell.
Uncompress the archive.
\begin{itemize}
  \item \verb"tar xfvz nfft3.x.x.tar.gz"
\end{itemize}
Change to the newly created directory.
\begin{itemize}
  \item \verb+cd nfft3.x.x+
\end{itemize}
Optional: Set environment variables (see also below).
Leave this step out, if you are unsure what to do.
Example:
\begin{itemize}
  \item \verb+export CPPFLAGS="-I/path/to/include/files"+\\
  \item \verb+export LDFLAGS="-L/path/to/libraries"+
\end{itemize}
Run the configure script.
\begin{itemize}
  \item \verb+./configure+
\end{itemize}
Compile the sources.
\begin{itemize}
  \item \verb+make+
\end{itemize}
This builds the NFFT library and all examples and applications. You can run
\verb"make install" to install the library 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 on fine grained
control over the installation directories, run the \verb"configure" script
with the option \verb+--help+, i.e. \verb+./configure --help+.
For example, to see what gets installed, run
\verb+./configure --prefix=some/path/to/a/temporary/directory+ in the line of
commands above. After \verb+make+ and \verb+make install+, you will find all
installed files in the directory structure created under the temporary directory.

The following options can also be used in conjunction with the configure script:
By default, NFFT 3 routines use the Kaiser-Bessel window function
(see Appendix \ref{sec.nfft}). You can
change the window function with the option \verb"--with-window=ARG",
where \verb"ARG" is replaced by \verb"kaiserbessel" (Kaiser-Bessel),
\verb"gaussian" (Gaussian), \verb"bspline" (B-spline),
\verb"sinc" (sinc power) or \verb"delta" (Dirac delta) to yield the
window function in parenthesis.

The NFFT transforms can be configured to measure elapsed
time for each step of Algorithm~\ref{algo:nfft} and \ref{algo:nfftH}.
This should help customise the library to one's needs.
You can enable this behaviour with
the options \verb"--enable-measure-time" and/or
\verb"--enable-measure-time-fftw".

%------------------------------------------------------------------------------
\subsection{Procedure for computing an NFFT}\label{sect:general_proc}
%------------------------------------------------------------------------------
One has to follow certain steps to write a simple program using the NFFT
library.
A complete example is shown in Section \ref{sect:first_trafo}.
The first argument of each function is a pointer to a application-owned
variable of type {\tt nfft\_plan}.
The aim of this structure is to keep interfaces small, it contains all
parameters and data.

\subsubsection*{Initialisation}
Initialisation of a plan is done by one of the {\tt nfft\_init}-functions.
The simplest version for the univariate case $d=1$ just specifies the number
of Fourier coefficients $N_0$ and the number of nonequispaced nodes $M$.
For an application-owned variable {\tt nfft\_plan my\_plan} the function call is
\begin{center}
  {\tt nfft\_init\_1d(\&my\_plan,N0,M);}
\end{center}
The first argument should be uninitialised. Memory allocation is completely
done by the init routine.

\subsubsection*{Setting nodes}
One has to define the nodes $\mathbf{x}_j \in \T^d$ for the
transformation in the member variable {\tt my\_plan.x}.
The $t$-th coordinate of the $j$-th node $\mathbf{x}_j$ is assigned by
\begin{center}
  {\tt my\_plan.x[d*j+t]= /* your choice in [-0.5,0.5] */;}
\end{center}

\subsubsection*{Precompute $\tilde\psi$}
The precomputation of the values $\tilde\psi\left({\mathbf{x}}_j - \mathbf{n}^{-1}\odot
\mathbf{l}\right)$, i.e. the entries of the matrix $\mathbf{B}$, depends on the choice
for {\tt my\_plan.x}.
If one of the proposed precomputation strategies ({\tt PRE\_LIN\_PSI}, {\tt
  PRE\_FG\_PSI}, {\tt PRE\_PSI}, or {\tt PRE\_FULL\_PSI}) is used, one has to call
the appropriate precomputation routine {\bf AFTER} setting the nodes.
For simpler usage this has been summarised by
\begin{center}
  {\tt if(my\_plan.nfft\_flags \& PRE\_ONE\_PSI)
    nfft\_precompute\_one\_psi(\&my\_plan);}
\end{center}
Note furthermore, that
\begin{enumerate}
\item {\tt FG\_PSI} and {\tt PRE\_FG\_PSI} rely on the Gaussian window
  function,
\item {\tt PRE\_LIN\_PSI} is actual node independent and in this case the
  precomputation can be done before setting nodes,
\item {\tt PRE\_FULL\_PSI} asks for quite a bit of memory.
\end{enumerate}

\subsubsection*{Doing the transform}
Algorithm \ref{algo:nfft} is implemented as
\begin{center}
  {\tt nfft\_trafo(\&my\_plan);}
\end{center}
takes {\tt my\_plan.f\_hat} as its input and overwrites {\tt my\_plan.f}.
One only needs one plan for several transforms of the same kind,
i.e. transforms with equal initialisation parameters.
For comparison, the direct calculation of \eqref{ndft} is done by
{\tt ndft\_trafo}.
The adjoint transforms are given by {\tt nfft\_adjoint} and {\tt
  ndft\_adjoint}, respectively.

All data with multi-indices is stored plain, in particular the Fourier
coefficient $\hat f_{\mathbf{k}}$ is stored in {\tt my\_plan.f\_hat[k]} with
${\text{\tt k}}:=\sum_{t=0}^{d-1} (k_t+\frac{N_t}{2}) \prod_{t'=t+1}^{d-1}
N_{t'}$.
Data vectors might be exchanged by
\begin{center}
  {\tt NFFT\_SWAP\_complex(my\_plan.f\_hat,new\_f\_hat);}.
\end{center}

\subsubsection*{Finalisation}
All memory allocated by the init routine is deallocated by
\begin{center}
  {\tt nfft\_finalize(\&my\_plan);}
\end{center}
Note, that almost all (de)allocation operations of the library are done by
{\tt fftw\_malloc} and {\tt fftw\_free}.
Additional data, declared and allocated by the application, have to be
deallocated by the user's program as well.

\subsubsection*{Data structure and functions}
 The library defines the structure {\tt nfft\_plan}, the most important
 members are listed in Table \ref{table:nfft_members}.
 Moreover, the user functions for the NFFT are collected in Table
 \ref{table:nfft_user_fun}.
 They all have return type {\tt void} and their first argument is of type {\tt
 nfft\_plan*}.


 \begin{table}[ht!]
   \begin{center}
     \begin{tabular}{|r l|l|l|}
       \hline Type                 & Name           & Size
        & Description                                              \\
       \hline
       \hline {\tt int}            & {\tt d}        & 1
        & Spatial dimension $d$                                    \\
       \hline {\tt int*}           & {\tt N}        & d
        & Multibandwidth $\mathbf{N}$                                   \\
       \hline {\tt int}            & {\tt N\_total} & 1
        & Number of coefficients $|I_{\mathbf{N}}|$                     \\
       \hline {\tt int}            & {\tt M\_total} & 1
        & Number of nodes $M$                                      \\
       \hline {\tt double complex*}& {\tt f\_hat}   &
        $|I_{\mathbf{N}}|$ & Fourier coefficients $\mathbf{\hat{f}}$ or     \\
                                   &                &
        & adjoint coefficients $\mathbf{\hat{h}}$                      \\
       \hline {\tt double complex*}& {\tt f}        & $M$
        & Samples $\mathbf{f}$                                          \\
       \hline {\tt double*}        & {\tt x}        & $dM$
        & Sampling set ${\cal X}$                                  \\
       \hline
     \end{tabular}
   \end{center}
   \caption{Interesting members of {\tt nfft\_plan}.}
   \label{table:nfft_members}
 \end{table}

 \begin{table}[ht!]
   \begin{center}
     \begin{tabular}{|l l|}
       \hline Name & Additional arguments \\
       \hline
       \hline {\tt ndft\_trafo} & \\
       \hline {\tt ndft\_adjoint} &\\
       \hline {\tt nfft\_trafo} &\\
       \hline {\tt nfft\_adjoint} &\\
       \hline {\tt nfft\_init\_1d} & {\tt int N0, int M}\\
       \hline {\tt nfft\_init\_2d} & {\tt int N0, int N1, int M}\\
       \hline {\tt nfft\_init\_3d} & {\tt int N0, int N1, int N2, int M}\\
       \hline {\tt nfft\_init} & {\tt int d, int* N, int M}\\
       \hline {\tt nfft\_init\_guru} & {\tt int d, int* N, int M, int* n, int m,} \\
       & {\tt unsigned nfft\_flags, unsigned fftw\_flags}\\
       \hline {\tt nfft\_precompute\_one\_psi} &\\
       \hline {\tt nfft\_check} &\\
       \hline {\tt nfft\_finalize} &\\
       \hline
    \end{tabular}
   \end{center}
   \caption{User functions of the NFFT.}
   \label{table:nfft_user_fun}
  \end{table}


%------------------------------------------------------------------------------
 \subsection{Generalisations and nomenclature}
%------------------------------------------------------------------------------
The generalisations to the nonequispaced FFT in time and frequency domain, the
nonequispaced fast (co)sine transform, and the nonequispaced sparse FFT are
based on the same procedure as the NFFT.
In the following we give some more details for the nonequispaced fast
spherical Fourier transform.

\subsubsection*{NFSFT}
Compared to the NFFT routines, the NFSFT module shares a similar basic
procedure for computing the transformations. It differs, however, slightly
in data layout and function interfaces.

The NFSFT routines depend on initially precomputed global data that limits the
maximum degree \verb"N" that can be used for computing transformations. This
precomputation has to be performed only once at program runtime regardless of
how many individual transform plans are used throughout the rest. This is done
by
\begin{center}
  \verb"nfsft_precompute(N,1000,0U,0U);"
\end{center}
Here, \verb"N" is the maximum degree that can be used in all subsequent transformations,
\verb"1000" is the threshold for the FPTs used internally, and
\verb"0U" and \verb"0U" are optional flags for the NFSFT and FPT. If you
are unsure which values to use, leave the threshold at \verb"1000" and
the flags at \verb"0U".

Initialisation of a plan is done by calling one of the
\verb"nfsft_init"-functions. We need no dimension parameter here, as in
the NFFT case. The simplest version is a call to \verb"nfsft_init"
specifying only the bandwidth $N$ and the number of nonequispaced nodes $M$.
For an application-owned variable  \verb"nfsft_plan my_plan" the function call
is
\begin{center}
  \verb"nfsft_init(&my_plan,N,M);"
\end{center}
The first argument should be uninitialised. Memory allocation is completely
done by the init routine.

After initialising a plan, one defines the nodes $(\vartheta_{j},\varphi_{j}) \in
[0,2\pi) \times [0,\pi]$ in spherical coordinates in the member variable
\verb"my_plan.x".
For consistency with the other modules and the conventions used there, we
use swapped and scaled spherical coordinates
\begin{align*}
  \tilde{\varphi} & :=
    \left\{
      \begin{array}{ll}
        \frac{\varphi}{2\pi}, & \text{if $0 \le \varphi < \pi$}, \\[1ex]
        \frac{\varphi}{2\pi}-1, & \text{if $\pi \le \varphi < 2\pi$},
      \end{array}
    \right. & \tilde{\vartheta} & := \frac{\vartheta}{2\pi}
\end{align*}
and identify a point from $\mathbb{S}^2$ with the vector
$(\tilde{\varphi}, \tilde{\vartheta}) \in
[-\frac{1}{2}, \frac{1}{2}) \times [0,\frac{1}{2}]$.
The angles $\tilde{\varphi}_{j}$ and $\tilde{\vartheta}_{j}$ for
the $j$-th node are assigned by
\begin{center}
  \begin{minipage}{10cm}
    \verb"my_plan.x[2*j]   = /* your choice in [-0.5,0.5] */;"
    \verb"my_plan.x[2*j+1] = /* your choice in [0,0.5] */;"
  \end{minipage}
\end{center}
After setting the nodes, the same node-dependent (see the exception above)
precomputation as for the NFFT has to be performed.
This is done by
\begin{center}
  \verb"nfsft_precompute_x(&my_plan);"
\end{center}
which itself calls \verb"nfft_precompute_one_psi(&my_plan)" as explained above. You might modify
the precomputation strategy by passing the appropriate NFFT flags to one of the
more advanced \verb"nfsft_init" routines.
Setting the spherical Fourier coefficients $\hat{f}_{k}^n$ in \verb+my_plan.f_hat+
should be done using the helper macro \verb+NFSFT_INDEX+. This reads
\begin{center}
  \begin{minipage}{11cm}
  \verb"for (k = 0; k <= N; k++)" \\
%  \verb"{" \\
  \verb"  for (n = -k; n <= k; n++)" \\
%  \verb"  {" \\
  \verb"    my_plan.f_hat[NFSFT_INDEX(k,n,&plan)] = /* your choice */;" \\
%  \verb"  }"\\
%  \verb"}"
  \end{minipage}
\end{center}
For the NFSFT, spherical Fourier coefficients are stored in a fashion different from
the NFFT. You should always use the helper macro \verb+NFSFT_INDEX+ for setting
the spherical Fourier coefficients. If, in place of a NFSFT transformation, you
would like to perform an adjoint NFSFT transformation, you set
the values $f\left(\vartheta_{j}, \varphi_{j}\right)$ in \verb+my_plan.f+ as follows:
\begin{center}
  \begin{minipage}{11cm}
  \verb"for (j = 0; j < M; j++)" \\
%  \verb"{" \\
  \verb"  my_plan.f[j] = /* your choice */;" \\
%  \verb"}"
  \end{minipage}
\end{center}

The actual transforms are computed by calling either \verb+nfsft_trafo(&my_plan)+
for the NFSFT (see Algorithm \ref{NFSFT:Algorithm:NFSFT}) or
\verb+nfsft_adjoint(&my_plan)+ for the adjoint NFSFT (see Algorithm
\ref{NFSFT:Algorithm:adjointNFSFT}).
Remember, that you must make sure that \verb+my_plan.x+ and \verb+my_plan.f_hat+
or \verb+my_plan.f+ have been initialised appropriately prior to calling the
corresponding transformation routine.
On execution, \verb+nfsft_trafo+ overwrites \verb+my_plan.f+ while
\verb+nfsft_adjoint+ overwrites \verb+my_plan.f_hat+.
One only needs one plan for several transforms of the same kind,
i.e. transforms with equal initialisation parameters.
For comparison, the direct calculation of \eqref{eq.ndsft} and \eqref{eq.adjointndsft}
are done by \verb+ndsft_trafo+ and \verb+ndsft_adjoint+, respectively.

All memory allocated by the init routine is deallocated by
\begin{center}
  \verb+nfsft_finalize(&my_plan);+
\end{center}
while the memory allocated by \verb+nfsft_precompute+ gets
freed by calling
\begin{center}
  \verb+nfsft_forget();+
\end{center}
Additional data, declared and allocated by the application, has to be
deallocated by the user's program as well.

The library defines the structure \verb+nfsft_plan+, the most important
members are listed in Table \ref{table:nfsft_members}.
The structure contains, public read-only (r) and public read-write (w)
members.
The user functions for the NFSFT are collected in Table
\ref{table:nfsft_user_fun}.

\begin{table}[ht!]
  \begin{center}
    \begin{tabular}{|r l|l|c|l|}
      \hline Type                 & Name           & Size           & Type & Description \\
      \hline
        \verb+int+                & \verb+N+       & 1              & r    & Bandwidth $N$ \\
        \verb+int+                & \verb+N_total+ & 1              & r    & Set to $4|\mathcal{I}_N|=4(N+1)^2$,\\
                            &                &                &      & quadruple of the number of Fourier coefficients\\
        \verb+int+                & \verb+M_total+ & 1              & r    & Total number of nodes $M$ \\
        \verb+double complex*+    & \verb+f_hat+   & $4|\mathcal{I}_N|$ & w    & Fourier coefficients $\hat{\mathbf{f}}$ or\\
                                  &                &                &      & adjoint coefficients $\hat{\mathbf{h}}$ \\
        \verb+double complex*+    & \verb+f+       & $M$            & w    & Samples $\mathbf{f}$ \\
        \verb+double*+            & \verb+x+       & $M$            & w    & Sampling set $\mathcal{X}$ \\
       \hline
      \end{tabular}
   \end{center}
   \caption{Interesting members of {\tt nfsft\_plan}.\label{table:nfsft_members}}
 \end{table}

 \begin{table}[ht!]
   \begin{center}
     \begin{tabular}{|l l|}
       \hline Name & Additional arguments \\
       \hline
       \hline \verb+nfsft_precompute+    & \verb+int N, double threshold,+\\
                                         & \verb+unsigned int nfsft_flags, unsigned int fpt_flags+\\
       \hline \verb+nfsft_init+          & \verb+int N, int M+ \\
       \hline \verb+nfsft_init_advanced+ & \verb+int N, int M, unsigned int nfsft_flags+ \\
       \hline \verb+nfsft_init_guru+     & \verb+int N, int M, unsigned int nfsft_flags,+ \\
                                         & \verb+int nfft_flags, int nfft_cutoff+ \\
       \hline \verb+nfsft_precompute_x+  & \\
       \hline \verb+ndsft_trafo+         & \\
       \hline \verb+ndst_adjoint+        & \\
       \hline \verb+nfsft_trafo+         & \\
       \hline \verb+nfsft_adjoint+       & \\
       \hline \verb+nfsft_finalize+      & \\
       \hline \verb+nfsft_forget+        & \\
       \hline
    \end{tabular}
   \end{center}
   \caption{User functions of the NFSFT.\label{table:nfsft_user_fun}}
  \end{table}

Some more things should be kept in mind when using the NFSFT module:
\begin{itemize}
\item The bandwidth $N_{\text{max}}$ up to which global precomputation is
  performed when calling \verb+nfsft_precompute+ is always chosen as the
  next power of two with respect to the specified maximum bandwidth $N$.
\item By default, the NDSFT and NFSFT transformations,
  \verb+ndsft_trafo+ and \verb+nfsft_trafo+, are allowed to destroy the
  input \verb+f_hat+ while the input \verb+x+ is preserved. On the contrary,
  the adjoint NDSFT and NFSFT transformations, \verb+ndsft_adjoint+ and
  \verb+nfsft_adjoint+, do not destroy the input \verb+f+ and \verb+x+ by
  default. The desired behaviour can be assured by using the
  \verb+NFSFT_PRESERVE_F_HAT+, \verb+NFSFT_PRESERVE_X+, \verb+NFSFT_PRESERVE_F+ and
  \verb+NFSFT_DESTROY_F_HAT+, \verb+NFSFT_DESTROY_X+, \verb+NFSFT_DESTROY_F+
  flags.
\end{itemize}

%------------------------------------------------------------------------------
 \subsection{Inversion and solver module} \label{sect:general_proc_infft}
%------------------------------------------------------------------------------
 The flow chart in Figure \ref{fig:struct} shows how to use the inverse
 transforms as implemented in the solver module of the library.
 There is no general stopping rule implemented, since this task is highly
 dependent on the particular application.
  \begin{figure}[ht!]
   \begin{center}
     \input{struct}
   \end{center}
   \caption{Flow chart of the inverse transforms.\label{fig:struct}}
 \end{figure}

 Each inverse transform basically wraps an already initialised Fourier
 transform.
 The user specifies one of the Algorithms
 \ref{algo:landweber}-\ref{algo:wCGNE(R)} by setting one of the flags
 ({\tt LANDWEBER}, {\tt STEEPEST\_DESCENT}, {\tt CGNR}, {\tt CGNE}).
 Weights and/or damping factors are used if the flags {\tt
 PRECOMPUTE\_WEIGHT}, {\tt PRECOMPUTE\_DAMP} are set and one has to initialise
 the members {\tt my\_iplan.w}, {\tt my\_iplan.w\_hat} in these cases.
 Default flags are {\tt CGNR}.

 In case of the NFFT, the library defines the structure {\tt infft\_plan}.
 The members of this plan are given by Table \ref{table:infft_members}.
 The user functions for the inverse NFFT are collected in Table
 \ref{table:infft_user_fun}.
 They all have return type {\tt void} and their first argument is of type {\tt
 infft\_plan*}.
 Replacing {\tt nfft} by any other Fourier transform gives the appropriate
 inverse to this transform.

 \begin{table}[ht!]
   \begin{center}
     \begin{tabular}{|r l|l|l|}
       \hline Type               & Name                & Size           & Description                       \\
       \hline
       \hline {\tt double*}      & {\tt w}             & {\tt M\_total} & weights $\mathbf{w}$ \\
       \hline {\tt double*}      & {\tt w\_hat}        & {\tt N\_total} & damping factors $\mathbf{\hat{w}}$ \\
       \hline {\tt (FLT\_TYPE)*} & {\tt y}             & {\tt M\_total} & right hand side $\mathbf{y}$ \\
       \hline {\tt (FLT\_TYPE)*} & {\tt f\_hat\_iter}  & {\tt N\_total} & actual solution \\
       \hline {\tt (FLT\_TYPE)*} & {\tt r\_iter}       & {\tt M\_total} & residual vector $\mathbf{r}_{l+1}$\\
       \hline {\tt double}       & {\tt dot\_r\_iter}  & 1 & $\|\mathbf{r}_{l+1}\|_{\mathbf{W}}^2$ \\
       \hline
     \end{tabular}
   \end{center}
   \caption{More important members of each inverse plan, where {\tt
   (FLT\_TYPE)} is {\tt double} for the NFCT and NFST and {\tt double complex}
   in all other cases.}
   \label{table:infft_members}
  \end{table}

  \begin{table}[ht!]
   \begin{center}
    \begin{tabular}{|l l|}
     \hline Name & Additional arguments \\
     \hline
     \hline {\tt infft\_init} & {\tt nfft\_plan *mv}\\
     \hline {\tt infft\_init\_advanced} & {\tt nfft\_plan *mv, unsigned infft\_flags}\\
     \hline {\tt infft\_before\_loop} & \\
     \hline {\tt infft\_loop\_one\_step} &\\
     \hline {\tt infft\_finalize} &\\
     \hline
    \end{tabular}
   \end{center}
   \caption{User functions of the inverse NFFT.}
   \label{table:infft_user_fun}
  \end{table}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\clearpage
\section{Examples}\label{sect:ex}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

The library was tested on a AMD Athlon(tm) XP 2700+, 1GB memory, SuSe-Linux,
kernel 2.4.20-4GB-athlon, gcc version 3.3.
In all tests with random input the nodes $\mathbf{x}_j$ and the Fourier
coefficients $\hat f_{\mathbf{k}}$ are chosen pseudo randomly  with  $\mathbf{x}_j \in
[-0.5,0.5]^d$ and $\hat f_{\mathbf{k}} \in [0,1] \times [0,1]{\rm i}$.

%------------------------------------------------------------------------------
\subsection{Computing your first transform}\label{sect:first_trafo}
%------------------------------------------------------------------------------

The following code summarises the steps of Section \ref{sect:general_proc} and
computes a univariate NFFT from 14 Fourier coefficients and 19 nodes.
Note that this routine is part of {\tt simple\_test.c} in {\tt examples/nfft/}
and uses additional routines as defined {\tt include/util.h} to set up and
show vectors.
\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{Computation time vs. problem size}
%------------------------------------------------------------------------------
The program {\tt nfft\_times} in the same directory compares the computation
time of the FFT (\cite{fftw}, {\tt FFTW\_MEASURE}), the straightforward
evaluation of \eqref{eq:nfft}, denoted by NDFT, and the NFFT for increasing
total problem sizes $|I_{\mathbf{N}}|$ and space dimensions $d=1,2,3$, where $
\mathbf{N}=(N,\hdots,N)^{\top},\;N\in\N$.
While the nodes for the FFT are restricted to the lattice $\mathbf{N}^{-1}\odot
I_{\mathbf{N}}$, we choose $M=N^d$ random nodes for the NDFT and the NFFT.
Within the latter, we use the oversampling factor $\sigma=2$, the cut-off
$m=4$, and the Kaiser-Bessel window function ({\tt PRE\_PSI}, {\tt
  PRE\_PHI\_HUT}).
This results in a fixed accuracy of $E_{\infty}:=\|\mathbf{f}-
\mathbf{s}\|_{\infty}/\|\mathbf{\hat{f}}\|_1\approx 10^{-8}$ for $d=1,2,3$.

\begin{table}[ht!]
  \centering
  \begin{tabular}{|r|r|r|r||r|r|r|r|}
    \hline
    $l_N$ & FFT & NDFT & NFFT &
    $l_N$ & FFT & NDFT & NFFT\\
    \hline
\multicolumn{4}{|c||}{$d=1$} & \multicolumn{4}{c|}{$d=2$}\\ \hline
$3$&    $1.3e-07$&      $8.7e-06$&      $4.6e-06$&      $6$&    $9.9e-07$&      $5.7e-04$&      $3.2e-04$ \\
$4$&    $2.0e-07$&      $3.5e-05$&      $8.7e-06$&      $8$&    $4.4e-06$&      $9.2e-03$&      $1.3e-03$ \\
$5$&    $4.0e-07$&      $1.4e-04$&      $1.7e-05$&      $10$&   $2.1e-05$&      $1.5e-01$&      $5.2e-03$ \\
$6$&    $8.9e-07$&      $5.6e-04$&      $3.6e-05$&      $12$&   $1.2e-04$&      $2.4e+00$&      $2.3e-02$ \\
$7$&    $2.2e-06$&      $2.2e-03$&      $7.2e-05$&      $14$&   $1.7e-03$&      $3.8e+01$&      $1.5e-01$ \\
$8$&    $4.8e-06$&      $9.0e-03$&      $1.4e-04$&      $16$&   $2.1e-02$&      *&              $6.8e-01$ \\
$9$&    $1.1e-05$&      $3.6e-02$&      $2.9e-04$&      $18$&   $8.4e-02$&      *&              $2.8e+00$ \\
$10$&   $2.4e-05$&      $1.4e-01$&      $6.0e-04$&      $20$&   $3.2e-01$&      *&              $1.2e+01$ \\
$11$&   $5.7e-05$&      $5.8e-01$&      $1.4e-03$&      $22$&   $1.4e+00$&
    *&              $5.3e+01$ \\ \cline{5-8}
$12$&   $1.5e-04$&      $2.3e+00$&      $3.2e-03$&
    \multicolumn{4}{c|}{$d=3$}\\ \cline{5-8}
$13$&   $5.5e-04$&      $9.4e+00$&      $8.2e-03$&      $9$&    $1.0e-05$&      $3.7e-02$&      $2.5e-02$ \\
$14$&   $1.7e-03$&      $3.8e+01$&      $2.0e-02$&      $12$&   $1.1e-04$&      $2.4e+00$&      $2.5e-01$ \\
$15$&   $3.8e-03$&      $1.5e+02$&      $4.9e-02$&      $15$&   $3.4e-03$&      $1.5e+02$&      $2.4e+00$ \\
$16$&   $8.2e-03$&      *&              $1.2e-01$&      $18$&   $5.2e-02$&      *&              $2.1e+01$ \\
$17$&   $1.9e-02$&      *&              $2.4e-01$&      $21$&   $9.0e-01$&      *&              $1.8e+02$ \\
$18$&   $4.5e-02$&      *&              $3.6e-01$&       & & & \\
$19$&   $9.2e-02$&      *&              $9.8e-01$&       & & & \\
$20$&   $1.9e-01$&      *&              $2.1e+00$&       & & & \\
$21$&   $4.2e-01$&      *&              $4.2e+00$&       & & & \\
$22$&   $1.0e-00$&      *&              $9.5e+00$&       & & & \\
    %\input{images/nfft/timings_nice_dagobert}
    \hline
  \end{tabular}
  \caption{Computation time in seconds with respect to $l_N=\log_2|I_{
  \mathbf{N}}|$.
    Note that we used accumulated measurements in case of small times and the
    times (*) are not displayed due to the large response time in comparison
    to the FFT time.}
  \label{tab:TimeSpace_1}
\end{table}

We conclude the following: The FFT and the NFFT show the expected
$\cO(|I_{\mathbf{N}}|\log|I_{\mathbf{N}}|)$ time complexity, i.e., doubling the total size $|I_{\mathbf{N}}|$
results in only slightly more than twice the computation time, whereas the
NDFT behaves as $\cO(|I_{\mathbf{N}}|^2)$.
Note furthermore, that the constant in the $\cO$-notation is independent of
the space dimension $d$ for the FFT and the NDFT, whereas the computation time
of the NFFT increases considerably for larger $d$.

%------------------------------------------------------------------------------
\subsection{Accuracy vs. window function and cut-off parameter $m$}
%------------------------------------------------------------------------------
The accuracy of the Algorithm \ref{algo:nfft}, measured by
\begin{equation*}
E_2=\frac{\|\mathbf{f}- \mathbf{s}\|_2}{\|\mathbf{f}\|_2}
=\left({\sum_{j=0}^{M-1}|f_j-s\left(
\mathbf{x}_j\right)|^2}/{\sum_{j=0}^{M-1}
    |f_j|^2}\right)^{\frac{1}{2}}
\end{equation*}
and
\begin{equation*}
E_{\infty}=\frac{\|\mathbf{f}- \mathbf{s}\|_{\infty}}{\|\mathbf{\hat{f}}\|_1}
={\max_{0\le j<M}|f_j-s\left(
\mathbf{x}_j\right)|}/{\sum_{\mathbf{k} \in I_{\mathbf{N}}} |\hat
  f_{\mathbf{k}}|}
\end{equation*}
is shown in Figure \ref{fig:accuracy}.

 \begin{figure}[ht!]
   \begin{center}
    \begin{tabular}{c c c}
     \includegraphics[width=4.8cm]{images/accuracy1} &
     \includegraphics[width=4.8cm]{images/accuracy2} &
     \includegraphics[width=4.8cm]{images/accuracy3} \\
     \includegraphics[width=4.8cm]{images/accuracy4} &
     \includegraphics[width=4.8cm]{images/accuracy5} &
     \includegraphics[width=4.8cm]{images/accuracy6}
    \end{tabular}
   \end{center}
   \caption{The error $E_2$ (top) and $E_{\infty}$ (bottom) with respect to $m$,
     from left to right $d=1,2,3$ ($N=2^{12},2^6,2^4,\, \sigma=2,\,M=10000$),
     for Kaiser Bessel- (circle), Sinc- (x), B-Spline- ($+$), and Gaussian window (triangle).}
   \label{fig:accuracy}
 \end{figure}

%------------------------------------------------------------------------------
\subsection{Computing an inverse transform}\label{sect:rec_example}
%------------------------------------------------------------------------------
The usage of the inverse NFFT is demonstrated by {\tt simple\_test} in {\tt
  examples/solver}.
Executing the {\sc MatLab} script {\tt glacier.m} in the same  directory
produces the following two plots.
Note that the corresponding C-file {\tt glacier.c} is called from the {\sc
  MatLab} script.
\begin{figure}[ht!]
  \begin{center}
    \begin{tabular}{cc}
      \includegraphics[width=6cm]{images/glacier1} &
      \includegraphics[width=6cm]{images/glacier2}
    \end{tabular}
  \end{center}
  \caption{Reconstruction of the glacier from samples at $M=8345$ nodes ({\tt
      vol87.dat} from \cite{Franke.Daten}) with $N_0=N_1=256$  and $40$
      iterations. \label{Fig:glacier}}
\end{figure}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\clearpage
\section{Applications}\label{sec:appl}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

In this section we describe important applications which are based on the NFFT
kernel.
One can find these programs in the directory {\tt applications}.

%------------------------------------------------------------------------------
\subsection{Summation of smooth and singular kernels}
%------------------------------------------------------------------------------

We are interested in the fast evaluation of linear combinations of radial
functions, i.e. the computation of
\begin{equation*}
  g\left(\mathbf{y}_j\right) := \sum_{k=1}^N \alpha_k K\left(\left\|\mathbf{y}_j-
  \mathbf{x}_k\right\|_2\right)
\end{equation*}
for $j=1,\hdots,M$ and nodes $\mathbf{x}_k, \mathbf{y}_j \in \R^d$.
For smooth kernels $K$ with an additional parameter $c>0$, 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}$ our algorithm requires $\cO(N+M)$
arithmetic operations.
In the case of 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*}
an additional regularisation procedure must be incorporated and the algorithm
has the arithmetic complexity  $\cO(N \log N +M)$ or $\cO(M\log M + N)$ if
either the target nodes $\mathbf{y}_j$ or the source nodes $\mathbf{x}_k$ are
``reasonably uniformly distributed''.

Note that the proposed fast algorithm \cite{post02,postni04,FeSt03}
generalises the diagonalisation of convolution matrices by Fourier matrices to
the setting of arbitrary nodes.
In particular, this yields nearly the same arithmetic complexity as the FMM
\cite{beatgreen97} while allowing for an easy change of various kernels.
A recent application in particle simulation is given in \cite{popori05}.
The directory {\tt applications/fastsum} contains C and {\sc MatLab} programs
that show how to use the fast summation method.

%------------------------------------------------------------------------------
\subsection{Fast Gauss transform}
%------------------------------------------------------------------------------

This is a special case of the fast summation method, we compute approximations
of the following sums.
Given complex coefficients $\alpha_k \in \C$ and source nodes
$x_k\in[-\frac{1}{4},\frac{1}{4}]$, our goal consists in the fast evaluation
of the sum
\begin{equation*}
  g\left(y\right)=\sum_{k=1}^N \alpha_k {\rm e}^{-\sigma|y-x_k|^2}
\end{equation*}
at the target nodes $y_j \in [-\frac{1}{4},\frac{1}{4}]$, $j=1,\ldots,M$,
where  $\sigma =a +{\rm i} b$, $a>0$, $b\in\R$, denotes a  complex parameter.
For details see \cite{KuPoSt06} and the related paper \cite{AnBe04} for
applications.
All numerical examples of \cite{KuPoSt06} are produced by the programs in {\tt
  applications/fastgauss}.

%------------------------------------------------------------------------------
\subsection{Summation of zonal functions on the sphere}
%------------------------------------------------------------------------------

Given $M,N \in \N$, arbitrary source nodes $\mathbf{\eta}_k \in \S^2$ and real
coefficients $\alpha_k\in\R$, evaluate the sum
\begin{equation*}
  g\left(\mathbf{\xi}\right) := \sum_{k = 1}^N \alpha_k K\left(\mathbf{\eta}_k \cdot
  \mathbf{\xi}\right)
\end{equation*}
at the target nodes $\mathbf{\xi}_j\in\S^2$, $j=1,\hdots,M$.
The naive approach for evaluating this sum takes $\cO(MN)$ floating point
operations if we assume that the zonal function $K$ can be evaluated in
constant time or that all values $K(\mathbf{\eta}_k \cdot \mathbf{\xi}_j)$ can be stored
in advance.

In contrast, our scheme is based on the nonequispaced fast spherical Fourier
transform, has arithmetic complexity $\cO(M+N)$, and can be easily adapted to
such different kernels $K$ as
  \begin{enumerate}
  \item the Poisson kernel
    $Q_{h}:[-1,1] \rightarrow \R$ with $h \in (0,1)$  given by
    \[
    Q_{h}(x) := \frac{1}{4\pi} \frac{1-h^2}{\left(1-2 h x+h^2\right)^{3/2}},
    \]
  \item  the singularity kernel
    $S_{h}:[-1,1] \rightarrow \R$ with $h \in (0,1)$  given by
    \[
    S_{h}(x) := \frac{1}{2\pi} \frac{1}{\left(1-2 h x +h^2\right)^{1/2}},
    \]
\item the locally supported kernel
 $L_{h,\lambda}:[-1,1] \rightarrow \R$ with $h \in (-1,1)$ and
 $\lambda \in \mathbb{N}_0$ given by
 \[
    L_{h,\lambda}(x) :=
    \left\{\begin{array}{l@{\quad \text{if} \quad}rcl}
        0 & -1 &\le x \le& h, \\
        \frac{\lambda+1}{2\pi(1-h)^{\lambda+1}}\left(x-h\right)^{\lambda} &
        h & <  x \le& 1,
      \end{array}\right.
    \]
or
\item the spherical Gaussian kernel
  $G_{\sigma}:[-1,1] \rightarrow
  \R$ with $\sigma>0$
  \[
  G_{\sigma}(x) := \e^{2\sigma x-2\sigma}\,.
  \]
  \end{enumerate}

For details see \cite{KeKuPo06}, all corresponding numerical examples can be
found in {\tt applications\ /fastsumS2}.

%------------------------------------------------------------------------------
\subsection{Iterative reconstruction in magnetic resonance imaging}
%------------------------------------------------------------------------------

In magnetic resonance imaging (MRI) the raw data is measured in k-space, the
domain of spatial frequencies.
Methods that use a non-Cartesian sampling grid in k-space, e.g. a spiral, are
becoming increasingly important.
Reconstruction is usually performed by resampling the data onto a Cartesian
grid and the usage of the standard FFT - often denoted by gridding.
Another approach, the inverse model, is based on an implicit discretisation.
Both discretisations are solved efficiently by means of the NFFT and the
inverse NFFT, respectively.
Furthermore, a unified approach to field inhomogeneity correction has been
included, see \cite{KnKuPo,EgKnPo06} for details.

%------------------------------------------------------------------------------
\subsection{Computation of the polar FFT}
%------------------------------------------------------------------------------

The polar FFT is a special case of the NFFT, where one computes the Fourier
transform on particular grids.
Of course, the polar as well as a so-called pseudo-polar FFT can be computed
very accurately and efficiently by 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{FeKuPo06} and for further applications
\cite{AvCoDoElIs}.

\begin{figure}[ht]
  \centering
  \subfigure
  {\includegraphics[width=0.3\textwidth]{images/polar_grid}}\hfill
  \subfigure
  {\includegraphics[width=0.3\textwidth]{images/modified_polar_grid}}\hfill
  \subfigure
  {\includegraphics[width=0.3\textwidth]{images/linogram_grid}}
 \caption{Left to right: polar, modified polar, and linogram grid of size
   $R=16$, $T=32$. \label{fig:grids}}
\end{figure}

%------------------------------------------------------------------------------
\subsection{Radon transform, computer tomography, and ridgelet transform}
%------------------------------------------------------------------------------

We are interested in efficient and high quality reconstructions of digital $N
\times N$ medical images from their Radon transform.
The standard reconstruction algorithm, the filtered backprojection, ensures a
good quality of the 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, the straightforward Fourier reconstruction algorithm suffers
from unacceptable artifacts so that it is useless in practice.
A better quality of the reconstructed images can be achieved by our algorithm
based on NFFTs.
For details see \cite{PoSt01IMA, PoSt00SPIE, PoSt02JFAA} and the directory
{\tt applications/radon}.

Another application of the discrete Radon transform is the discrete Ridgelet
transform, see e.g. \cite{CaDeDoYi}.
A simple test program for denoising an image by hard thresholding the ridgelet
coefficients can be found in {\tt applications/radon}.
It uses the NFFT-based discrete Radon transform and the translation-invariant
discrete Wavelet transform based on {\sc MatLab} toolbox WaveLab850
\cite{WaveLab}.
See \cite{FeMa05} for details.

\bibliographystyle{abbrv}
%\bibliographystyle{alpha}
\bibliography{../references/references}
\end{document}
