## Adapted from kzft::kzp example 2
t <- 1:2000
y <- 1.1*sin(2*pi*0.0339*t)+7*sin(2*pi*0.0366*t)+5*rnorm(length(t),0,1)
y[sample(t,100,replace=FALSE)] <- NA
ft <- kz.ft(y, f=c(0.0339, 0.0366), k=2, m=1000, n=10)
# It may take 10 ~ 20 seconds
# system.time(ft <- kz.ft(y, k=2, m=1000, n=10))
plot(y=log(ft$pg+1), x=ft$f, type="l", xlim=c(0.02,0.05))
abline(v=c(0.0339, 0.0366), lty=21, col="red")
## recover signal
t <- 1:2000
y <- 1.1*sin(2*pi*0.01*t)+2*sin(2*pi*0.03*t)
y[sample(t,500,replace=FALSE)] <- NA
noise <- rnorm(2000,0,1)
system.time(ft <- kz.ft(y + 1.0*noise, k=5, f=c(0.01,0.03), m=100))
yr <- 2*Re(rowSums(ft$tf))
cor(yr, y[1:length(yr)], use="pairwise.complete.obs")
plot(yr, type="p", cex=0.5)
points(y[1:length(yr)],type="l", col="red")
## Example for kz.ftc
t <- runif(2000)*2000
f <- c(0.15, 0.1)
x <- sin(2*pi*f[1]*t + pi/4)
y <- sin(2*pi*f[2]*t + pi/12)
y <- y[order(t)]
x <- x[order(t)]
tr <- t[order(t)]
noise <- rnorm(length(tr),0,1)
plot(y=y+x, x=tr, type="l")
ft <- kz.ftc(x+y+2*noise, xt=tr, k=2, m=1000)
plot(y=ft$pg, x=ft$f, type="l")
abline(v=f, col="grey", lty=21)
mtext("Spectrum of Longitudinal Data, Selected f")
ft <- kz.ftc(x+y+noise, xt=tr, f=f, k=1, m=1900)
yr <- rowSums(2*Re(ft$tf))
iv <- 0:60
plot(y=(x+y+noise)[iv], x=tr[iv], type="p", col="grey")
xt <- (0:8000)/100
yt <- sin(2*pi*f[1]*xt+pi/4) + sin(2*pi*f[2]*xt+pi/12)
y2 <- sin(2*pi*f[1]*iv+pi/4) + sin(2*pi*f[2]*iv+pi/12)
points(yt, x=xt, col="grey", cex=0.5, lwd=1, type="l")
points(y2, x=iv, col="blue", cex=0.75, lwd=1, type="p")
points(y=yr, x=0:(length(yr)-1), type="p", cex=0.5, lwd=1, col="red")
mtext("k=1, m=1900, x+y+noise, red for reconstruction, grey points for data input", cex=0.75)
Run the code above in your browser using DataLab