1 Introduction

The aim of spectral pre-treatment is to improve signal quality before modelling as well as remove physical information from the spectra. Applying a pre-treatment can increase the repeatability/reproducibility of the method, model robustness and accuracy, although there are no guarantees this will actually work. The pre-processing functions currently available in the package are listed in Table 1.

prospectr version 0.2.9 -- proxy

check the package repository at: https://l-ramirez-lopez.github.io/prospectr/,
https://github.com/l-ramirez-lopez/prospectr

We show below how they can be used, using the Near-Infrared (NIR) dataset (NIRsoil) included in the package (Fernandez-Pierna and Dardenne, 2008). Observations should be arranged row-wise.

library(prospectr)
data(NIRsoil)
# NIRsoil is a data.frame with 825 observations and 5 variables:
# Nt (Total Nitrogen), Ciso (Carbon), CEC (Cation Exchange Capacity),
# train (0/1 indicator for training (1) and validation (0) samples),
# spc (spectral matrix)
str(NIRsoil)

'data.frame':   825 obs. of  5 variables:
 $ Nt   : num  0.3 0.69 0.71 0.85 NA ...
 $ Ciso : num  0.22 NA NA NA 0.9 NA NA 0.6 NA 1.28 ...
 $ CEC  : num  NA NA NA NA NA NA NA NA NA NA ...
 $ train: num  1 1 1 1 1 1 1 1 1 1 ...
 $ spc  : num [1:825, 1:700] 0.339 0.308 0.328 0.364 0.237 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:825] "1" "2" "3" "4" ...
  .. ..$ : chr [1:700] "1100" "1102" "1104" "1106" ...

2 Noise removal

Noise represents random fluctuations around the signal that can originate from the instrument or environmental laboratory conditions. The simplest solution is to perform \(n\) repeated measurements and average the individual spectra; noise decreases by a factor of \(\sqrt{n}\). When this is not possible, or if residual noise remains, it can be removed mathematically.

2.1 Moving average

A moving average filter is a column-wise operation that averages contiguous wavelengths within a given window size, as illustrated in Figure 1.

# Add some noise
noisy <- NIRsoil$spc + rnorm(length(NIRsoil$spc), 0, 0.001)

# Plot the first spectrum
plot(
  x = as.numeric(colnames(NIRsoil$spc)),
  y = noisy[1, ],
  type = "l",
  lwd = 1.5,
  xlab = "Wavelength",
  ylab = "Absorbance"
)

# Window size of 11 bands; note the 5 first and last bands are lost
X <- movav(X = noisy, w = 11)
lines(x = as.numeric(colnames(X)), y = X[1, ], lwd = 1.5, col = "red")
grid()
legend(
  "topleft",
  legend = c("Raw", "Moving average"),
  lty = 1,
  col = c("black", "red")
)

2.2 Savitzky-Golay filtering

Savitzky-Golay filtering (Savitzky and Golay, 1964) is a widely used pre-processing technique. It fits a local polynomial regression on the signal and requires equidistant bandwidth. Mathematically, it operates as a weighted sum of neighbouring values:

\[x_j^* = \frac{1}{N}\sum_{h=-k}^{k} c_h x_{j+h}\]

where \(x_j^*\) is the smoothed value, \(N\) is a normalising coefficient, \(k\) is the number of neighbouring values on each side of \(j\), and \(c_h\) are pre-computed coefficients that depend on the chosen polynomial order, window size, and derivative degree (0 = smoothing, 1 = first derivative, 2 = second derivative).

# p = polynomial order
# w = window size (must be odd)
# m = m-th derivative (0 = smoothing)
# Accepts vectors, data frames, or matrices (observations row-wise)
sgvec <- savitzkyGolay(X = NIRsoil$spc[1, ], p = 3, w = 11, m = 0)
sg <- savitzkyGolay(X = NIRsoil$spc, p = 3, w = 11, m = 0)

# Bands at the edges of the spectral matrix are lost
dim(NIRsoil$spc)

[1] 825 700

dim(sg)

[1] 825 690

3 Derivatives

Taking numerical derivatives of the spectra can remove both additive and multiplicative effects and has further consequences summarised in Table 2.

First and second derivatives can be computed with the finite difference method (difference between two subsequent data points), provided that bandwidth is constant:

