Non Uniform DFT Algorithm by Dutt & Rokhlin

Question

I have implemented Dutt & Rokhlin's FFT algorithm [1] for non-equispaced data but for some reason I am getting very large errors ($E_2=0.04$ and $E_\infty=0.1$). I was wondering if there are any freely available implementation of this method so I can compare my implementation.

[1] Dutt, Alok, and Vladimir Rokhlin. "Fast Fourier transforms for nonequispaced data." SIAM Journal on Scientific computing 14.6 (1993): 1368-1393.

Answer 1

You may download the NFFT Library - Non Equispaced Fast Fourier Transform.
It has interfaces for MATLAB, Octave and Julia.

Answer 2

After more testing I think I found a problem with the way the algorithm is formulated. According to Eq. (60)

$$ \sum_{k=0}^N\alpha_k \sum_{j=-q/2}^{q/2}P_{jk}\, e^{i(\mu_k+j)x/m}=\sum_{j=-mN/2}^{mN/2-1}\tau_j\, e^{ijx/m} $$.

The idea is to compute the r.h.s. using an inverse FFT. The problem with this is is that some terms of the sum on the l.h.s. are left out: I have verified that when these terms are included by explicitly computing the sum (no FFT involved), the relative approximation error goes down to 1e-3. So the question is: how does one compute the l.h.s. sum using an $m N$ points inverse DFT with kernel $e^{2\pi i k j/(m N)}$ (Eq. (62)) ?

Answer 3

The correct way to implement the sum in the previous answer is

for $l=-m N/2$ to $m N/2-1$

$\quad$ for $k=0$ to $N$

$\quad\quad$ for $j=-q/2$ to $q/2$

$\quad\quad\quad$ if $l==\mod(\mu_k+j+m N/2,m N)-m N/2$,

$\quad\quad\quad\quad$ $\tau_l=\tau_l+\alpha_k P_{jk}$

$\quad\quad\quad$ end

$\quad\quad$ end

$\quad$ end

In this way all the terms $\alpha_k P_{jk}$ are accounted for in the sum.