library("Cairo") getQuasiPeriod <- function(To,Tsd,Ao,Asd,No,Nsd,tlen,savg) { w <- NULL T <- NULL A <- NULL wx <- NULL periods = (tlen/To) * 2 # for wiggle room for (i in 1:periods) { Ti = To + rnorm(1,sd=Tsd) idx = Ti-1 T = c(T,rep(Ti,idx)) Ai = Ao + rnorm(1,sd=Asd) A = c(A,rep(Ai,idx)) w = 1 wx = c(wx,rep((length(wx)+1),idx)) } x <- 1:length(wx) y = A*sin(((x-wx)/T)*(2 * pi)) + (No + rnorm(length(wx),sd=Nsd)) y = filter(y, rep(1/savg,savg), sides=2) y = y[!is.na(y)] y = y[1:tlen] return(y) } testsp <- function(x) { tau = 1 # threshold sp <- spectrum(x,log="no",ylim=c(0,10)) return(1/sp$freq[sp$spec>= tau]) } a = 120 # length of record # simple period y = getQuasiPeriod(40,0,1,0,0,0,a,1) CairoPNG("img/trend-period-noise-fig1.png") plot(y1,type="l",ylab="",xlab="", main="Sine Wave / No Noise") dev.off() # noisy amplitude y = getQuasiPeriod(40,0,1,0.5,0,0,a,1) y1 = getQuasiPeriod(40,0,1,0,0,0,a,1) CairoPNG("img/trend-period-noise-fig2.png") plot(y1,type="l",ylab="",xlab="",col="gray",ylim=c(min(y),max(y)), main="Sine Wave / Noisy Amplitude") lines(y,col="black") dev.off() # noisy period y = getQuasiPeriod(40,10,1,0,0,0,a,1) y1 = getQuasiPeriod(40,0,1,0,0,0,a,1) CairoPNG("img/trend-period-noise-fig3.png") plot(y1,type="l",ylab="",xlab="",col="gray",ylim=c(min(y),max(y)), main="Sine Wave / Noisy Period") lines(y,col="black") dev.off() # noisy series y = getQuasiPeriod(40,0,1,0,0,1,a,1) y1 = getQuasiPeriod(40,0,1,0,0,0,a,1) CairoPNG("img/trend-period-noise-fig4.png") plot(y1,type="l",ylab="",xlab="",col="gray",ylim=c(min(y),max(y)), main="Sine Wave / Noisy Signal") lines(y,col="black") dev.off() # noisy everything y = getQuasiPeriod(40,4,1,0.5,0,1,a,1) y1 = getQuasiPeriod(40,0,1,0,0,0,a,1) CairoPNG("img/trend-period-noise-fig5.png") plot(y1,type="l",ylab="",xlab="",col="gray",ylim=c(min(y),max(y)), main="Sine Wave / Noisy Period, Amplitude, and Signal") lines(y,col="black") dev.off() sp <- spectrum(y1,log="no") CairoPNG("img/trend-period-noise-fig6.png") plot(sp,log="no") dev.off() 1/sp$freq[sp$spec>2] sp <- spectrum(y,log="no") CairoPNG("img/trend-period-noise-fig7.png") plot(sp,log="no") dev.off() 1/sp$freq[sp$spec>2] # noisy everything y = getQuasiPeriod(40,4,1,0.5,0,1,a,5) y1 = getQuasiPeriod(40,0,1,0,0,0,a,1) CairoPNG("img/trend-period-noise-fig8.png") plot(y,type="l",ylab="",xlab="",col="black",ylim=c(min(y),max(y)), main="Sine Wave / Noisy Period, Amplitude, and Signal") dev.off() sp <- spectrum(y,log="no") CairoPNG("img/trend-period-noise-fig9.png") plot(sp,log="no") dev.off() 1/sp$freq[sp$spec>2] n = 1000 mc <- matrix(NA,ncol=n,nrow=a) #for (i in 1:n) { mc[,i] = getQuasiPeriod(11.25,1,1,0.05,1,1,a,3) } for (i in 1:n) { mc[,i] = getQuasiPeriod(7.5,0.5,1,0.01,1,1,a,10) } mct <- unlist(apply(mc,2,testsp)) h <- hist(mct,100) h$breaks[h$counts>0] h$counts[h$counts>0]