# R code for background removal

# load BkgdRemoval.RData and examine the seven objects within:
#   conc: concentrations for five standard solutions of analyte
#   wavelength: wavelengths used for spectra
#   y1-y5: the spectra for the five standards
load("BkgdRemoval.RData")

# first, let's look at the data
plot(wavelength, y1, ylim = c(0, 120), type = "l", col = "blue")
lines(wavelength, y2, lty = 2, col = "blue")
lines(wavelength, y3, lty = 3, col = "blue")
lines(wavelength, y4, lty = 4, col = "blue")
lines(wavelength, y5, lty = 5, col = "blue")
legend(x = "topleft", legend = c("std 1", "std 2", "std 3", "std 4", "std 5"), col = "blue", lty = c(1, 2, 3, 4, 5), bty = "n")

# calibration curve using data as provided
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)

# replot spectra after removing background
plot(wavelength, y1.fd5, ylim = c(-7, 7), type = "l", col = "blue")
lines(wavelength, y2.fd5, lty = 2, col = "blue")
lines(wavelength, y3.fd5, lty = 3, col = "blue")
lines(wavelength, y4.fd5, lty = 4, col = "blue")
lines(wavelength, y5.fd5, lty = 5, col = "blue")
legend(x = "topleft", legend = c("std 1", "std 2", "std 3", "std 4", "std 5"), col = "blue", lty = c(1, 2, 3, 4, 5), bty = "n")

# replot the calibration curve
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")