[R] FFT, frequs, magnitudes, phases

[Michael A. Miller]

>>>>> "Wolfgang" == Wolfgang Waser <wolfgang.waser at rz.hu-berlin.de> writes:

    > I would be most obliged for any comments and help.

Wolfgang,

I've used R's fft to filter ECG signals and will comment on your
commentary based on my experience.  First, as an easily
accessible reference, I suggest "The Scientist and Engineer's
Guide to Digital Signal Processing," which is available in pdf
form at http://www.dspguide.com.  It includes several chapters on
the discrete Fourier transform and the fast Fourier transform
algorithm (which is what R's fft implements) and a chapter on
applications that contains info on spectral analysis.

    > What does R's fft() deliver?

    > fft() is calculated with a single one-dimensional
    > vector. Information on data acquisition frequency and block
    > length (in sec or whatever) can not be included into the
    > fft()-call.

    > Confusingly, if you calculate fft() on a sample vector
    > consisting of 2 pure sin frequencies, you get 4 peaks, not
    > 2.

That is the nature of the fft algorithm.  It returns values of
the discrete Fourier transform for both positive and negative
frequencies.

    > As stated above, fft() gives only "meaningful" frequency up
    > to half the sampling frequency. R, however, gives you
    > frequencies up to the sampling frequency. 

It is important to remember that the fft algorithm doesn't return
any frequency data at all.  It returns values of the fft that
correspond to frequencies from -f_Nyquist to +f_Nyquist.  It is
up to the user to calculate the frequency values.

Here's an example:

## Read some sample ecg data
ecg <- read.table('http://www.indyrad.iupui.edu/public/mmiller3/sample-ecg-1kHz.txt')
names(ecg) <- c('t','ecg')

ecg$t <- ecg$t/1000  # convert from ms to s

par(mfrow=c(2,2))

## Plot the ecg:
plot(ecg ~ t, data=ecg, type='l', main='ECG data sampled at 1 kHz', xlab='Time [s]')

## Calculate fft(ecg):
ecg$fft <- fft(ecg$ecg)

## Plot fft(ecg):
#plot(ecg$fft, type='l')

## Plot Mod(fft(ecg)):
plot(Mod(ecg$fft), type='l', log='y', main='FFT of ecg vs index')

## Find the sample period:
delta <- ecg$t[2] - ecg$t[1]

## Calculate the Nyquist frequency:
f.Nyquist <- 1 / 2 / delta

## Calculate the frequencies.  (Since ecg$t is in seconds, delta
## is in seconds, f.Nyquist is in Hz and ecg$freq is in Hz)
## (Note: I may be off by 1 in indexing here ????)
ecg$freq <- f.Nyquist*c(seq(nrow(ecg)/2), -rev(seq(nrow(ecg)/2)))/(nrow(ecg)/2)

## Plot fft vs frequency
plot(Mod(fft) ~ freq, data=ecg, type='l', log='y', main='FFT of ECG vs frequency', xlab='Frequency [Hz]')

## Now let's look at some artificial data:
x <- seq(100000)/1000  # pretend we're sampling at 1 kHz

## We'll put in two frequency components, plus a dc offset
f1 <- 5  # Hz
f2 <- 2  # Hz
y <- 0.1*sin(2*pi*f1*x) + sin(2*pi*f2*x) + 50
fft.y <- fft(y)
delta <- x[2] - x[1]
f.Nyquist <- 1 / 2 / delta
f <- f.Nyquist*c(seq(length(x)/2), -rev(seq(length(x)/2)))/(length(x)/2)

par(mfrow=c(2,2))
plot(x,y, type='l', xlim=c(0,20))
plot(f, Mod(fft.y), type='l', log='y')

## Now let's zoom in and mark the points were I expect to see peaks:
plot(f, Mod(fft.y), type='l', log='y', xlim=c(-10,10))
rug(c(-f1, -f2, 0, f1, f2), col='red', side=3)




Hope this is helpful, Mike

-- 
Michael A. Miller                               mmiller3 at iupui.edu
  Imaging Sciences, Department of Radiology, IU School of Medicine