# 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