Fft.R

Code snipits and programs written in R, S or S-PLUS Fs <- 5000 # Sampling frequency [Hz] Ts <- 20 # Length of sample [s] f1 <- 2 # Frequency of signal [Hz] f2 <- 2.7

A1 <- 1.3 # Amplitude of sine functions A2 <- 0.5 Time <- seq(0, Ts, length=Fs*Ts)

y <- A1*sin(2*pi*f1*Time) + A2*sin(2*pi*f2*Time+20*pi/180) + 0.01 * rnorm(length(Time)) par(mfrow=c(1,1)) plot(Time, y, type="l")

fty <- fft(y)/length(y)*2

Frequency <- Fs/2*seq(0,2, length=length(fty))

par(mfrow=c(2,1))

xmin <- 1.9 xmax <- 2.8

tolerance <- 0.005

1. Phasors with an amplitude of zero make it hard to determine phase

ok <- (Mod(fty) >= tolerance) Colours <- ifelse(ok, "black","red")

plot(Frequency, Mod(fty), type="b", xlim=c(xmin, xmax)) # Close up of peaks grid(NULL, NULL, lwd = 2) points(Frequency, Mod(fty), col=Colours) plot(Frequency, Arg(fty), type="b", xlim=c(xmin, xmax)) # Close up of peaks grid(NULL, NULL, lwd = 2) points(Frequency, Arg(fty), col=Colours)

1. See also Nyquist frequency http://en.wikipedia.org/wiki/Nyquist_frequency

plot(Frequency, Mod(fty), type="l") plot(Frequency, Arg(fty), type="l")

X11() par(mfrow=c(2,1)) angles <- Arg(fty) angles[!ok] <- 0 plot(Frequency, Mod(fty), type="l") plot(Frequency, angles, type="l")

graphics.off()