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
I see now the trick... the square wave clarified it for me. It is indeed faster than re-arranging the data if putting the zero frequency in the middle of the data is your goal. Since I only do that for the purposes of teaching I probably won't be using it, but it may well be an interesting "trick" to put in the examples. On Tue, 3 Feb 2015, Franklin Bretschneider wrote:> 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 > > >--------------------------------------------------------------------------- Jeff Newmiller The ..... ..... Go Live... DCN:<jdnewmil at dcn.davis.ca.us> Basics: ##.#. ##.#. Live Go... Live: OO#.. Dead: OO#.. Playing Research Engineer (Solar/Batteries O.O#. #.O#. with /Software/Embedded Controllers) .OO#. .OO#. rocks...1k
Dear Jeff, Franklin, Eike and other interested parties, I agree that the help page for fft() deserves a bit more detail; its source is https://svn.r-project.org/R/trunk/src/library/stats/man/fft.Rd For instance, it has the comment %% %% Here, we should really have a nice \deqn{}{} giving the definition %% of the DFT ! %% and indeed, the example is minimal, depicting in which sense the fft(*, inverse=TRUE) is the inverse of fft(). Both a definition of the DFT (in LaTeX or directly as \deqn{}{}) which matches *our* fft(), and an extended commented example (ca 20 lines, rather than ca 100 lines!) would be welcome additions. If you are willing to donate a bit of time to the R project, we (the R users world!) will be happy for extended versions of the file fft.Rd either in public, here, or as 'wishlist' bug report with an attachment if you want -- or by e-mail to me. You may also chose to discuss an "optimal example" by e-mail between you, as you've already started to do so, implicitly. Martin Maechler, ETH Zurich>>>>> Jeff Newmiller <jdnewmil at dcn.davis.ca.us> >>>>> on Tue, 3 Feb 2015 19:43:25 -0800 writes:> I see now the trick... the square wave clarified it for > me. It is indeed faster than re-arranging the data if > putting the zero frequency in the middle of the data is > your goal. > Since I only do that for the purposes of teaching I > probably won't be using it, but it may well be an > interesting "trick" to put in the examples. > On Tue, 3 Feb 2015, Franklin Bretschneider wrote: >> 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 >> >> >> > --------------------------------------------------------------------------- > Jeff Newmiller The ..... ..... Go Live... > DCN:<jdnewmil at dcn.davis.ca.us> Basics: ##.#. ##.#. Live > Go... Live: OO#.. Dead: OO#.. Playing Research Engineer > (Solar/Batteries O.O#. #.O#. with /Software/Embedded > Controllers) .OO#. .OO#. rocks...1k > ______________________________________________ > R-help at r-project.org mailing list -- To UNSUBSCRIBE and > more, see https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide > http://www.R-project.org/posting-guide.html and provide > commented, minimal, self-contained, reproducible code.