# V.CAN4_02: Biopsychological Theories of Personality =========================

# To issue commands, highlight them with the mouse cursor and hit 'Run' above.
# You can highlight more than one line to run several commands at once. To learn 
# more about the functions used, type help(function) or ?function! 

# 02.1: Eysenck/EEG -----------------------------------------------------------

# The following code generates the figure illustrating that a raw EEG signal is  
# composed of several frequencies.

# Required function for generating sine waves of a given frequency f with a 
# certain power (default = 1)) and duration (default = 1000); explanation:
# sin(seq(0,2*pi,length.out=...)) generates a sine wave of length 1/f data 
# points that is repeated 2*f times via rep(), of which the first 1000 data 
# points are returned.

freq <- function(f,duration=1000,power=1) {
  return(rep(sin(seq(0,2*pi,length.out = floor(duration/f))),2*f)[1:duration]*power)
}

# From which frequencies to sample? We want to have a 32 trials with different 
# compositions of delta, theta, alpha, beta, and gamma frequencies as well as of a 
# theta burst that will give rise to an raw.signals. Again, we use the function seq() that 
# via length.out = 32 generates 32 evenly spaced valueswithin the range given by the 
# first two arguments (here: the classical definitions of the frequency bands)   

delta = seq( 0.5, 4,length.out=32)
theta = seq( 4.1, 8,length.out=32)
alpha = seq( 8.1,13,length.out=32)
beta  = seq(13.1,30,length.out=32)
gamma = seq(30.1,99,length.out=32)

# We put them together in a data.frame

bands = data.frame(delta,theta,alpha,beta,gamma)

# Creating an raw.signals 
# The following lines generate via dnorm() a normal distribution vector of length 1000 
# with means and standard deviations as given by the last two arguments. The division by 
# this vectors maximum makes the maxmimum equal one and the multiplier in the end scales 
# the distribution.

# This generates a vector that defines the time poits of the raw.signalss for plotting. We start 
# with -200 ms, because usually, one defines some time window before a stimulus onset as
# baseline. We will need this for plotting only, I provide this here to make it 
# comprehensible that I set the means of the normal distributions 200 ms later than one 
# would expect them to occur. 

time  = seq(-200,799)

LPP   = dnorm(1:1000,800,80)/max(dnorm(1:1000,800,80))*2.5
N1    = dnorm(1:1000,300,25)/max(dnorm(1:1000,300,25))*5
N2    = dnorm(1:1000,400,25)/max(dnorm(1:1000,400,25))*5
P3    = dnorm(1:1000,500,50)/max(dnorm(1:1000,500,50))*20

# Let's see how this looks like via the plot() function! Execute the following commands 
# line by line:

# We plot the largest signal first and provide the plot(function with the x values 
# (first argument) and the y values (second argument) as well as with the plot typ ("l" = lines),
# y limits of the plot, the color of the line, and a label for the y axis.
plot(time, P3,type="l", ylim=c(-5,20),col=8, ylab="raw.signals")

# We want to have lines at time = 0 and voltage = 0, and as we need to do this later on 
# several times, we create a custom function:
zerolines <- function() {
  # par("usr) returns the x and y limits of the current plot, the first two are the x limits. 
  # We need them as limits for our lines. The rep() function repeats the 0 two times. We want
  # solid lines, thus line type is set to lty=1 (default, we could omit this statement, but 
  # you may want to have dashed  (lty=2) or dotted (lty=3) lines, just try ...)
  lines(par("usr")[1:2],rep(0,2),lty=1) # x limits first  -> horizontal line
  lines(rep(0,2),par("usr")[3:4],lty=1) # y limits second -> vertical lines
}
zerolines() # This function needs no arguments, as we always want it to do the same job.

# We now add the remaining lines, the first four of which in grey color (col=8).
lines(time,N1,col=8)
lines(time,N2,col=8)
lines(time,LPP,col=8)
lines(time,-N1-N2+P3+LPP, lwd=2) # the last one ist black (col=1 per default) and has a line width of lwd=2

# The combination of these raw.signalss will serve as a vector that "filters" a theta signal.
# as we want to have different theta "bursts" for every trial, we generate then 32 
# times via a for loop. We first need to setup an empty variable burst, and then - for
# 32 values between the lower theta frequencies 4.1 and 5 - we create a "filtered" theta 
# signal (via the power= argument to the function freq() defined above) and add it to 
# burst via rbind() that binds the generated vectors row by row. We want to see how the 
# results look line and thus setup an empty plot (type = "n") first:

plot(c(1,1000),c(-20,20),type="n",xlab="Time",ylab="Signal")

burst = NULL
for (i in seq(4.1,5,length.out=32)) { 
  burst = rbind(burst,freq(i,power=-N1-N2+P3+LPP)) 
  lines(burst[nrow(burst),],col=8) # this plots the last line of burst (i.e., current output of freq())
  }

