# R code for signal averaging

# load SignalNoise.RData and examine the four objects within:
#   noise: random noise from population with a mean of 0 and a 
#          standard deviation of 10
#   time: an index representing increments of time
#   signal: the pure signal, which consists of three Gaussian 
#           peaks at times of 256, 513, and 641
#   sig.tot: the sum of signal and time, or the noisy signal
load("SignalNoise.RData")
mean(noise)
sd(noise)

# plot the signal with noise and overlay the pure signal
plot(time, sig.tot, type = "l", col = "blue")
lines(time, signal, type = "l", col = "black", lwd = 4)

# estimate standard deviation of the noise using a signal-free
# area; the theoretical standard deviation of the noise is 10
plot(time[800:1000], sig.tot[800:1000], type = "l", col = "blue")
lines(time[800:1000], signal[800:1000], type = "l", col = "black", lwd = 4)
sd.n1 = sd(sig.tot[800:1000])
sd.n1

# estimate signal for the first peak; theoretical value is 100
plot(time[200:300], sig.tot[200:300], type = "l", col = "blue")
lines(time[200:300], signal[200:300], type = "l", col = "black", lwd = 4)
abline(v = 256, lty = 2, lwd = 2, col = "red")
m.n1 = mean(sig.tot[254:258])
m.n1

# estimate the S/N ratio; theoretical value is 100/10 = 10
sn1 = m.n1/sd.n1
sn1

# simulate signal averaging for n = 2 to 4 scans by adding random
# noise 2-4 times
avg2 = (2 * signal + rnorm(1024, 0, 10) +  rnorm(1024, 0, 10))/2
avg3 = (3 * signal + rnorm(1024, 0, 10) +  rnorm(1024, 0, 10)  
        + rnorm(1024, 0, 10))/3
avg4 = (4* signal + rnorm(1024, 0, 10) +  rnorm(1024, 0, 10)  
        + rnorm(1024, 0, 10) +  rnorm(1024, 0, 10))/4

# plot results for n = 1 to 4 scans
old.par = par(mfrow = c(2, 2))
plot(time, sig.tot, type = "l", col = "blue", ylim = c(-30,120), 
     main = "Number of Scans: 1")
plot(time, avg2, type = "l", col = "blue", ylim = c(-30,120), 
     main = "Number of Scans: 2")
plot(time, avg3, type = "l", col = "blue", ylim = c(-30,120), 
     main = "Number of Scans: 3")
plot(time, avg4, type = "l", col = "blue", ylim = c(-30,120), 
     main = "Number of Scans: 4")
par(old.par)

# estimate the S/N ratio for n = 1 to 4 scans; theoretical value
# is (n)^0.5 * 10
sd.n2 = sd(avg2[800:1000])
m.n2 = mean(avg2[254:258])

sd.n3 = sd(avg3[800:1000])
m.n3 = mean(avg3[254:258])

sd.n4 = sd(avg4[800:1000])
m.n4 = mean(avg4[254:258])

sn2 = m.n2/sd.n2
sn3 = m.n3/sd.n3
sn4 = m.n4/sd.n4

# create data frame to hold results for n = 1 to 4 scans
s.n_df = data.frame(c(1,2,3,4), 
                    c(m.n1, m.n2, m.n3, m.n4), 
                    c(sd.n1, sd.n2, sd.n3, sd.n4), 
                    c(sn1, sn2, sn3, sn4), 
                    10 * c(sqrt(1), sqrt(2), sqrt(3), sqrt(4)))

colnames(s.n_df) = c("number scans","mean signal", 
                     "std dev noise", "S/N found", "S/N Theory")

s.n_df