From: Wolfgang Waser <wolfgang.waser_at_rz.hu-berlin.de>

Date: Tue 30 Aug 2005 - 19:35:04 EST

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

# R-script

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 lengthResult:

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