Nonequispaced fast Fourier transform
Interpolation

# Interpolation

Let $${\boldsymbol{N}}=(N_1,\ldots,N_d)^{\top}\in 2\mathbb{N}^d$$. For given samples $$({\boldsymbol{x}}_j,y_j)\in \mathbb{T}^d\times\mathbb{C},\; j=0,\dots,M-1$$, the index set $$I_{{\boldsymbol{N}}}^d:= \big\{-\frac{N_1}{2},\dots,\frac{N_1}{2}-1\big\}\times \dots\times\big\{-\frac{N_d}{2},\dots,\frac{N_d}{2}-1\big\}$$, of frequencies, we construct a $$d$$-variate trigonometric polynomial $f\left({\boldsymbol{x}}\right):=\sum\limits_{{\boldsymbol{k}} \in I_{{\boldsymbol{N}}}^d} \hat f_{{\boldsymbol{k}}} {\rm e}^{ 2\pi{\mathrm i} {\boldsymbol{k}} {\boldsymbol{x}}}$ such that $$f({\boldsymbol{x}}_j) \approx y_j,\;j=0,\dots,M-1$$. Turning this into matrix vector notation, we aim to solve the system of linear equations ${\boldsymbol{A}} {\boldsymbol{\hat f}} \approx {\boldsymbol{y}}\tag{1}$ for the unknown vector of Fourier coefficients $${\boldsymbol{\hat f}}:=(\hat f_{{\boldsymbol{k}}} )_{{\boldsymbol{k}}\in I_{{\boldsymbol{N}}}^d} \in \mathbb{C}^{N^d}$$. We denote the vector of the given sample values by $${\boldsymbol{y}}:=(y_{j})_{j=0,\ldots,M-1}\in \mathbb{C}^{M}$$ and the nonequispaced Fourier matrix by ${\boldsymbol{A}}:=\left({\rm e}^{ 2\pi{\mathrm i} {\boldsymbol{k}} {\boldsymbol{x}}_j}\right)_{j=0,\dots,M-1;{\boldsymbol{k}} \in I_{{\boldsymbol{N}}}^d} \in \mathbb{C}^{M\times N^d}.$

We focus on the under-determined and consistent linear system $${\boldsymbol{A}} {\boldsymbol{\hat f}} = {\boldsymbol{y}}$$, i.e., we expect to interpolate the given data $$y_j\in\mathbb{C}$$, $$j=0, \dots,M-1$$, exactly.

In particular, we incorporate damping factors $$\hat w_{{\boldsymbol{k}}}> 0$$, $${\boldsymbol{k}} \in I_{{\boldsymbol{N}}}^d$$, and consider the optimal interpolation problem $\|{\boldsymbol{\hat f}}\|_{{\boldsymbol{\hat W^{-1}}}}^2 =\sum_{{\boldsymbol{k}}\in I_{{\boldsymbol{N}}}^d} \frac{|\hat f_{{\boldsymbol{k}}}|^2}{\hat w_{{\boldsymbol{k}}}} \stackrel{{\boldsymbol{\hat f}}}{\rightarrow} \min \quad \text{subject to} \quad {\boldsymbol{A}} {\boldsymbol{\hat f}} = {\boldsymbol{y}}, \tag{2}$ where $${\boldsymbol{\hat W}}:= {\rm diag}(\hat w_{{\boldsymbol{k}}})_{{\boldsymbol{k}}\in I_{{\boldsymbol{N}}}^d}$$.

This library of C functions computes approximations of (2) with the CGNE method.

The algorithms are implemented by Stefan Kunis in ./solver. Related paper are
• Kunis, S. and Potts, D.
Stability Results for Scattered Data Interpolation by Trigonometric Polynomials.
SIAM J. Sci. Comput. 29, 1403 - 1419,   (full paper ps, pdf)   2007