FFT

The Fast Fourier Transform or FFT, is the efficient algorithm for the numerical evaluation of the Discrete Fouler Tansform defined with
 * \(\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\displaystyle (1) ~ ~ ~

B_k=\sum_{m=0}^{N-1} A_m \exp\Big(-\mathrm{i} \frac{2 \pi}{N}~ k~ m\Big)\)

where \(N\) is natural number; usually, \(2^n\) for some natural \(n\); \(A\) is the input array of length \(N\), numbered beginning from 0; \(B\) is the output array of length \(N\), numbered beginning from 0. In such a way, \(k\) is integer number, \(0\le k < N\).

There are many ways of evaluation of the expression (1). The trivial one would requires of order of \(N^2\) operations. The efficient way requires of order on \(N~ \log_2 N\) operations and called FFT. There are several efficient algorithms to perform the FFT .

One of the simplest algorithms realized in the implementation of the C++ routine fft(*complex(double), int, int) is stored in file fafo.cin. This routine is used in the implementation fafo(*complex(double), int, int) of the Fourier operator.

FFT and fafo
The routine "fafo" is used in calculations described in articles and for plotting of some pictures in TORI.

Since year 2011, fafo.cin is available with examples. This routine seems to be simplest and the most portable among those available in the Web; it is less than one screen of code, and it seems to be compatible with various operational systems, the only C++ is required. The reason of use of the fafo.cin is not the performance (perhaps some factor of order of 0.5 can be achieved at the optimization of the CPU time), but its compactness, readability and portability; no any "installing" is necessary after the downloading.

Operator FFT defined with (1) is often confused with the implementation of the Fourier operator; both replace the input array to the output one. But the result is not the same. The simple way to see the difference is to perform the transform of the array of values of the Gaussian exponential, let
 * \(\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\displaystyle (2) ~ ~ ~ A_k=\exp(- x_k^2/2~ )\)
 * \(\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\displaystyle (3) ~ ~ ~ x_k= (k-N/2) d\)

where \(d=\sqrt{2 \pi/N}\) is step of the equidistant grid. It is assumed that \(k\) is integer; \(0\!\le\!k\!<\!N\). The expression (2) defines the discrete analogy of the simplest self-Fourier function. With such an input array, The fafo returns almost the same array as that at the input (at least, beginning with \(N=8\); the FFT returns some set of values that are difficult to interpret in terms of the Fourier operator.

The case of the Gaussian self-Fourier function by (2) is shown at the figure at right for \(N\!=\!16\). The dotted curve shows \(y=\exp(-x^2/2)\); this is initial exponential. The blue segmented line shows its discrete representation by (2). This representation practically coincide with the numerical estimate for its Fourier transform The Red curve shows the Discrete Fourier Transform, that is equivalent to equation (1).

If for the Gaussian exponential, the routine fft returns the array similar to that shown at the figure with red line, then this indicates that, perhaps, the routine works correctly and returns the approximation for the sum by equaiton (1).

The fft is used for the implementation of fafo; and it can be easy extracted from there. The algorithm of expression of fafo through the FFT is realized through the displacement of the elements of array. This can be done by various ways; the most straightforward (as simplest to understand) is realized.

Matlab and other languages
For the translation of fafo to other languages (for example, Matlab), where FFT is already implemented, there is no need to translate FFT routine. In particular, the conventional Matlabs's fft can be used "as is".

The numeric implementation of the fafo through the matlab built-in function fft can be even simplified using the built-in functions "fftshift" and "ifftshift"; For the array X, the correct call may have form
 * ifftshift(fftshift(X))

The way of realization of these functions may be dependent of the matlab version; one can play with these functions and experimentally find the way of calling, that approximates the Fourier operator
 * \(\!\!\!\!\!\!\!\!\!\!\!\!\! (4) ~ ~ ~ \displaystyle \mathcal{F}(A)(x)= \frac{1}{\sqrt{2\pi}} \int_{- \infty}^{\infty} \exp( - \mathrm{i} x y) A(y) ~ \mathrm{d} y \)