5.8 Fourier transforms

The fft and ifft commands take Fourier transforms and inverse Fourier transforms respectively of data. As with other commands, data can be supplied from a data file, from functions, or from a colon-separated list of vectors (see Section 6.5.3). In each case, a regular grid of abscissa values must be specified on which to take the discrete Fourier transform, which can extend over an arbitrary number of dimensions. The following example demonstrates the syntax of these commands as applied to a two-dimensional top-hat function:

step(x,y) = tophat(x,0.2) * tophat(y,0.4)
fft  [  0: 1:0.01][  0: 1:0.01] f() of step()
ifft [-50:49:1   ][-50:49:1   ] g() of f()

In the fft command above, $N_ x=100$ equally-spaced samples of the function step $(x,y)$ are taken between limits of $x_\mathrm {min}=0$ and $x_\mathrm {max}=1$ for each of $N_ y=100$ equally-spaced values of $y$ on an identical raster, giving a total of $10^4$ samples. These are converted into a rectangular grid of $10^4$ samples of the Fourier transform f $(\omega _ x,\omega _ y)$ at

	$\displaystyle \omega _ x = \frac{j}{\Delta x}$	$\displaystyle \textrm{for}$	$\displaystyle -\frac{N_ x}{2}\leq j <\frac{N_ x}{2} \; \left(\textrm{equivalently, for} -\frac{N_ x}{2\Delta x}\leq \omega _ x <\frac{N_ x}{2\Delta x} \right), \nonumber$
	$\displaystyle \omega _ y = \frac{k}{\Delta y}$	$\displaystyle \textrm{for}$	$\displaystyle -\frac{N_ y}{2}\leq k <\frac{N_ y}{2} \; \left(\textrm{equivalently, for} -\frac{N_ y}{2\Delta y}\leq \omega _ y <\frac{N_ y}{2\Delta y} \right). \nonumber$

where $\Delta x=x_\mathrm {max}-x_\mathrm {min}$ and $\Delta y$ is analogously defined. These samples are interpolated stepwise, such that an evaluation of the function f $(\omega _ x,\omega _ y)$ for general inputs yields the nearest sample, or zero outside the rectangular grid where samples are available. In general, even the Fourier transforms of real functions are complex, and their evaluation when complex arithmetic is not enabled (see Section 4.5) is likely to fail. For this reason, a warning is issued if complex arithmetic is disabled when a Fourier transform function is evaluated.

In the example above, we go on to convert this set of samples back into the function with which we started by instructing the ifft command to take $N_ x=100$ equally-spaced samples along the $\omega _ x$ -axis between $\omega _{x,\mathrm{min}}=-{N_ x}/{2\Delta x}$ and $\omega _{x,\mathrm{max}}=(N_ x-1)/{2\Delta x}$ , with similar sampling along the $\omega _ y$ -axis.

Taking the simpler example of a one-dimensional Fourier transform for clarity, as might be calculated by the instructions

step(x) = tophat(x,0.2)
fft  [  0: 1:0.01] f() of step()

the fft and ifft commands compute, respectively, discrete sets of samples $F_ m$ and $I_ n$ of the functions $F(\omega _ x)$ and $I(x)$ , which are given by

$F_ j = \sum _{m=0}^{N-1} I_ m e^{-2\pi ijm/N} \, \delta x,\; \textrm{for}\; -\frac{N}{2}\leq j <\frac{N}{2} ,$

and

$I_ j = \sum _{m=0}^{N-1} F_ m e^{ 2\pi ijm/N} \, \delta \omega _ x,\; \textrm{for}\; -\frac{N}{2}\leq j <\frac{N}{2} ,$

where:

$I(x)$	=	Function being Fourier transformed.
$F(\omega _ x)$	=	Fourier transform of $I()$ .
$N$	=	The number of values sampled along the abscissa axis.
$\delta x$	=	Spacing of values sampled along the abscissa axis.
$\delta \omega _ x$	=	Spacing of abscissa values sampled along the $\omega _ x$ axis.
$i$	=	$\sqrt {-1}$ .

It may be shown in the limit that $\delta x$ becomes small – i.e. when the number of samples taken becomes very large – that these sums approximate the integrals

$\begin{equation} F(\omega _ x) = \int I(x) e^{-2\pi ix\omega _ x} \, \mathrm{d}x , \end{equation}$

(5.1)

and

$\begin{equation} I(x) = \int F(\omega _ x) e^{ 2\pi ix\omega _ x} \, \mathrm{d}\omega _ x . \end{equation}$

(5.2)

Fourier transforms may also be taken of data stored in data files using syntax of the form

fft [-10:10:0.1] f() of 'datafile.dat'

In such cases, the data read from the data file for an $N$ -dimensional FFT must be arranged in $N+1$ columns¹, with the first $N$ containing the abscissa values for each of the $N$ dimensions, and the final column containing the data to be Fourier transformed. The abscissa values must strictly match those in the raster specified in the fft or ifft command, and must be arranged strictly in row-major order.

Example: The Fourier transform of a top-hat function

It is straightforward to show that the Fourier transform of a top-hat function of unit width is the function ${\rm sinc}(x^\prime =\pi x)=\sin (x^\prime )/x^\prime$ . If

$f(x)=\left\{ \begin{array}{l}1\; |x|\leq \nicefrac {1}{2}\\ 0\; |x|>\nicefrac {1}{2}\end{array}\right. ,$

the Fourier transform $F(\omega )$ of $f(x)$ is

	$\displaystyle F(\omega )$	$\displaystyle =$	$\displaystyle \int _0^\infty f(x) \exp \left(-2\pi ix\omega \right) \, \mathrm{d}x = \int _{-\nicefrac {1}{2}}^{\nicefrac {1}{2}} \exp \left(-2\pi ix\omega \right) \, \mathrm{d}x$
	$\displaystyle$	$\displaystyle =$	$\displaystyle \frac{1}{2\pi \omega }\left[ \exp \left(\pi i\omega \right) - \exp \left(-\pi i\omega \right) \right] = \frac{\sin (\pi \omega )}{\pi \omega } = {\rm sinc}(\pi \omega ).$

In this example, we demonstrate this numerically by taking the Fourier transform of such a step function, and comparing the result against the function sinc(x) which is pre-defined within Pyxplot:
set numerics complex
step(x) = tophat(x,0.5)
fft [-1:1:0.01] f() of step()
plot [-10:10] Re(f(x)), sinc(pi*x)
Note that the function Re(x) is needed in the plot statement here, since although the Fourier transform of a symmetric function is in theory real, in practice any numerical Fourier transform will yield a small imaginary component at the level of the accuracy of the numerical method used. Although the calculated numerical Fourier transform is defined throughout the range $-50\leq x<50$ , discretised with steps of size $\Updelta x=0.5$ , we only plot the central region in order to show clearly the stepping of the function:
$\includegraphics{examples/eps/ex_fft}$
In the following steps, we take the square of the function ${\rm sinc}(\pi x)$ just calculated, and then plot the numerical inverse Fourier transform of the result:
g(x) = f(x)**2
ifft [-50:49.5:0.5] h(x) of g(x)
plot [-2:2] Re(h(x))
$\includegraphics{examples/eps/ex_fft2}$
As can be seen, the result is a triangle function. This is the result which would be expected from the convolution theorem, which states that when the Fourier transforms of two functions are multiplied together and then inverse transformed, the result is the convolution of the two original functions. The convolution of a top-hat function with itself is, indeed, a triangle function.

Footnotes

The using, every, index and select modifiers can be used to arrange it into this form.