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

f=0:500/1000
lfp1detrend=lfp1-splinefit1
lfp30detrend=lfp30-splinefit30

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(smoothI30d, type = "l", lwd = 3, xlab = "Frequency", ylab = "", xaxt = "n", yaxt = "n", cex.lab = 1.8)
lines(smoothI1d, lwd = 3, 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.7)
axis(2, at = seq(from = 0, to = 8000, by = 2000), labels = seq(from = 0, to = 8000, by = 2000), cex.axis = 1.7)

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