R help - Feb 2015 - FFT Normalization Documentation

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 <bretschr at
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 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
>
>
>
	[[alternative HTML version deleted]]

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.

FWIW, my take on this example:

# increment of time
dT <- 0.001 # assumed here to be seconds
# number of samples
N <- 1000
# sampled time axis
t <- seq( 1, N ) * dT
# Nyquist frequency
Fnyq <- 0.5/dT
# increment of frequency
dF <- 1 / ( N * dT )
# sampled frequency axis
f <- seq( 0L, N/2-1 ) * dF
# example data, whole number of cycles within the input window!
y <- 1 + sin( 2 * pi * 20 * t ) + 1.5 * sin( 2 * pi * 30 * t )
# Fourier Series scaling
Y <- fft( y ) / length( y )
# double all sinusoidal components, don't look beyond Nyquist frequency
Ymag <- abs( c( Y[ 1 ], 2 * Y[ 2:(N/2) ] ) )
# plot Fourier Series coefficients, note correspondence with equation
plot( f, Ymag, type='h', ylab="Amplitude",
xlab="Frequency (Hz)" )
# compute inverse FFT (if original data was real, inverse will be real also)
yr <- fft( Y, inverse = TRUE )
# check real components ... differs by numerical noise
max( Re( yr ) - y )
# check imaginary components ... differs by numerical noise
max( Im( yr ) )

# now try including a frequency component that exceeds Fnyq
# leads to aliasing
y2 <- 1 + sin( 2 * pi * 20 * t ) + 1.5 * sin( 2 * pi * 970 * t )
# Fourier Series scaling
Y2 <- fft( y2 ) / length( y2 )
# double all sinusoidal components, don't look beyond Nyquist frequency
Y2mag <- abs( c( Y2[ 1 ], 2 * Y2[ 2:(N/2) ] ) )
# plot Fourier Series coefficients, note unfortunate similarity to
# previous plot
plot( f, Y2mag, type='h', ylab="Amplitude",
xlab="Frequency (Hz)" )
# moral is not to apply fft when significant amplitudes are present above 
# the Nyquist frequency (know your data)

On Tue, 3 Feb 2015, Eike Petersen wrote:
> 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 <bretschr at
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 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
>>
>>
>>
> 	[[alternative HTML version deleted]]
>
> ______________________________________________
> 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.
---------------------------------------------------------------------------
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 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