\[x_i' = x_i - x_{i-1}\]

\[x_i'' = x_{i-1} - 2 \cdot x_i + x_{i+1}\]

In R this can be achieved with the diff function from base, as shown in Figure 2:

d1 <- t(diff(t(NIRsoil$spc), differences = 1))  # first derivative
d2 <- t(diff(t(NIRsoil$spc), differences = 2))  # second derivative

plot(
  as.numeric(colnames(d1)), d1[1, ],
  type = "l", lwd = 1.5,
  xlab = "Wavelength", ylab = ""
)
lines(as.numeric(colnames(d2)), d2[1, ], lwd = 1.5, col = "red")
grid()
legend(
  "topleft",
  legend = c("1st derivative", "2nd derivative"),
  lty = 1, col = c("black", "red")
)

Derivatives tend to amplify noise. The gap derivative and the Savitzky-Golay algorithm both address this. The gap derivative is computed as:

\[x_i' = x_{i+k} - x_{i-k}\]

\[x_i'' = x_{i-k} - 2 \cdot x_i + x_{i+k}\]

where \(k\) is the gap size. This can be achieved using the lag argument of diff (see Figure 3):

# m = derivative order; w = gap size; s = segment size
gsd1 <- gapDer(X = NIRsoil$spc, m = 1, w = 11, s = 5)

plot(
  as.numeric(colnames(d1)), d1[1, ],
  type = "l", lwd = 1.5,
  xlab = "Wavelength (nm)", ylab = "1st derivative"
)
par(new = TRUE)
plot(
  as.numeric(colnames(gsd1)), gsd1[1, ],
  type = "l", lwd = 1.5, col = "red",
  xaxt = "n", yaxt = "n",
  xlab = "", ylab = ""
)
axis(4, col = "red")
mtext("Gap-segment derivative", side = 4, line = 3, col = "red")
grid()
legend(
  "topleft",
  legend = c("1st derivative", "Gap-segment 1st derivative"),
  lty = 1, col = c("black", "red")
)
par(new = FALSE)

For greater control over smoothing, the Savitzky-Golay (savitzkyGolay) and gap-segment derivative (gapDer) functions are preferable. The gap-segment algorithm first smooths the signal over a given segment size, then computes a derivative of a given order over a given gap size. See Luo et al. (2005) and Hopkins (2001) for details on these methods.

4 Scatter and baseline corrections

Undesired spectral variations due to light scatter effects and variations in effective path length can be removed using scatter corrections.

4.1 Standard Normal Variate (SNV)

The Standard Normal Variate (SNV) transformation (Barnes et al., 1989) is a row-wise standardisation that corrects for light scatter by centering and scaling each spectrum individually (Figure 4):

\[\text{SNV}_i = \frac{x_i - \bar{x}_i}{s_i}\]

where \(\bar{x}_i\) is the mean and \(s_i\) the standard deviation of the \(i\)-th spectrum (i.e. computed across all bands of that observation).

snv <- standardNormalVariate(X = NIRsoil$spc)
wav <- as.numeric(colnames(NIRsoil$spc))

par(mfrow = c(1, 2))

matplot(
  wav, t(NIRsoil$spc[1:5, ]),
  type = "l", lty = 1, lwd = 1.5,
  col = c("black", "grey20", "grey40", "grey60", "grey80"),
  xlab = "Wavelength (nm)", ylab = "Absorbance",
  main = "Raw"
)
grid()

matplot(
  wav, t(snv[1:5, ]),
  type = "l", lty = 1, lwd = 1.5,
  col = c("red", "tomato", "coral", "salmon", "lightsalmon"),
  xlab = "Wavelength (nm)", ylab = "SNV (Absorbance)",
  main = "SNV"
)
grid()

par(mfrow = c(1, 1))

According to Fearn (2008), SNV should be applied after filtering (e.g. Savitzky-Golay) rather than before.

4.2 Multiplicative Scatter Correction (MSC)

Along with SNV, MSC (Geladi et al., 1985) is one of the most widely used pre-processing techniques in NIR spectroscopy (Rinnan et al., 2009). It is a scatter correction method that accounts for additive and multiplicative effects by regressing each spectrum (\(x_i\)) against a reference spectrum (\(x_r\), typically the mean spectrum, Figure 5):

\[x_i = m_i x_r + a_i\]

\[\text{MSC}(x_i) = \frac{x_i - a_i}{m_i}\]

where \(a_i\) and \(m_i\) are the additive and multiplicative scatter coefficients, respectively, estimated by ordinary least squares for each spectrum.

# Reference spectrum is the column mean of the full matrix
msc_spc <- msc(X = NIRsoil$spc, ref_spectrum = colMeans(NIRsoil$spc))

plot(
  as.numeric(colnames(NIRsoil$spc)), NIRsoil$spc[1, ],
  type = "l", lwd = 1.5,
  xlab = "Wavelength (nm)", ylab = "Absorbance"
)
lines(as.numeric(colnames(msc_spc)), msc_spc[1, ], lwd = 1.5, col = "red")
grid()
legend(
  "topleft",
  legend = c("Raw", "MSC"),
  lty = 1, col = c("black", "red")
)

Because MSC correction depends on a reference spectrum, that reference must be stored and reused when processing new spectra. The following example shows how to apply the correction estimated on a training set to a separate validation set:

# Training and validation partitions
Xr <- NIRsoil$spc[NIRsoil$train == 1, ]
Xu <- NIRsoil$spc[NIRsoil$train == 0, ]

# Fit MSC on the training set (reference = mean of Xr by default)
Xr_msc <- msc(Xr)

# Apply the same reference spectrum to the validation set
Xu_msc <- msc(Xu, ref_spectrum = attr(Xr_msc, "Reference spectrum"))

4.3 SNV-Detrend

The SNV-Detrend transformation (Barnes et al., 1989) extends SNV by additionally correcting for wavelength-dependent scattering effects, i.e. residual curvature remaining after SNV. After the SNV transformation, a second-order polynomial is fitted to each spectrum and subtracted from it (Figure 6).

wav <- as.numeric(colnames(NIRsoil$spc))
dt  <- detrend(X = NIRsoil$spc, wav = wav)

# Dual y-axis: raw and detrended spectra are on different scales
plot(
  wav, NIRsoil$spc[1, ],
  type = "l", lwd = 1.5,
  xlab = "Wavelength (nm)", ylab = "Absorbance"
)
par(new = TRUE)
plot(
  wav, dt[1, ],
  type = "l", lwd = 1.5, col = "red",
  xaxt = "n", yaxt = "n",
  xlab = "", ylab = ""
)
axis(4, col = "red")
mtext("Detrended", side = 4, line = 3, col = "red")
grid()
legend(
  "topleft",
  legend = c("Raw", "SNV-Detrend"),
  lty = 1, col = c("black", "red")
)
par(new = FALSE)

4.4 Baseline removal

The baseline function estimates and removes the baseline of each spectrum in three steps, with the result illustrated in Figure 7:

1. The convex hull points of each spectrum are identified using grDevices::chull.
2. The convex hull points are linearly interpolated to the wavelength positions of the original spectrum to produce the baseline.
3. The baseline is subtracted from the original spectrum.

wav <- as.numeric(colnames(NIRsoil$spc))
bs <- baseline(NIRsoil$spc, wav)
fitted_baselines <- attr(bs, "baselines")

par(mfrow = c(1, 2))

# Raw spectra with fitted baselines
matplot(
  wav, t(NIRsoil$spc[1:3, ]),
  type = "l", lty = 1, lwd = 1.5,
  col = c("black", "grey40", "grey70"),
  xlab = "Wavelength (nm)", ylab = "Absorbance",
  main = "Raw and baseline"
)
matlines(wav, t(fitted_baselines[1:3, ]), lty = 2, col = c("black", "grey40", "grey70"))
grid()

# Baseline-corrected spectra
matplot(
  wav, t(bs[1:3, ]),
  type = "l", lty = 1, lwd = 1.5,
  col = c("red", "tomato", "salmon"),
  xlab = "Wavelength (nm)", ylab = "Absorbance",
  main = "Baseline corrected"
)
grid()

par(mfrow = c(1, 1))

5 Centering and scaling

Centering and scaling transform a given matrix so that columns have zero mean (centering), unit variance (scaling), or both (auto-scaling):

\[Xc_{ij} = X_{ij} - \bar{X}_{j}\]

\[Xs_{ij} = \frac{X_{ij} - \bar{X}_{j}}{s_{j}}\]

where \(Xc\) and \(Xs\) are the mean-centred and auto-scaled matrices, \(\bar{X}_j\) and \(s_j\) are the mean and standard deviation of variable \(j\). In R these operations are obtained with the scale function.

Spectroscopic models can often be improved by incorporating ancillary data (e.g. temperature) alongside the spectra (Fearn, 2010). However, due to the high dimensionality of spectral data, ancillary variables risk being dominated by the spectral block and contributing negligibly to the model (Eriksson et al., 2006). Block scaling addresses this by applying different scaling weights to different blocks of variables. With soft block scaling, each block is divided by a constant such that the sum of its column variances equals \(\sqrt{p}\), where \(p\) is the number of variables in the block. With hard block scaling, the block is scaled so that the sum of its column variances equals 1.

# type = "soft" or "hard"
# Returns a list with the scaled matrix (Xscaled) and the divisor (f)
bsc <- blockScale(X = NIRsoil$spc, type = "hard")$Xscaled
sum(apply(bsc, 2, var))  # should equal 1

[1] 1

A limitation of block scaling is that all variables within a block are brought to the same variance. When this is undesirable, sum-of-squares block weighting (blockNorm) provides an alternative: the spectral matrix is multiplied by a scalar to achieve a pre-specified sum of squares.

# targetnorm = desired Frobenius norm (sum of squares) for X
bn <- blockNorm(X = NIRsoil$spc, targetnorm = 1)$Xscaled
sum(bn^2)  # should equal 1

[1] 1

6 Resampling

To match the spectral response of one instrument to another, a signal can be resampled to new band positions by interpolation (resample) or by convolution using full-width at half-maximum (FWHM) values (resample2). Figure 8 shows an example of resampling a NIR spectrum to a coarser resolution of 10 nm using resample.

wav <- as.numeric(colnames(NIRsoil$spc))
wav_coarse <- seq(min(wav), max(wav), by = 10)

rs <- resample(X = NIRsoil$spc, wav = wav, new.wav = wav_coarse)

plot(
  wav, NIRsoil$spc[1, ],
  type = "l", lwd = 1.5,
  xlab = "Wavelength (nm)", ylab = "Absorbance"
)
lines(
  wav_coarse, rs[1, ], lwd = 1.5, col = "red", type = "b", pch = 19, cex = 0.5
)
grid()
legend(
  "topleft",
  legend = c("Original", "Resampled (10 nm)"),
  lty = 1, col = c("black", "red")
)

7 Other transformations

7.1 Continuum removal

The continuum removal technique was introduced by Clark and Roush (1984) to highlight absorption features by removing the effect of the overall spectral shape (albedo). The method identifies the convex-hull envelope (continuum) of each spectrum and normalises the original values against it (Figure 9):

\[\phi_{i} = \frac{x_{i}}{c_{i}}, \quad i = \{1, \ldots, p\}\]

where \(x_i\) and \(c_i\) are the original and continuum values at the \(i\)-th wavelength of a set of \(p\) wavelengths. The result \(\phi_i\) is dimensionless and bounded in \([0, 1]\) for reflectance spectra, where values close to 1 indicate no absorption relative to the continuum.

The continuumRemoval function handles both reflectance (type = "R", default) and absorbance (type = "A") spectra. For absorbance data, spectra are first converted to reflectance (\(1/X\)) before computing the continuum, and the result is back-transformed afterwards. This means that for absorbance data, continuumRemoval and baseline are not equivalent: they compute the convex-hull envelope on different scales. For reflectance data the two functions are more directly comparable, differing only in the final step: baseline subtracts the continuum (\(x_i - c_i\)), whereas continuumRemoval divides by it (\(x_i / c_i\)).

# type = "A" for absorbance, "R" for reflectance (default)
cr  <- continuumRemoval(X = NIRsoil$spc, type = "A")
wav <- as.numeric(colnames(NIRsoil$spc))

par(mfrow = c(1, 2))

matplot(
  wav, t(NIRsoil$spc[1:3, ]),
  type = "l", lty = 1, lwd = 1.5,
  col = c("black", "grey40", "grey70"),
  xlab = "Wavelength (nm)", ylab = "Absorbance",
  main = "Raw"
)
grid()

matplot(
  wav, t(cr[1:3, ]),
  type = "l", lty = 1, lwd = 1.5,
  col = c("red", "tomato", "salmon"),
  xlab = "Wavelength (nm)", ylab = "Continuum-removed",
  main = "Continuum removal"
)
grid()

par(mfrow = c(1, 1))

References