## STAT 720 sp 2019 Lec 07 Spectral

# The discrete Fourier transform (DFT)

n <- 100
t <- 1:n
X <- sin(2*pi*t*(1/n)) + cos(2*pi*t*(2/n)) + rnorm(n,0,.5)
plot(X~t)

lambda <- (-floor((n-1)/2):floor(n/2))/n*2*pi

E <- matrix(NA,n,n)
for(i in 1:n)
{
	E[i,] <- 1/sqrt(n) * exp(1i*i*lambda)
}		

# see that E is orthonormal
t(Conj(E)) %*% E

# compute discrete Fourier transform of X		
D <- as.vector(t(Conj(E)) %*% X)

# see that reconstructed series is same as original
as.vector( Re(E %*% D) )
X

# plot the periodogram, which is the squared modulus of D
I <- Mod(D[lambda>=0])^2
plot(I~lambda[lambda>=0])

# keep K frequencies:
K <- 2
j.keep <- order(-Mod(D)^2)[1:(2*K)]
X.approx <- Re(E[,c(j.keep)] %*% D[c(j.keep)])

plot(X~t)
lines(X.approx~t)