The Discrete Fourier Transform is used to convert an equally-space sampled signal from its representation as a function of time to its representation as a function of frequency. The DFT calculates \(n/2\) coefficients from \(n\) samples covering \(D\) seconds. These coefficients (complex numbers) are the magnitude and phase at frequencies \(\frac{0}{D}\)Hz, \(\frac{1}{D}\)Hz, ..., \(\frac{n/2-1}{D}\)Hz of cosine functions. An efficient algorithm to calculate these cofficients is the Fast Fourier Transform. This algorithm works only if the number of samples is a power of two.

Sampling of the function:

fun: the function to sample, one input, one output range: the sampling range (min time and max time, inclusive) nbSample: the number of sample, the first sample is taken at min time, the last sample is taken at (min time + (max time - min time) * (nbSample-1)/nbSample) SampleFunction(fun, range, nbSample): samples: vector of nbSample real values rangeFrom = [0, nbSample-1] rangeTo = [ range.min, range.min + (range.max - range.min) * (nbSample-1)/nbSample] for iSample in [0..nbSample[: samples[i] = fun(lerp(i, rangeFrom, rangeTo)) return samples

Discrete Fourier Transform:

samples: the samples value range: the sampling range (min time and max time, inclusive) nbSample: the number of sample, the first sample is taken at min time, the last sample is taken at (min time + (max time - min time) * (nbSample-1)/nbSample) DFT(samples, range, nbSample): x: vector coeffs: structure to memorise the nbSample complex coefficients and the sampling range FFT(x, nbSample, coeffs.vals) coeffs.range = range coeffs.nbSample = nbSample return coeffs FFT(x, n, y): if n = 1: y = x else xe, xo, ye, yo: vectors of n/2 complex values xe = [x[0], x[2], ..., x[(n-1)/2]] xo = [x[1], x[3], ..., x[(n-1)/2+1]] FFT(xe, n/2, ye) FFT(xo, n/2, yo) omega = exp(2.0 * PI * I / n) w = 1 for i in [0, n/2[: t = w * yo[i] y[i] = ye[i] + t y[i + n/2] = ye[i] - t; w *= omega

Spectrum of the function from its DFT coefficients:

coeffs: DFT coefficients as returned by the function DFT() above Spectrum(coeffs): bins: array of coeffs.nbSample/2 real values for i in [0, coeffs.nbSample/2[: bins[i] = sqrt(real_part(coeffs.vals[iBin])^2 + imaginary_part(coeffs.vals[iBin])^2) if i > 0: bins[i] *= 2 bins[iBin] /= coeffs.nbSample return bins

Reconstruction of the function from its DFT:

coeffs: DFT coefficients as returned by the function DFT() above t: real value, time at which the function is evaluated EvalFromDFT(coeffs, t): x = (t - coeffs.range.min) / (coeffs.range.max - coeffs.range.min); y = real_part(coeffs.vals[0]) / coeffs.nbSample loop iCoeff in [1, coeffs.nbSample / 2[: amp = 2.0 * sqrt(real_part(coeffs.vals[iCoeff])^2 + imaginary_part(coeffs.vals[iCoeff])^2) / coeffs.nbSample theta = atan( imaginary_part(coeffs.vals[iCoeff]) / real_part(coeffs.vals[iCoeff])) phi = iCoeff * x * 2.0 * PI - theta y += amp * cos(phi) return y

Example of a square wave function (in red) DFT converted and reconstructed by sampling from t=1s to t=5s using 128 samples. First graph shows the magnitude for each of the 64 frequency bins (spectrum), second graph shows the square wave function and its reconstruction (in blue) from the DFT coefficients.