# Now, we generate a set of signals corresponding to the frequencies defined above and 
# stored in  the data.frame bands. Again, we setup two empty variables, the first one of which 
# eventually will end up as a data.frame containing the generated matrices trial by trial and the 
# second one (raw.signals) ending up as one single matrix containing the averaged raw EEG signals for 
# every trial.
raw.data=raw.signals=NULL

# And now, for every line (trial) i in bands that has nrow=32 lines ...
for (i in 1:nrow(bands)) {
  # a new variable is generated with the first entry being the burst signal 
  # defined above for that trial ...
  raw=burst[i,]
  # and for every frequency j stored in the ncol=5 columns of bands ...
  for (j in 1:ncol(bands)) {
      # the respective sine wave for trial i and frequency band j is generated by freq()
      # and added to the raw via rbind().
      raw=rbind(raw,freq(bands[i,j]))
  }
  # Before going to the next trial, we create the raw EEG signal for the current
  # trial by averaging the sine waves over every column in raw via colMeans() 
  avg=colMeans(raw)
  # and add the average to raw.
  raw=rbind(raw,avg)
  # Also, we add the average to raw.signals to ease later plotting.
  raw.signals=rbind(raw.signals,avg)
  # We now name the rows of the raw matrix ...
  rownames(raw) = c("burst",colnames(bands),"avg")
  # And now a bit more complicated code: we want to have names for the entries in 
  # the raw.data data.frame corresponding to i. We achieve this by composing a command via
  #  eval(parse(text=paste(...,sep=""))). Type ?eval for more details ...
  eval(parse(text=paste("raw.data$id",i,"=raw",sep="")))
}

# Now for the actual plotting, we want to plot for three of the trials ...
idx=c(4,8,16)
# the sine waves of the five frequency bands, the theta "burst" and the averaged signal.
# This gives seven rows and three columns for plotting. These are set up by par(mfrow=c())
par(mfrow=c(7,3))
# Also, to have narrower margins, we use par(mar=c()), that has a default of about c(5,4,4,2)
# Otherwise, we would habe very tiny plots surrounded by large white margins.
par(mar=c(2,2,0.5,1))

# As plotting is done rowwise, we need to plot every for frequency band every trial first
# before going to the next band. This can be achieved using a for loop, and again, we need eval(parse(text= ...))
# BTW: xaxt="n" omits x axes for the first six lines and gives the x axis for the last line only.
for (i in idx) {  eval(parse(text=paste("plot(time,raw.data$id",i,"['delta',],type='l',ylim=c(-1.0,1.0),xaxt='n')",sep=""))) }
for (i in idx) {  eval(parse(text=paste("plot(time,raw.data$id",i,"['theta',],type='l',ylim=c(-1.0,1.0),xaxt='n')",sep=""))) }
for (i in idx) {  eval(parse(text=paste("plot(time,raw.data$id",i,"['alpha',],type='l',ylim=c(-1.0,1.0),xaxt='n')",sep=""))) }
for (i in idx) {  eval(parse(text=paste("plot(time,raw.data$id",i,"['beta', ],type='l',ylim=c(-1.0,1.0),xaxt='n')",sep=""))) }
for (i in idx) {  eval(parse(text=paste("plot(time,raw.data$id",i,"['gamma',],type='l',ylim=c(-1.0,1.0),xaxt='n')",sep=""))) }
for (i in idx) {  eval(parse(text=paste("plot(time,raw.data$id",i,"['burst',],type='l',ylim=c(-10.0,20),xaxt='n')",sep=""))) }
for (i in idx) {  eval(parse(text=paste("plot(time,raw.data$id",i,"['avg',  ],type='l',ylim=c(-2.0,4.0))",sep="")));zerolines() }
# We set mfrow to default and mar to have now very large lower and upper margins to make the 
# width to height ratio resemble the ones in the above plots.  
par(mfrow=c(1,1))
par(mar=c(10,2,10,1))

# We now compute the average over all trials, i.e. the ERP time course, and
# setup the plot frame.
raw.signals=colMeans(raw.signals)
plot(time,raw.signals,type="n",ylim=c(-3,4),xlab="Time")
# We now plot the individual trials in light grey ...
for(i in 1:nrow(raw.signals)) {
  lines(time,raw.signals[i,],col=grey(.8))
}
# add the zero lines
zerolines()
#and now add the ERP over trials as thicker black line 
lines(time,raw.signals,lwd=2)

# Last, we also set mar to default.
par(mar=c(5,4,4,2))

# Finally, we export the generated figures by clicking export in the plot window.
# I chose to export them as EPS and edited the figures a little bit in Adobe Illustrator,
# but as you might not own this tool (or some other vector graphics program), just save 
# the plots as PNG ... they will have no very good resolution, though.

# I you have questions about the code, just send me an email (alexander.strobel@tu-dresden.de)