library(EMD)



readCRU<-function(source="http://www.cru.uea.ac.uk/cru/data/temperature/hadcrut3vgl.txt" , temps=14,plot=T) {
	f<-read.table(source,fill=TRUE)
	d<-as.vector(t(as.matrix(f[seq(1,length(f[,1]),by=2),temps])))
	d<-d[d!= 0]
	na.omit(ts(d,frequency=length(temps),start=1850))
}

#Comment out this line after first run
CRU<-readCRU(temps=14)

myemd2<-function(Y,xlim=F,N=2,predict=F) {
	x<-time(Y)
	x1<-c(x,seq(end(Y),end(Y)+100),by=1/frequency(Y))
	if (!xlim) xlim=c(x[1],x[length(x)]+(x[length(x)]-x[1])/2)
	par(mfcol=c(1,1))
	E<-emd(Y,plot.imf=FALSE,boundary="periodic")
	H<-hilbertspec(E$imf)
	f<-apply(H$inst,MARGIN=2,FUN="mean")
	s<-apply(H$inst,MARGIN=2,FUN="sd")
	p<-1/f
	d<-diff(c(0,H$energy))
	r<-rbind(f,p,1/s,d)
	P<-NULL
	Sr<-ts(E$residue,start=start(Y))
	for (i in 1:E$nimf) {
		P<-c(P,m<-mean(H$instantfreq[,i]))
	}
# Plot of each imf
	par(mfrow=c(E$nimf+1, 1), mar=c(1,3,1,3)) 
	rangeimf <- range(E$imf)
	for(i in 1:E$nimf) {
		imf <- ts( E$imf[,i],start=start(Y))
		plot(imf, xlab="", ylab="", xaxt="n",main="",col="blue")
		text(start(Y),0.05,sprintf("IMF %d P=%.1f", i, 1/P[i]),pos=4)
		abline(h=0,lty=2)
	}
	par(mar=c(1,3,1,3)) 
	plot(ts(E$residue,start=start(Y)), xlab="", ylab="",xaxt="n", main="",col="blue")
	axis(1,pos=0)
	abline(h=0,lty=2)
	
	browser()
	P<-c(P,mean(hilbertspec(as.matrix(E$residue))$instantfreq))
	print(1/P)
	r<-range(Y);r[2]<-r[2]*3
	if (N==1) S<-E$imf[,E$nimf] else {
		ra<-(E$nimf-N+1):E$nimf
		S<-apply(E$imf[,ra],MARGIN=1,FUN="sum")
	}
		S<-ts(S,start=start(Y))

	
#plotting
	par(mfcol=c(1,1),mar=c(4,4,3,3))
	plot(Y,type="l",ylab="Anomaly Degrees C")
	lines(S,col=2,lty=2)	
	lines(Sr,col=2)
	lines(S+Sr,col=2,lwd=2)
	abline(h=0,lty=3,col=2)
	lines(c(1975,1975),c(-1,1),col="blue",lty=2)
	lines(c(2000,2000),c(-1,1),col="blue",lty=2)
	s1<-window(S,start=1950,end=2000)
	sr<-window(Sr,start=1950,end=2000)
	lines(s1,col="blue",lwd=2)
	lines(sr,col="blue",lwd=2)
	dr<-diff(range(sr))
	d1<-diff(range(s1))
	print(c(dr,d1,dr/(d1+dr)))
}

myemd2(CRU,N=2)