Numerical Approximation of the Continuous Fourier Transform Given a function $F(k)$ in frequency space (sufficiently nice enough, eg. a Gaussian), I would like to compute its Fourier inverse 
\begin{equation}f(x) = \frac{1}{2\pi}\int_{-\infty}^{\infty}e^{ikx}F(k)dk\end{equation}
numerically (there is no explicit analytical formula), at some specified evenly distributed points $x_n$. 
Assume $F(k)$ is symmetric about $k=0$ and "essentially zero" outside the interval $-a/2\leq k \leq a/2$.
Naively, one way I would proceed is by dividing the interval $(-a/2,a/2)$ into $N$ subintervals (choose $N$ to be a power of two) and then piecewise approximate the integral: 
$$ f(x_n) \approx \frac{1}{2\pi}\frac{a}{N}\sum_{m=0}^{N-1} e^{ik_mx_n}F(k_m)$$
where, for example, we might take $k_m = (m-N/2)\frac{a}{N}$for $0\leq m\leq N-1$.
For efficiency I would like to use the Fast Fourier Transform (FFT) to approximate $f(x_n)$. So I need to put $f(x)$ into the correct format for the Discrete Fourier Transform (DFT) $$L_k = \frac{1}{N}\sum_{m=0}^{N-1}G_me^{(2\pi i) km/N},\,\,\,\,k = 0,\ldots,N-1$$
where $G_m$ is the $m$th sampling of some given function $G$.
To proceed, we divide the interval $(-a/2,a/2)$ into $N$ pieces ($N$ a power of two) and let $k_m = (m-N/2)\frac{a}{N}$for $0\leq m\leq N-1$, as before. Then
$$f(x) \approx \frac{1}{2\pi}\frac{a}{N}\sum_{m=0}^{N-1} e^{ik_mx}F(k_m)$$
However, in order for this formula to be in the correct format to use the DFT, we cannot evaluate f at any $x$: I believe they must be of the form $x_k = \frac{2\pi}{a}(k-N/2), 0\leq k \leq N-1$. 
So in order to use the FFT, the points at which I can "evaluate" the function $f$ is determined by how I sample my function $F(k)$ and the number of sample points.
But in the naive method, I can, in principle, find $f(x)$ for any given x, regardless of how many sample points I take of the function $F(k)$.
As you can tell my numerics knowledge about this is shaky at best.
My questions are:


*

*Is it even possible to use the FFT, if I need to compute f at specified points $x_n$?

*If not, are there other numerical methods that would be more suitable?
 A: $\def \o {\omega}$
(Here I use $t$ and $\o$ instead $x$ and $k$, but is equivalent)
We want to calculate numerically the inverse Fourier transform of a given 
function $F(\o)$ . By definition
$f(t)=\frac{1}{2\pi} \int_{-\infty}^{\infty} F(\o) e^{-i\o t} d\o $
We can split it into two integrals
$f(t)=\frac{1}{2\pi} \left( \int_{0}^{\infty} F(\o) e^{-i\o t} d\o + \int_{0}^{\infty} F(-\o) e^{i\o t} d\o \right) $ 
If the original function was real then $F(-\o)=F(\o)^{*}$, and therefore
$f(t)=\frac{1}{\pi}  \textbf{real} \left( \int_{0}^{\infty} F(\o) e^{-i\o t} d\o \right)$
Then we use the Riemann integral definition for $0$ to a very great frequency $W$.
$I=\int_{0}^{W} F(\o) e^{-i\o t} d\o\approx  \sum_{k=0}^{N-1} F(w_k) e^{-i \Delta \o k \Delta t n}  \Delta\o $
where $\Delta w =W/N$ so that $w_k= \Delta w \ \ k$  with  $k=0,1,2,..., N-1$ . We discretized the time by 
$t_n= \Delta t \ \ n$ with $n=0,1,2,...,N-1$. Next we make the basic assumption that
$ W \Delta t = 2 \pi $
so we end up with
$I =   \Delta\o  \sum_{k=0}^{N-1} F(w_k) e^{- 2 \pi i  \frac{n k}{N} }   $
which is the basic definition of the Fast Fourier Transform. The problem was when I used it for a square function
$F(\o)=\frac{1}{i\o}\left( e^{i\o T_0} -1 \right) $
This method gives an offset value.  Then I modified the integral as
$I=  \Delta\o \left( \sum_{k=1}^{N-1} F(w_k)  e^{- 2 \pi i  \frac{n k}{N} }  + \frac{1}{2}(F(w_0)+F(w_N)) \right) $
Which is the trapezoidal rule instead of the Riemman integral, and gave a better result.
Best regards.
E. Gutierrez-Reyes 
A: The following two articles may be of use.  They don't let you do the transform at any old x values.  But it does allow transformation from m equally spaced points in the w domain to m equally spaced points in the x domain.  The Bailey and Swarztrauber article also uses a gaussian function like yours as an example.  Sorry I can't be more specific but I'm learning all this stuff myself.
Bailey, D., and P. Swarztrauber. 1994. ‘A Fast Method for the Numerical Evaluation of Continuous Fourier and Laplace Transforms’. SIAM Journal on Scientific Computing 15 (5): 1105–10. doi:10.1137/0915067.
Inverarity, G. 2002. ‘Fast Computation of Multidimensional Fourier Integrals’. SIAM Journal on Scientific Computing 24 (2): 645–51. doi:10.1137/S106482750138647X.
A: One way that comes to mind would be to compute both approximations for $f(x_n)$ and for $f'(x_n)$ one some grid $x_n$ that allows you to use FFT. To compute the approximations for $f'(x_n)$, you can use the identity $\mathcal{F}(f')(k) = ik\mathcal{F}(f)(k)$ (i.e., you get the fourier transform of the derivative by multiplying with $ik$).
To evaluate $f(x)$ at $y$, you'd then find the closest $x_n,x_{n+1}$ with $x_n \leq y < x_{n+1}$ and use polynomial interpolation (with polynomials of degree 3) to find $f(y)$ from $f(x_n)$, $f(x_{n+1})$, $f'(x_n)$, $f'(x_{n+1})$. 
A: As long as F(k) is limited and zero outside of $-2\pi A$ to $2\pi A$ the Shannon interpolation can be used
$$f(x)=\sum_{n=-\infty}^{n=\infty} f(x_n) sinc\left(\frac{\pi}{T}(x-nT)\right) $$
where $T$ is the sampling period used to determine $x_n$ and $sinc(t)=sin(t)/t$
This formula is exact as long as function spectrum ($F(k)$) is limited.
Implementation of fast sinc interpolation exist
