NFSFT - NFFT on the sphere

This library of C functions computes evaluates a function $f\in {\mathrm{L}}^2({\mathbb{S}}^2)$ with finite orthogonal expansion $$f(\vartheta ,\phi )=\sum _{k=0}^M\sum _{n=-k}^k{a}_k^n{Y}_k^n(\vartheta ,\phi )$$ in terms of spherical harmonics $Y_k^n$ on a set of arbitary nodes $({\vartheta }_d,{\phi }_d)$, $d=1,\dots ,D$, $D\in \mathbb{N}$, in spherical coordinates. Furthermore, the fast evaluation of sums $${\tilde{a}}_k^n:=\sum _{d=1}^{D}f({\vartheta }_d,{\phi }_d)\overline{Y_k^n({\vartheta }_d,{\phi }_d)}$$ for given function values $f({\vartheta }_d,{\phi }_d)\in \mathbb{C}$ and all indices $k=0,\dots ,M$, $n=-k,\dots ,k$ is possible. The algorithm is also kwnon as fast spherical harmonic transform for nonequispaced data.

Detailed Description

Fourier analysis on the sphere has practical relevance in tomography, geophysics, seismology, meteorology and crystallography. In analogy to the complex exponentials ${\mathrm{e}}^{\mathrm{i}kx}$ on the torus, the spherical harmonics form the orthogonal Fourier basis with respect to the usual inner product on the sphere.

Spherical coordinates

Every point in ${\mathbb{R}}^3$ can be described in spherical coordinates by a vector $(r,\vartheta ,\phi {)}^{\mathrm{\top }}$ with the radius $r\ge 0$ and two angles $\vartheta \in [0,\pi ]$, $\phi \in [0,2\pi )$. We denote by ${\mathbb{S}}^2$ the two-dimensional unit sphere embedded into ${\mathbb{R}}^3$, i.e. $${\mathbb{S}}^2:=\{\mathbf{x}\in {\text{R}}^3:\Vert \mathbf{x}{\Vert }_2=1\}$$ and identify a point from ${\mathbb{S}}^2$ with the corresponding vector $(\vartheta ,\phi {)}^{\mathrm{\top }}$.

Legendre polynomials and associated Legendre functions

