### Lecture 18 ################################################# # # Examples of DFT # ################################################################ reset <- function() { par(mfrow=c(1,1)); } # --- constant x <- rep(1,8) ft <- fft(x) par(mfrow=c(1,2)) plot(x); plot(abs(ft)); reset(); # --- slow ft (n is even) F <- function(n) { f <- matrix(0, nrow=n, ncol=n); i <- 0:(n-1) lim <- floor(n/2)-1 if (lim >= 1) { for (j in 1:(floor(n/2)-1)) { col <- 2 * j; f[,col] <- cos(2*pi*j*i/n); # not most efficient f[,col+1] <- sin(2*pi*j*i/n); } } # handle the ends (no sine terms) f[,1] <- 1; f[,n] <- c(1,-1); f } F(2) round( F(8) , digits=4) # only one sine/cos round( t(F(8)) %*% F(8), digits = 8) sft <- t(F(8)) %*% x # --- boxcar n <- 64 t <- 0:(n-1) f1 <- fft(y1 <- c(rep(1,n/2),rep(0,n/2))); f1; par(mfrow=c(1,2)) plot(y1); plot(abs(f1), type="b"); reset(); f2 <- fft(y2 <- c(rep(1,n/4),rep(0,n/2),rep(1,n/4)));f2; par(mfrow=c(1,2)) plot(y2); plot(abs(f2), type="b"); reset(); plot(abs(f2)); lines(abs(f1), col="blue") Mod(f1[1:5]) Mod(f2[1:5]) Arg(f1[1:5]) Arg(f2[1:5]) # --- sinusoid, location at FF n <- 64 t <- 0:(n-1) la <- 10/n f1 <- fft(y1 <- sin(2*pi*t*la)); f1 <- fft(y1 <- sin(2*pi*(t+24)*la)); plot(Mod(f1)) xi <- la + 1/(5*n) f2 <- fft(sin(2*pi*t*xi)); Mod(f2); lines(Mod(f2), col="lightblue") xi <- la + 1/(3*n) f3 <- fft(sin(2*pi*t*xi)); Mod(f3); lines(Mod(f3), col="pink") xi <- la + 1/(2*n) f4 <- fft(sin(2*pi*t*xi)); Mod(f4); lines(Mod(f4), col="black") # --- sinusoid + outlier n <- 64 t <- 0:(n-1) la <- 10/n y <- sin(2*pi*t*la) outlier <- rep(0,n); outlier[6] <- 1 points(y+outlier) plot(y, type="b"); points(y+outlier, col="red") f <- fft(y + outlier); plot(Mod(f)) f1 <- fft(y); plot(Mod(f1)) f2 <- fft(outlier); plot(Mod(f2)) # --- repeating pattern (like seasonal) plus noise n <- 128 t <- 0:(n-1) per<- 16 la <- 10.5/per x <- sin(2 * pi * la * 1:per) s <- 0 y <- rep(x,n/per) + s * rnorm(n) par(mfrow=c(1,2)) plot(y); plot(abs(fft(y)), type="b"); reset();