# 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