####多个频率成分
x1 = 2*cos(2*pi*1:1000*6/100) + 3*sin(2*pi*1:1000*6/100)
x2 = 4*cos(2*pi*1:1000*10/100) + 5*sin(2*pi*1:1000*10/100)
x3 = 6*cos(2*pi*1:1000*40/100) + 7*sin(2*pi*1:1000*40/100)
x = x1 + x2 + x3
par(mfrow=c(2,2))
plot.ts(x1[1:100], ylim=c(-10,10), main=expression(omega==6/100~~~A^2==13))
plot.ts(x2[1:100], ylim=c(-10,10), main=expression(omega==10/100~~~A^2==41))
plot.ts(x3[1:100], ylim=c(-10,10), main=expression(omega==40/100~~~A^2==85))
plot.ts(x[1:100], ylim=c(-16,16), main="sum")

spectrum(x,span=5)
abline(v=(6/100))
abline(v=10/100)
abline(v=40/100)

par(mfrow=c(2,2))
P = abs(2*fft(x)/1000)^2; Fr = 0:999/1000
plot(Fr, P, type="o", xlab="frequency", ylab="periodogram",sub="x")
P = abs(2*fft(x1)/1000)^2; Fr = 0:999/1000
plot(Fr, P, type="o", xlab="frequency", ylab="periodogram",sub="x1")
P = abs(2*fft(x2)/1000)^2; Fr = 0:999/1000
plot(Fr, P, type="o", xlab="frequency", ylab="periodogram",sub="x2")
P = abs(2*fft(x3)/1000)^2; Fr = 0:999/1000
plot(Fr, P, type="o", xlab="frequency", ylab="periodogram",sub="x3")



fig.dir = ""
pdf(paste(fig.dir,"PeriodicSum_Obs.pdf",sep=""))
par(mfrow=c(2,1))
plot.ts(x, ylim=c(-16,16), main="x",ylab="Obs")
abline(v=500,col="blue",lwd=2);abline(v=600,col="blue",lwd=2)
plot.ts(x[500:600], ylim=c(-16,16), main="x",ylab="Obs")
dev.off()

pdf(paste(fig.dir,"PeriodicSum_FreqFFT.pdf",sep=""))
P = abs(2*fft(x)/1000)^2; Fr = 0:999/1000
plot(Fr, P, type="o", xlab="frequency", ylab="periodogram",sub="x")
dev.off()

pdf(paste(fig.dir,"PeriodicSum_Obs_Decompose.pdf",sep=""))
par(mfrow=c(2,2))
plot.ts(x1[500:600], ylim=c(-10,10), main=expression(omega==6/100~~~A^2==13))
plot.ts(x2[500:600], ylim=c(-10,10), main=expression(omega==10/100~~~A^2==41))
plot.ts(x3[500:600], ylim=c(-10,10), main=expression(omega==40/100~~~A^2==85))
plot.ts(x[500:600], ylim=c(-16,16), main="sum")
dev.off()

pdf(paste(fig.dir,"PeriodicSum_FreqFFT_Decompose.pdf",sep=""))
par(mfrow=c(2,2))
P = abs(2*fft(x1)/1000)^2; Fr = 0:999/1000
plot(Fr, P, type="o", xlab="frequency", ylab="periodogram",sub="x1")
P = abs(2*fft(x2)/1000)^2; Fr = 0:999/1000
plot(Fr, P, type="o", xlab="frequency", ylab="periodogram",sub="x2")
P = abs(2*fft(x3)/1000)^2; Fr = 0:999/1000
plot(Fr, P, type="o", xlab="frequency", ylab="periodogram",sub="x3")
P = abs(2*fft(x)/1000)^2; Fr = 0:999/1000
plot(Fr, P, type="o", xlab="frequency", ylab="periodogram",sub="x")
dev.off()

pdf(paste(fig.dir,"PeriodicSum_Spectral.pdf",sep=""))
spectrum(x,span=5)
abline(v=(6/100))
abline(v=10/100)
abline(v=40/100)
dev.off()




###imotor data
www <- "http://www.math.pku.edu.cn/teachers/xirb/Courses/TimeSeries/IntroTimeSeriesRData/imotor.txt"
imotor.dat <- read.table(www, header = T)
attach (imotor.dat)

pdf(paste(fig.dir,"imotor_obs.pdf",sep=""))
layout (1:2)
plot(good[4400:5600],type="l",ylab="good")
plot(broken[4400:5600],type="l",ylab="bad",col="red")
dev.off()

pdf(paste(fig.dir,"imotor_spectral.pdf",sep=""))
layout (1:3)
xg.spec <- spectrum(good, span = 9)
xb.spec <- spectrum(broken, span = 9,col="red")
freqg <- 400 * xg.spec$freq [4400:5600]
freqb <- 400 * xb.spec$freq [4400:5600]
plot(freqg, 10*log10(xg.spec$spec[4400:5600]), main = "",
xlab = "Frequency (Hz)", ylab = "Current spectrum (dB)", type="l")
lines(freqb, 10 * log10(xb.spec$spec[4400:5600]), lty = "dashed",col="red")
dev.off()



####generate AR(1)


set.seed(1)
x <- w <- rnorm(1024,sd=0.5)
x[1] = 1
for (t in 2:1024) x[t]<- 0.9 * x[t-1] + w[t]
plot(as.ts(x))



set.seed(1)
x <- w <- rnorm(1024,sd=0.5)
x[1] = 1
for (t in 2:1024) x[t]<- 1 * x[t-1] + w[t]
plot(as.ts(x))


set.seed(1)
x <- w <- rnorm(1024,sd=0.5)
x[1] = 1
for (t in 2:1024) x[t]<- 1.1 * x[t-1] + w[t]
plot(as.ts(x))