library(jpeg)
implot <- function(x){
    xx <- t(x[nrow(x):1,ncol(x):1])
    xx = -xx
    xx <- xx-min(xx)
    xx <- xx/max(xx)
    image( xx,axes=FALSE)
}
    
ftrans <- function(x){
    N <- length(x)
    a <- fft(x)
    freq <- (0:(N/2))/N
    a <- a[ 1:length(freq)]
    list( freq= freq, ampl = abs( a), phase=Arg(a))
}

plotfourier <- function( x){
    x <- ftrans(x)
    plot( x$freq, x$ampl, type="h", xlab=" Frequency", ylab="Amplitude")
}

## X-ray pulsar
## data from "http://xweb.nrl.navy.mil/timeseries/herx1.diskette"
dat <- read.table("pulsar.txt",header=FALSE)
x <- as.vector(t(as.matrix(dat)))
plot(x[1:2000],type="l")
plot(abs(fft(x)))
plotfourier(x[1:200])






## simulated cartoon version of exoplanet/pulsar time series
n <- 10000
x <- 0.1*rnorm(n) + 0.1*sin(2*pi* (0:(n-1))*1/100)^100
plot(x[1:2000],type="l")
plotfourier(x) 






###white noise
n <- 1000
W <- rnorm(n)
plot(W,type="l")
plot(abs(fft(W))) 
plotfourier(W)








## MA-process
n = 1000
W = rnorm(n)
theta = c(1,1,0.8,0.5,0.3,0.2)
theta = c(0,0,1,0,0,0)
theta = c(1,-1)
X = W
for (i in 1:n){
    X[i] = sum(W[pmax(1,i-(0:5))] * theta)
}





par(mfrow=c(1,3))
plot(X)
acf(X)
plotfourier(X)

fftheta = fft(c(theta, rep(0, length(W)-length(theta))))
ffW = fft(W)
ffXc = ffW * fftheta

plot(Re(fft( ffXc, inv= TRUE))/n)
lines(X, type="l")






#### AR-process
x <- filter(W, filter= c(1.9,-0.995),method="recursive", init=rep(0,2) ) 
plot(x,type="l")
plot(abs(fft(x)))
plotfourier(x)






## timeshift
par(mfrow=c(3,2))
x <- filter(W, filter= c(1.9,-0.995),method="recursive", init=rep(0,2) ) 
plot(x)

z <- rep(0,length(x))
z[200] <- 1
plot(z)

ffx <- fft(x)
plot(abs(ffx),type="l")

ffz <- fft(z)
plot(abs(ffz))

ffxz <- ffx*ffz
xz <- Re(fft( ffxz, inverse=TRUE))
plot(x,type="l")
plot(xz,type="l")






## images
im <- readJPEG("tukey.jpg")
ffi <- fft(im)

ims <- 0*im
ims[100,100] <- 1
ffs <- fft(ims)
print(abs(ffs)) ## just a phase-shift

par(mfrow=c(1,2))
implot(im)
implot( Re(fft(ffi*ffs, inverse=TRUE)))




## mixing phase and amplitude
n <- 300
W <- rnorm(n)

x1 <- filter(W, filter= rep(1,5),sides=1,circular=TRUE ) ## MA-process
x2 <- filter(W, filter= c(1.99,-0.995),method="recursive", init=rep(0,2) )  ## AR-process

par(mfrow=c(4,2))
plot(x1,type="l")
plot(x2,type="l")

acf(x1)
acf(x2)

f1 <- fft(x1)
f2 <- fft(x2)
fn1 <- abs(f1)*exp(1i * Arg(f2))
fn2 <- abs(f2)*exp(1i * Arg(f1))

xn1 <- as.numeric(fft(fn1, inverse=TRUE))
xn2 <- as.numeric(fft(fn2, inverse=TRUE))
plot(xn1,type="l")
plot(xn2,type="l")

acf(xn1)
acf(xn2)