Difference between revisions of "Fft.R"

From Organic Design wiki
(FFT of sin waves in R)
 
(like matlab example)
 
Line 1: Line 1:
# {{R}}
+
Fs <- 1000  # {{R}} Sampling frequency [Hz]
Fs <- 5000  # Sampling frequency [Hz]
+
Ts <- 1/Fs  # Sample time
Ts <- 20  # Length of sample [s]
+
L <- 10    # Length of signal [s]
 +
 
 
f1  <- 2    # Frequency of signal [Hz]
 
f1  <- 2    # Frequency of signal [Hz]
 
f2 <- 2.7
 
f2 <- 2.7
  
A1  <- 1.3 # Amplitude of sine functions
+
A1  <- 2.7 # Amplitude of sine functions
 
A2  <- 0.5  
 
A2  <- 0.5  
Time <- seq(0, Ts, length=Fs*Ts)
+
Time <- seq(0, L, length=Fs*L)
 
 
  
 
y    <- A1*sin(2*pi*f1*Time) + A2*sin(2*pi*f2*Time+20*pi/180) + 0.01 * rnorm(length(Time))
 
y    <- A1*sin(2*pi*f1*Time) + A2*sin(2*pi*f2*Time+20*pi/180) + 0.01 * rnorm(length(Time))

Latest revision as of 21:59, 2 August 2009

Fs <- 1000 #

Code snipits and programs written in R, S or S-PLUS Sampling frequency [Hz] Ts <- 1/Fs # Sample time L <- 10 # Length of signal [s]

f1 <- 2 # Frequency of signal [Hz] f2 <- 2.7

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

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()