# R code for Savitzky-Golay 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)

# create the SG filters
sg5 = c(-3, 12, 17, 12, -3)/35
sg9 = c(-21, 14, 39, 54, 59, 54, 39, 14, -21)/231
sg13 = c(-11, 0, 9, 16, 21, 24, 25, 24, 21, 16, 9, 0, -11)/143

# apply the SG filters to the noisy data
sim.sg5 = filter(sig.tot, sg5)
sim.sg9 = filter(sig.tot, sg9)
sim.sg13 = filter(sig.tot, sg13)

# 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.sg5, type = "l", col = "blue", 
     main = "Size of Filter: 5")
plot(time, sim.sg9, type = "l", col = "blue", 
     main = "Size of Filter: 9")
plot(time, sim.sg13, type = "l", col = "blue", 
     main = "Size of Filter: 13")
par(old.par)

# calculate the signal-to-noise ratios
sd.sg1 = sd(sig.tot[800:1000])
m.sg1 = mean(sig.tot[254:258])

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

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

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

sn1 = m.sg1/sd.sg1
sn5 = m.sg5/sd.sg5
sn9 = m.sg9/sd.sg9
sn13 = m.sg13/sd.sg13

sg_df = data.frame(c(1,5,9,13), 
                   c(m.sg1, m.sg5, m.sg9, m.sg13), 
                   c(sd.sg1, sd.sg5, sd.sg9, sd.sg13), 
                   c(sn1, sn5, sn9, sn13))

colnames(sg_df) = c("size SG filter","mean signal", 
                    "std dev noise", "S/N")

sg_df