lfp=read.table("data/lfp-ryan.dat")[,1]
library("splines")
time=1:1000
lfp1=lfp[1:1000]
lfp30=lfp[29001:30000]
splinefit1=lm(lfp1 ~ ns(time,knots=c(200,400,600,800)))$fit
splinefit30=lm(lfp30 ~ ns(time,knots=c(200,400,600,800)))$fit

I30=spec.pgram(lfp30,plot=FALSE)$spec
f=0:500/1000
lfp1detrend=lfp1-splinefit1
lfp30detrend=lfp30-splinefit30
plot(time,lfp1detrend,type="l",lwd=2)
plot(time,lfp30detrend,type="l",lwd=2)
I1detrend=spec.pgram(lfp1detrend,plot=FALSE,detrend=FALSE)$spec
I30detrend=spec.pgram(lfp30detrend,plot=FALSE,detrend=FALSE)$spec
smoothI1d=ksmooth(f[1:50],I1detrend[1:50],"normal",bandwidth=.002)
smoothI30d=ksmooth(f[1:50],I30detrend[1:50],"normal",bandwidth=.002)

plot(f[1:50],I1detrend[1:50],type="l",xlab="frequency",ylab=" ",lwd=1)
lines(smoothI1d,col="blue",lwd=2)

par(oma = rep(2, 4))
par(mar = c(5, 5, 1, 1))
par(mfrow = c(2, 1))

plot(f[1:50], I1detrend[1:50], type = "l", xlab = "Frequency", ylab = "",
     xaxt = "n", yaxt = "n", main = "", cex.lab = 1.6, lwd = 2)
lines(smoothI1d, lwd = 2, col = 4)

axis(1, at = seq(from = 0, to = 0.05, by = 0.01),
     labels = seq(from = 0, to = 0.05, by = 0.01), cex.axis = 1.5)
axis(2, at = seq(from = 0, to = 6000, by = 3000),
     labels = seq(from = 0, to = 6000, by = 3000), cex.axis = 1.5)

plot(f[1:50], I30detrend[1:50], type = "l", xlab = "Frequency", ylab = "",
     xaxt = "n", yaxt=  "n", main = "", cex.lab = 1.6, lwd = 2)
lines(smoothI30d, lwd = 2, col = 4)
axis(1, at = seq(from = 0, to = 0.05, by = 0.01),
     labels = seq(from = 0, to = 0.05, by = 0.01), cex.axis = 1.5)
axis(2, at = seq(from = 0, to = 12000, by = 6000),
     labels = seq(from = 0, to = 12000, by = 6000), cex.axis = 1.5)

dev.print(device = postscript, "18.6.eps", horizontal = TRUE)