# R code for moving average smoothing

# 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)

# using R for a moving average filter
help(filter)

# default is convolution and sides = 2, which takes an odd number
# of points and replaces the middle point with the average value;
# the filter is a vector created using the command rep(1/w, w)
# where w is the filter's width; here we create filters with
# sizes of 5, 9, and 13
ma5 = rep(1/5, 5)
ma9 = rep(1/9, 9)
ma13 = rep(1/13, 13)

# now we apply the filter to the noisy signal
sim.ma5 = filter(sig.tot, ma5)
sim.ma9 = filter(sig.tot, ma9)
sim.ma13 = filter(sig.tot, ma13)

# plot the results
old.par = par(mfrow = c(2, 2))
plot(time, sig.tot, type = "l", col = "blue", 
     main = "Size of Filter: 0")
plot(time, sim.ma5, type = "l", col = "blue", 
     main = "Size of Filter: 5")
plot(time, sim.ma9, type = "l", col = "blue", 
     main = "Size of Filter: 9")
plot(time, sim.ma13, type = "l", col = "blue", 
     main = "Size of Filter: 13")
par(old.par)

# determine S/N rations
sd.ma1 = sd(sig.tot[800:1000])
m.ma1 = mean(sig.tot[254:258])

sd.ma5 = sd(sim.ma5[800:1000])
m.ma5 = mean(sim.ma5[254:258])

sd.ma9 = sd(sim.ma9[800:1000])
m.ma9 = mean(sim.ma9[254:258])

sd.ma13 = sd(sim.ma13[800:1000])
m.ma13 = mean(sim.ma13[254:258])

sn1 = m.ma1/sd.ma1
sn5 = m.ma5/sd.ma5
sn9 = m.ma9/sd.ma9
sn13 = m.ma13/sd.ma13

ma_df = data.frame(c(1,5,9,13), 
                   c(m.ma1, m.ma5, m.ma9, m.ma13), 
                   c(sd.ma1, sd.ma5, sd.ma9, sd.ma13), 
                   c(sn1, sn5, sn9, sn13))

colnames(ma_df) = c("number scans","mean signal", 
                    "std dev noise", "S/N")

ma_df