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.
You can download it here NFFT Library (Nonequispaced Fast Fourier transform):
http://www-user.tu-chemnitz.de/~potts/nfft/
Enjoy.
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)) ?
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.
