Discrete Fourier transform

DFT.R

# Simulate timeseries with period 10
time <- 1:200
y <- sin(time * 2 * pi / 10)


FT <- rep(NA, length(y) - 1)

for(k in 1:(length(y) - 1)){
  FT[k] <- sum(y * (cos(-2 * pi * k * 0:(length(y) - 1) / length(y)) + 1i * sin(-2 * pi * k * 0:(length(y) - 1) / length(y))))
}


amp <- Re(FT)

phase <- Im(FT)

plot(amp)




# Add some noise this time

y <- sin(time * 2 * pi / 10) + rnorm(200, 0, 0.4)


FT <- rep(NA, length(y) - 1)

for(k in 1:(length(y) - 1)){
  FT[k] <- sum(y * (cos(-2 * pi * k * 0:(length(y) - 1) / length(y)) + 1i * sin(-2 * pi * k * 0:(length(y) - 1) / length(y))))
}


amp <- Re(FT)

phase <- Im(FT)

plot(amp)










# Do some sine waves without noise

time <- 1:200
y <- sin(time * 2 * pi / 10) 

# Make vectors, but the first two elements will stay as NA. 1st one doesn't work, second one just gives a weird answer.
amp <- rep(NA, length(y) - 1)
phase <- rep(NA, length(y) - 1)

for(k in 3:(length(y) - 1)){
  amp[k] <- (lm(y ~ sin(time * 2 * pi / k) + cos(time * 2 * pi / k)))$coefficients[2]
  phase[k] <- (lm(y ~ sin(time * 2 * pi / k) + cos(time * 2 * pi / k)))$coefficients[3]
}

plot(amp)



# Do period 20 with missing data this time.

time <- 1:200
y <- sin(time * 2 * pi / 20) 

y[sample(1:200, 8)] <- NA



amp <- rep(NA, length(y) - 1)
phase <- rep(NA, length(y) - 1)

for(k in 3:(length(y) - 1)){
  amp[k] <- (lm(y ~ sin(time * 2 * pi / k) + cos(time * 2 * pi / k)))$coefficients[2]
  phase[k] <- (lm(y ~ sin(time * 2 * pi / k) + cos(time * 2 * pi / k)))$coefficients[3]
}

#Still works fine.
plot(amp)