[R] FFT Normalization Documentation

[Franklin Bretschneider]

Dear Elke; Jeff,


Re:

> Eike: Understanding Discrete Fourier Transform theory is not trivial... while a vignette added to the stat package has the potential help a lot of users, it is a bit ambitious to try to supplant the extensive published material on using and interpreting the DFT (particularly as there is "more than one way to do it" and the R fft() function is very typical of fft implementations). (Similar arguments could be applied to most of the stat package... note the absence of vignettes there.) It might be more practical to propose to R-devel some patches to the fft() help file references and examples sections. Alternatively, you could write YAB (Yet Another Blog) for people to search for.
> 
> Frank: While folding is an important concept to know about when interpreting DFT results, I think something went rather wrong in your example with your "mask" variable since folding applies to f (for forward fft) or t (for inverse fft), not to the corresponding magnitudes. In addition to that, it is simply not necessary to pre-fold your data before applying the fft... the folding is assumed by the math to exist in the input outside the input window, and there is nothing you can do to the data to affect that assumption. Folding in the output is more visibly evident, but presenting it as a symmetric plot is entirely optional and is not done in most cases.



Maybe I didn't use the proper terminology, but what I called 'folding' is a modification of the input signal used only to present the amplitude spectrum in a convenient way. The FFT ("butterfly algorithm") yields a complex array where the highest frequencies (pos and neg) are in the middle, the lowest (and DC and fNyq) are at the ends. To display this same array with the DC value in the middle, the neg frequencies increasing to the left and the pos frequencies to the right, the trick with the +1/-1 mask is performed. This mask function is in fact a "square wave" at the Nyquist frequency.
In Matlab, it is in a routine called "fftshift", see here:

> 
> Y = fftshift(X) rearranges the outputs of fft, fft2, and fftn by moving the zero-frequency component to the center of the array. It is useful for visualizing a Fourier transform with the zero-frequency component in the middle of the spectrum.
> 

This is from the MathWorks web site:     http://nl.mathworks.com/help/matlab/ref/fftshift.html.

In addition, in my example I forgot to scale the amplitude. This must indeed be divided by n (the number of data points).
So, change my line YY <- fft(yy) into YY <- fft(yy)/n. Now the amplitudes of the spectral line are numerically the same as given in the composition of y.
These values must indeed be regarded with caution, since with real-world signals the energy will most often be spread among several spectral "lines".
Windowing (Hann, Hanning, Blackman etc.) then improves the spectrum, but that's a different story.

Best wishes,


Frank
---




Franklin Bretschneider
Dept of Biology
Utrecht University
bretschr at xs4all.nl