FFTC

Returns a three-column array containing the complex form of the data array's Discrete Fourier Transform (DFT), computed using the radix-2 FFT (Fast Fourier Transform) algorithm.

Harmonics are listed in ascending order, with the zeroth harmonic in the first row. The number of rows is N div 2 + 1, where N is the nearest power of 2 equal to or higher than the number of elements in the data array.

For each harmonic, column #1 contains the real part of the Fourier coefficient, column #2 contains the imaginary part and column #3 contains the harmonic's angular frequency (listed for convenience, although implied by the row number and not a part of the transform proper).

All elements in data are converted to floats (with exclusion if conversion fails) and interpreted as a single, one-dimensional sequence. If the data array is two-dimensional, it's read in the usual order (i.e. row by row: left to right, top to bottom).

The radix-2 FFT algorithm requires the number of input elements to be a power of 2. In order to accomodate other data sizes, the fftc function uses zero padding: the sequence length is incremented up to the nearest power of 2 by appending zeros after the actual data. For instance, a 100-element array would be padded with 28 zeros to create an input sequence containing 2^7 = 128 elements. A 128-element array would not be padded since its length is already a power of 2.

The Discrete Fourier Transform is a decomposition of the input sequence as a sum of complex exponential sequences with frequencies that are integer multiples of the input sequence's fundamental frequency, 2 * PI / N (the input sequence is implicitly assumed to repeat itself every N elements). The zeroth harmonic has zero frequency (it's a constant), the first harmonic has the fundamental frequency (2 * PI / N), the second harmonic has twice the fundamental frequency (4 * PI / N) and so on.

Mathematically, the input sequence x[n] (with n in the range [0, N - 1]) is decomposed as the sum

		N - 1		i (2 p / N) k n
x[n] =	1	å	c[k] e
	¾
	N
		k = 0

where the complex number c[k] is the k:th Fourier coefficient, e is Euler's number and i is the imaginary unit. When the input sequence x[n] is real, only the first N div 2 + 1 coefficients are unique: coefficients number N div 2 + 1 to N - 1 are complex conjugates (same real part, opposite sign for imaginary part) of coefficients N div 2 down to 1.

Since the fftc function works with real input values, it only needs to return the first N div 2 + 1 coefficients. In terms of the above expression, the first column of the result array contains Re(c[row number - 1]), the second column contain Im(c[row number - 1]) and the third column contains the angular frequency = (row number - 1) * 2 * PI / N.

The 100-point sine wave with wavelength 100, amplitude 1 and DC component (i.e. constant vertical offset) 10,

which equals 7.8125. This is not 10 (the DC component of the input data) because of zero padding: as far as the fftc function is concerned, the 100 input elements containing the sine wave are followed by 28 elements containing zeros, so the average distance from the x axis is actually (10 * 100 + 0 * 128) / 128 = 7.8125.

Used as graph sources and plotted together, the above expressions look like this:

The black dots are the original sine wave, the blue line is the real part of the transform and the red line is the imaginary part. Note how the transform is smeared out over many harmonics. The presence of strong periodicity in the input sequence doesn't exactly leap out at the viewer. This happens because a wave with period 100 is being decomposed as a sum of harmonics of the fundamental period 128: no single harmonic is a good fit to the input sequence, so many harmonics are needed to do the job.

yields a very different result: now the amplitude of the zeroth harmonic is 10, since the DC component of the input sequence is not diluted by zero padding, and the only other significant element in the result is the imaginary component of the first harmonic, which equals -0.5. Graphically,

The lesson to take away from this is that even strong periodicities in the data are easy to miss if their period is close to the length of the input sequence being analyzed.

To see this clearly, go back to a sine period of 100 but use a 1024-point sequence instead:

_data = mop("sin()", makevector(1024, 0, 2 * PI / 100)) + makevector(1024, 10, 0)

Now the amplitude of the zeroth harmonic is 10.014, and there is a strong peak around the tenth harmonic. Excluding the input sequence and zooming in on the peak, the graph of the DFT now looks like this:

Keep in mind that the last example is highly idealized, with a perfect sine wave running through more than ten full periods and no background noise. In real life data, periodicities tend to have a limited lifespan and to play out against a noisy background. Consider trying the Maximum Entropy Method (functions mem, memdom, memf, memp) when looking for periodicities in such data.