The Legendre polynomials $P_k:[-1,1]\to \text{R},k\ge 0$, as classical orthogonal polynomials are given by their corresponding Rodrigues formula $$P_k(x):=\frac{1}{2^{k}k!}\frac{{\text{d}}^k}{\text{d}{x}^k}{(x^2-1)}^k.$$ The associated Legendre functions $P_k^n:[-1,1]\to \text{R},k\ge n\ge 0$ are defined by $$P_k^n(x):={\left(\frac{(k-n)!}{(k+n)!}\right)}^{1/2}{(1-x^2)}^{n/2}\frac{{\text{d}}^n}{\text{d}{x}^n}{P}_k(x).$$ 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 $$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)$$ for $k\ge n\ge 0,P_{n-1}^n(x)=0,P_n^n(x)=\frac{\sqrt{(2n)!}}{2^{n}n!}{(1-x^2)}^{n/2}$. For fixed $n$, the set $\{P_k^n:k\ge n\}$ forms a set of orthogonal functions, i.e., $$⟨P_k^n,P_l^n⟩={\int }_{-1}^1{P}_k^n(x)P_l^n(x)\text{d}x=\frac{2}{2k+1}{\delta }_{k,l}.$$ Again, we denote by ${\overline{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.

Spherical harmonics

The spherical harmonics $Y_k^n:{\mathbb{S}}^2\to \mathbb{C},k\ge |n|,n\in \mathbb{Z}$, are given by $$Y_k^n(\vartheta ,\phi ):=P_k^n(\cos \vartheta ){\mathrm{e}}^{\mathrm{i}n\phi }.$$ They form an orthogonal basis for the space of square integrable functions on the sphere, i.e., $$⟨Y_k^n,Y_l^m⟩={\int }_0^{2\pi }{\int }_0^{\pi }{Y}_k^n(\vartheta ,\phi )\overline{Y_l^m(\vartheta ,\phi )}\sin \vartheta \mathrm{d}\vartheta \mathrm{d}\phi =\frac{4\pi }{2k+1}{\delta }_{k,l}{\delta }_{n,m}.$$ The orthonormal spherical harmonics are denoted by ${\overline{Y}}_k^n=\sqrt{\frac{2k+1}{4\pi }}{Y}_k^n$. Hence, any square integrable function $f:{\mathbb{S}}^2\to \mathbb{C}$ has the expansion $$f=\sum _{k=0}^{\mathrm{\infty }}\sum _{n=-k}^k\hat{f}(k,n){\overline{Y}}_k^n,$$ with the spherical Fourier coefficients $\hat{f}(k,n)=⟨f,{\overline{Y}}_k^n⟩$. The function $f$ is called bandlimited, if $\hat{f}(k,n)=0$ for $k>N$ and some $N\in \text{N}$. In analogy to the NFFT on the Torus $\mathbb{T}$, we define the sampling set $$\mathcal{X}:=\{({\vartheta }_j,{\phi }_j)\in {\mathbb{S}}^2:j=0,\dots ,M-1\}$$ of nodes on the sphere ${\mathbb{S}}^2$ and the set of possible ''frequencies'' $${\mathcal{I}}_N:=\{(k,n):k=0,1,\dots ,N;n=-k,\dots ,k\}.$$ The nonequispaced discrete spherical Fourier transform (NDSFT) is defined as the evaluation of the finite spherical Fourier sum $$f_j=f({\vartheta }_j,{\phi }_j)=\sum _{(k,n)\in {\mathcal{I}}_N}{\hat{f}}_k^n{Y}_k^n({\vartheta }_j,{\phi }_j)=\sum _{k=0}^N\sum _{n=-k}^k{\hat{f}}_k^n{Y}_k^n({\vartheta }_j,{\phi }_j)$$ for $j=0,\dots ,M-1$. In matrix vector notation, this reads $\mathbf{f}=\mathbf{Y}\hat{\mathbf{f}}$ with $$\mathbf{f}:={\left(f_j\right)}_{j=0}^{M-1}\in {\mathbb{C}}^M,f_j:=f({\vartheta }_j,{\phi }_j),$$ $$\mathbf{Y}:={\left(Y_k^n({\vartheta }_j,{\phi }_j)\right)}_{j=0,\dots ,M-1;(k,n)\in {\mathcal{I}}_N}\in {\mathbb{C}}^{M\times (N+1{)}^2},$$$$\hat{\mathbf{f}}:={\left({\hat{f}}_k^n\right)}_{(k,n)\in {\mathcal{I}}_N}\in {\mathbb{C}}^{(N+1{)}^2}.$$ The corresponding adjoint nonequispaced discrete fast spherical Fourier transform (adjoint NDSFT) is defined as the evaluation of $${\hat{h}}_k^n=\sum _{j=0}^{M-1}f({\vartheta }_j,{\phi }_j)\overline{Y_k^n({\vartheta }_j,{\phi }_j)}$$ for all $(k,n)\in {\mathcal{I}}_N$. Again, in matrix vector notation, this reads $\hat{\mathbf{h}}={\mathbf{Y}}^H\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,\dots ,M-1$, and sufficient degree of exactness is available, one might recover the Fourier coefficients $\hat{f}(k,n)$ by evaluating $${\hat{f}}_k^n={\int }_0^{2\pi }{\int }_0^{\pi }f(\vartheta ,\phi )\overline{Y_k^n(\vartheta ,\phi )}\sin \vartheta \mathrm{d}\vartheta \mathrm{d}\phi =\sum _{j=0}^{M-1}{w}_{j}f({\vartheta }_j,{\phi }_j)\overline{Y_k^n({\vartheta }_j,{\phi }_j)}$$ for all $(k,n)\in {\mathcal{I}}_N$. Direct algorithms for computing the NDSFT and adjoint NDSFT transformations need $\mathcal{O}(M{N}^2)$ arithmetic operations. A combination of the fast polynomial transform and the NFFT leads to approximate algorithms with $\mathcal{O}(N^2{\log }^{2}N+M)$ arithmetic operations. These are denoted NFSFT and adjoint NFSFT, respectively.
 

Algorithm

The algorithms are implemented by Jens Keiner in ./kernel/nfsft. The OpenMP parallelization was implemented by Toni Volkmer. The Matlab interface was implemented by Stefan Kunis in ./matlab/nfsft. Related paper are
 

• Gr�f, M., Kunis, S., Potts, D.
  On the computation of nonnegative quadrature weights on the sphere.
  Appl. Comput. Harm. Anal. 27(1),   (full paper ps, pdf),   2009
   
• Keiner, J., Potts, D.
  Fast evaluation of quadrature formulae on the sphere
  Math. Comput. 77, 397 - 419,   (full paper ps, pdf),   2008
   
• Kunis, S. and Potts, D.
  Fast spherical Fourier algorithms.
  J. Comput. Appl. Math. 161, 75-98. (full paper ps, pdf),   2003
   
• Potts, D., Steidl G., and Tasche M.
  Fast and stable algorithms for discrete spherical Fourier transforms.
  Linear Algebra Appl. 275, 433-450. (full paper ps.Z, pdf),   1998