# load SignalNoise.RData and examine objects load("SignalNoise.RData") plot(time, sig.tot, type = "l", col = "blue") lines(time, signal, type = "l", col = "black", lwd = 2) # estimate noise signal-free area (theoretical 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 = 2) sd.n = sd(sig.tot[800:1000]) sd.n # estimate signal for first peak (theoretical 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 = 2) abline(v = 256, lty = 2, lwd = 2, col = "red") m.s = mean(sig.tot[254:258]) m.s # estimate S/N ratio (theoretical is 100/10 = 10) s.n = m.s/sd.n s.n # simulating signal averaging by adding together four individual trials sim = (4* signal + rnorm(1024, 0, 10) + rnorm(1024, 0, 10) + rnorm(1024, 0, 10) + rnorm(1024, 0, 10))/4 old.par = par(mfrow = c(2, 1)) plot(time, sig.tot, type = "l", col = "blue") plot(time, sim, type = "l", col = "blue") par(old.par) # estimate S/N ratio (theoretical is 2 * 100/10 = 20) sd.n = sd(sim[800:1000]) m.s = mean(sim[254:258]) s.n = m.s/sd.n m.s; sd.n; s.n # 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; filter is a vector # created using rep(1/w, w) where w is the filter's width ma5 = rep(1/5, 5) sim.ma5 = filter(sig.tot, ma5) old.par = par(mfrow = c(2, 1)) plot(time, sig.tot, type = "l", col = "blue") plot(time, sim.ma5, type = "l", col = "blue") par(old.par) sd.n = sd(sim.ma5[800:1000]) m.s = mean(sim.ma5[254:258]) s.n = m.s/sd.n m.s; sd.n; s.n # using R for a Savitzky-Golay filter sg5 = c(-3, 12, 17, 12, -3)/35 sim.sg5 = filter(sig.tot, sg5) old.par = par(mfrow = c(2, 1)) plot(time, sig.tot, type = "l", col = "blue") plot(time, sim.sg5, type = "l", col = "blue") par(old.par) sd.n = sd(sim.sg5[800:1000]) m.s = mean(sim.sg5[254:258]) s.n = m.s/sd.n m.s; sd.n; s.n # using R for Fourier filtering help(fft) test = c(1, 2, 3, 4) test # complete the FT to move from one domain to the other ft.test = fft(test) ft.test # complete the inverse FT to return to the original domain ift.test = fft(ft.test, inverse= TRUE) ift.test ift.test = fft(ft.test, inverse= TRUE)/length(test) ift.test Re(ift.test) # apply fft to noisy signal sim.ft = fft(sig.tot) head(sim.ft) plot(time, Re(sim.ft), type = "l", col = "blue") plot(time[1:512], Re(sim.ft)[1:512], type = "l", col = "blue") sim.ft[101:924] = 0 + 0i plot(time[1:512], Re(sim.ft)[1:512], type = "l", col = "blue") plot(time, Re(sim.ft), type = "l", col = "blue") sim.ift = fft(sim.ft, inverse = TRUE)/length(sim.ft) old.par = par(mfrow = c(2, 1)) plot(time, sig.tot, type = "l", col = "blue") plot(time, Re(sim.ift), type = "l", col = "blue") par(old.par) sd.n = sd(Re(sim.ift[800:1000])) m.s = mean(Re(sim.ift[254:258])) s.n = m.s/sd.n m.s; sd.n; s.n # using R to remove background signals rm(list = ls()) load("BkgdRemoval.RData") plot(wavelength, y1, ylim = c(0, 120), type = "l", col = "blue") lines(wavelength, y2, type = "l", col = "blue") lines(wavelength, y3, type = "l", col = "blue") lines(wavelength, y4, type = "l", col = "blue") lines(wavelength, y5, type = "l", col = "blue") abline(v = wavelength[207], lty = 2, col = "red") y.max = c(y1[207], y2[207], y3[207], y4[207], y5[207]) plot(conc, y.max, pch = 19, col = "blue", ylim = c(0, 120), xlim = c(0, 1)) lm.raw = lm(y.max ~ conc) summary(lm.raw) abline(lm.raw, col = "blue") # apply a Savitzky-Golay first-derivative filter fd5 = c(-2, -1, 0, 1, 2)/10 y1.fd5 = filter(y1, fd5) y2.fd5 = filter(y2, fd5) y3.fd5 = filter(y3, fd5) y4.fd5 = filter(y4, fd5) y5.fd5 = filter(y5, fd5) plot(wavelength, y1.fd5, ylim = c(-10, 10), type = "l", col = "blue") lines(wavelength, y2.fd5, type = "l", col = "blue") lines(wavelength, y3.fd5, type = "l", col = "blue") lines(wavelength, y4.fd5, type = "l", col = "blue") lines(wavelength, y5.fd5, type = "l", col = "blue") abline(v = wavelength[213], lty = 2, col = "red") abline(v = wavelength[201], lty = 2, col = "red") y.pp = c(y1.fd5[213] - y1.fd5[201], y2.fd5[213] - y2.fd5[201], y3.fd5[213] - y3.fd5[201], y4.fd5[213] - y4.fd5[201], y5.fd5[213] - y5.fd5[201]) plot(conc, y.pp, pch = 19, col = "blue", ylim = c(0, 15), xlim = c(0, 1)) lm.bkgd = lm(y.pp ~ conc) summary(lm.bkgd) abline(lm.bkgd, col = "blue")