[R] FFT Normalization Documentation

[Franklin Bretschneider]

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 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.


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
bretschr at xs4all.nl