# 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 $$f(x) = \frac{1}{2\pi}\int_{-\infty}^{\infty}e^{ikx}F(k)dk$$ 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:

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

2. If not, are there other numerical methods that would be more suitable?

• I came across your very same problem a week ago and I was quite astonished in realizing that there are no routines to do the job in known (at least those that I know) computer algebra programs (e.g. Mathematica/Matlab). The best thing I found is numerical recipes. There is a whole section on Fourier integration. Basically I just copied the routines suggested there and now I'm an happy person. I think that's the best what the market offers at least according to my experience. – lcv May 11 '13 at 0:48
• Bythe way, the routine in numerical recipes implements fgp's trick and more. – lcv May 11 '13 at 22:27
• @lcv: can you share the link to the routine you found? I also came across this problem and I'm still shocked that this isn't some standard Matlab routine. It took some time to realize that the maximum range of the x_k numbers is increasing with larger N. Increasing N will therefore make sure that the x_k's cover the entire x-region of interest but will not increase the density inside the x-region. – Frank Meulenaar Dec 5 '13 at 8:02
• Of course, one can increase the densitiy of the x_k's is by increasing a. – Frank Meulenaar Dec 5 '13 at 8:29
• @FrankMeulenaar It's numerical recipes chapter 13.9 "Computing Fourier integrals using FFT" . In the 'C' edition is on page 584. You can find online and downloadable versions of NR. The routine you can pretty much copy-paste it from the pdf. You'll have to work a bit to remove the comments. Otherwise the routines are available online for purchase on the NR site. Good luck! – lcv Dec 6 '13 at 20:30

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.

$$\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

• It should say -infinity to infinity, it was a mistake. – Edahi Aug 26 '19 at 19:24

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})$.

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