Re: [R] FFT, frequs, magnitudes, phases

From: Michael A. Miller <mmiller3_at_iupui.edu>
Date: Wed 31 Aug 2005 - 03:35:25 EST

>>>>> "Wolfgang" == Wolfgang Waser <wolfgang.waser@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@iupui.edu
  Imaging Sciences, Department of Radiology, IU School of Medicine

______________________________________________
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Received on Wed Aug 31 03:40:14 2005

This archive was generated by hypermail 2.1.8 : Sun 23 Oct 2005 - 16:09:08 EST