Dear Frank, thanks a lot for your reply. I was aware of the need for FFT normalization and folding, but you're making another good point here, stating that there should be examples on this (folding) in the documentation as it will pose quite an obstacle for newcomers attempting to practically use the FFT. When writing yesterday's mail, I falsely was under the impression that there is such a thing as a "correct" normalization of an FFT. In the meantime, I have come to realize that when considering "real", i.e., physical, signals, there is no "true" way of normalizing the result of an inverse FFT, since everything depends on the sampling characterstics and on the properties of the frequency domain representation of the signal. Just imagine two differently sampled versions of the frequency representation of some physical signal, both attaining the same number of sample points NFFT, but different sampling frequencies. There is no way to scale these two signals coherently, based only on the number of sample points. Assuming Fs1 < Fs2 and assuming the physical signal to be band-limited with Fband < Fs1, the inverse transform of the signal sampled at Fs1 will naturally be larger (due to the larger number of sinces/cosines being summed up in the time signal). When comparing differently sampled versions of the same signal, it might be desirable to scale the result based on the energy content of the signal, i.e. by multiplying the signal with the frequency bin width 'fbin = Fs/(NFFT-1)'. So this is where your calculation of the attained frequencies comes into play! Note that this is equivalent to calculating the energy in the frequency domain by integration of a piecewise constant interpolation of the frequency domain FFT coefficients. I want to repeat my question concerning proposals on where to publish such considerations for the use of as many R users as possible. Might it be a good idea to create and propose an official vignette to the FFT function? Kind regards, Eike
> On February 3, 2015 at 12:39 AM Franklin Bretschneider <brets...@xs4all.nl> > wrote: > > > Dear Eike Petersen, > > > Re: > > > Hello everyone, > > > > the docpage for the fft function states: > > > > “Description: Performs the Fast Fourier Transform of an array.” > > > > and > > > > “Arguments – inverse: if ‘TRUE’, the unnormalized inverse transform is > > computed > > (the inverse has a ‘+’ in the exponent of e, but here, we do _not_ divide by > > ‘1/length(x)’).” > > > > Judging from this, I would expect ‘fft(x)’ to yield the correct FFT of x, > > and > > ‘fft(X, inverse = TRUE) / length(X)’ to yield the correct inverse FFT of X. > > > > However, it seems to me that actually the result of ‘fft(x)’ should be > > scaled by > > ‘1/length(x)’, while ‘fft(X, inverse=TRUE)’ seems to yield a correct result > > by > > default: > > > > t <- seq(0.001, 1, 0.001) > > y <- 1 + sin(2*pi*20*t) + 1.5 * sin(2*pi*30*t) > > Y <- fft(y) > > dev.new() > > plot(abs(Y)) ## Shows peaks at amplitudes 1000, 500 and 750, while they > > should > > be at amplitude 1, 0.5 and 0.75. > > y2 <- Re(fft(Y / length(Y), inverse = TRUE)) > > max(abs(y-y2)) ## The IFFT yields a correctly scaled result by default, if > > applied to a correctly scaled FFT. > > > > Did I get something wrong? If not, having spent quite some time figuring > > this > > out, I would like to see the documentation clearly pointing this out. I find > > the > > current text rather confusing. > > > > On another note: I have spent some time working on demo files that showcase > > some > > of the properties of the FFT and their implementation in R. I have done this > > primarily for myself, as I keep forgetting how these things work, but I > > thought > > that it might be helpful to others as well. Any hints on where/how I should > > publish such a thing? > > > > Kind regards and many thanks in advance, > > > > Eike > > > > ______________________________________________ > > R-help@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. > > > As far as I know an FFT must be normalized and folded to obtain a spectrum in > the form we like, so this would be my version: > > # your time signal: > > t <- seq(0.001, 1, 0.001) > y <- 1 + sin(2*pi*20*t) + 1.5 * sin(2*pi*30*t) > > # find time and frequency calibration: > > n = length(y) > dt = t[2]-t[1] > fNyq = 1/(2*dt) > tmax = max(t) > df = 1/tmax > > # make frequency vector to display as x-values of the spectrum rather than the > index. > > f = seq(-fNyq, fNyq-df, by=df) > > # make folding mask > > mask=rep(c(1, -1),length.out=n) > > # fold the spectrum around the Nyquist frequency; so the DC value (f=0) is in > the middle; the - and + max frequency at the ends. > > yy = y * mask > > # # Then do the FFT > > YY <- fft(yy) > > Plot the amplitude spectrum vector against the freq. vector > > plot(f,abs(YY), type='h') > > > > -------- > > It would be a good idea to put such an example in the help pages indeed. > The short example given in the manual isn't of much help for the usual > (periodic!) time signals. > > Hoping this helps, I remain > > With best wishes, > > Frank > ---- > > > > > > > Franklin Bretschneider > Dept of Biology > Utrecht University > brets...@xs4all.nl > > > [[alternative HTML version deleted]] ______________________________________________ R-help@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.