# 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