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

From: Wolfgang Waser <wolfgang.waser_at_rz.hu-berlin.de>
Date: Tue 30 Aug 2005 - 19:35:04 EST


Hi,

here is some info about the first part of my "homework", for those, who want to break down their signal (heart beat or whatever) into a collection of pure sin waves to analyse "main" frequency magnitudes and phases.

First some very un-mathematical "applied" theory:

If you sample a waveform signal (heart beat pressure pulses, ECG, doppler flow signals, etc.) with a certain data acquisition frequency, an fft of your data gives you the decomposition/breakdown of the waveform signal into a series of pure sin waves of different frequencies. Each sin wave in that list has: a) a certain "magnitude", i.e. a measure of how much that particular frequency participates in the generation of your signal, and b) a phase, i.e. the starting point of each sin wave.

Two characteristics of an fft have to be considered: 1) the highest meaningful sin wave frequency of your fft-analysis of the original waveform signal is half the data acquisition frequency (actually, R's fft gives you a list of frequencies up to the acquisition frequency, but you can only use the first half of it, see below) 2) the frequency resolution of your fft-analysis depends on the sampling time. The longer the sampling/analysis interval, the finer the resolution. Frequency resolution is actually 1 divided by sampling time (sec).

An example:

- some complicated waveform signal
- 1000 Hz data acquisition frequency (going on for hours)
- fft-analysis of data blocks of 1 sec length
Result:
- vector of frequencies from 1 to 500 Hz with a resolution of 1 Hz, corresponding vector of magnitudes (one for each frequency) and phases (dito).
You can now e.g. pick the frequency with the highest magnitude within that 1 sec block and continue the fft analysis in 1 sec blocks for the complete data set, analysing the time course of the "main" frequency of your waveform signal.

If you need higher frequency resolution, increase the block length. Analysis of a 5 sec block will give you a list of frequencies from 0.2 to 500 Hz with a resolution of 0.2 Hz. However, increasing analysis-block length decreases temporal resolution, i.e. "main" frequency are now calculated only every 5 sec and not 1 sec.

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.

R delivers a single one-dimensional vector of the same length as the data vector containing a list of imaginary numbers. To extract the "magnitudes" use Mod(fft()). The magnitudes can also be calculated using the formula: magnitude = square root (real * real + imaginary * imaginary) real: Re(fft()), imaginary: Im(fft())

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

As stated above, fft() gives only "meaningful" frequency up to half the sampling frequency. R, however, gives you frequencies up to the sampling frequency. The point is, that sampling a signal in discret time intervals causes aliasing problems. E.g. when sampling a 50 Hz sin wave and 950 Hz sin wave with 1000 Hz, the results will be identical. An fft can not distinguish between the two frequencies. Therefore, the sampling frequency should always be at least twice as high as the expected signal frequency. So for each actual frequency in the signal, fft() will give 2 peaks (one at the "actual" frequency and one at sampling frequency minus "actual" frequency), making the second half of the magnitude vector a mirror image of the first half.
As long as the sampling frequency was at least twice as high as the expected signal frequency, all "meaningful" information is contained in the the first half of the magnitude vector. A peak in the low frequency range might nevertheless still be caused by a high "noise" frequency.

The vector of magnitudes extraced so far only has an index an no associated frequencies.

To calculated the frequencies, simply take (or generate) the index vector (1 to length(magnitude vector) and divide by the length of the data block (in sec).

That's it for now. The second half of my "homework" will be delivered as soon as I understand what to make out of the phases given by R. I again would expect a vector of the same length as the magnitude vector with the phases (0 to 2*pi or -pi to +pi) of each frequency. However, I do not know yet what R calculates.
I would be most obliged for any comments and help.

Wolfgang



# R-script
acq.freq <- 4000       # data acquisition frequency (Hz)
sig1.freq <- 50           # frequency of 1st signal component (Hz)
sig2.freq <- 130        # frequency of 2nd signal component (Hz)
time <- 5                    # measuring time interval (s)

# vector of sampling time-points (s)
smpl.int <- (1:(time*acq.freq))/acq.freq

# data vector containing two frequencies (2nd frequ with phase shift) data <- sin(sig1.freq*smpl.int*2*pi)+sin(sig2.freq*smpl.int*2*pi+pi/2)

plot(data,type="l")

# calculate fft of data
test <- fft(data)

# extract magnitudes and phases
magn <- Mod(test) # sqrt(Re(test)*Re(test)+Im(test)*Im(test)) phase <- Arg(test) # atan(Im(test)/Re(test))

# select only first half of vectors
magn.1 <- magn[1:(length(magn)/2)]
#phase.1 <- Arg(test)[1:(length(test)/2)]

# plot various vectors

# plot magnitudes as analyses by R
x11()
plot(magn,type="l")

# plot first half of magnitude vector
x11()
plot(magn.1,type="l")

# generate x-axis with frequencies
x.axis <- 1:length(magn.1)/time

# plot magnitudes against frequencies
x11()
plot(x=x.axis,y=magn.1,type="l")

>>> Hi,
>>>
>>> I'm in dire need of a fast fourier transformation for me
>>> stupid biologist, i.e. I have a heartbeat signal and
>>> would like to decompose it into pure sin waves, getting
>>> three vectors, one containing the frequencies of the sin
>>> waves, one the magnitudes and one the phases (that's what
>>> I get from my data acquisition software's FFT function).
>>> I'd be very much obliged, if someone could point out
>>> which command would do the job in R.

>> fft(), but notice that it gives the complex
>> transform. You need to do a little homework to get at
>> the magnitude/phase values. (Basically, you just have to
>> take Mod() and Arg(), but there some conventions about
>> the frequencies and multipliers that one can get wrong).

> Once you've finished the "homework", others might be interested
> in your result... so it will be found in the future using > RSiteSearch().
-- 
Dr. Wolfgang Waser
Humbolt-Universitšt zu Berlin
Institute of Biology
Department of Animal Physiology
Philippstrasse 13, Abderhaldenhaus
10115 Berlin
Germany
Tel: +49 (0)30 2093 6173
Fax: +49 (0)30 2093 6375

______________________________________________
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 Tue Aug 30 19:46:10 2005

This archive was generated by hypermail 2.1.8 : Fri 03 Mar 2006 - 03:40:00 EST