feemscatter {albatross}

Handle scattering signal in FEEMs

Description

Remove or interpolate scattering signal in individual FEEM objects, FEEM cube objects, or lists of them.

Usage

  feemscatter(x, ...)
  ## S3 method for class 'list'
feemscatter(x, ..., cl, progress = TRUE)
  ## S3 method for class 'feemcube'
feemscatter(x, ..., cl, progress = TRUE)
  ## S3 method for class 'feem'
feemscatter(
    x, widths, method = c("omit", "pchip", "loess", "kriging", "whittaker"),
    add.zeroes = 30, Raman.shift = 3400, ...
  )

Arguments

x

An individual FEEM object, FEEM cube object, or a list of them, to handle the scattering signal in.

widths

A numeric vector or a list containing the half-widths of the scattering bands, in nm. Rayleigh scattering is followed by Raman scattering, followed by second diffraction order for Rayleigh and Raman, and so on. (Typically, there's no need for anything higher than third order, and even that is rare.) For example:

1. Rayleigh scattering

2. Raman scattering

3. Rayleigh scattering, 2\lambda

4. Raman scattering, 2\lambda

5. ...

For higher diffraction orders, the peak widths are proportionally scaled, making it possible to provide the same number for all kinds of scattering visible in the EEM. Set a width to 0 if you don't want to handle this particular kind of scattering signal.

It's possible to specify the bands asymmetrically. If the area to be corrected should range from x nm to the left of the scattering peak to y nm to the right of it, pass a list instead of a vector, and put a two-element vector c(x, y) for the appropriate kind of scattering. For example, passing widths = list(c(30, 20), 20) means “handle -30 nm to the left and +20 nm to the right of Rayleigh peak and \pm 20 nm around Raman peak”.

To sum up, given two half-widths W_1 and W_2, the test for being inside a kth diffraction order of a scattering band is as follows:

-W_1 < \frac{\lambda_\mathrm{center}}{k} - \lambda_\mathrm{em} < +W_2

In this example, a much larger portion of the anti-Stokes area is removed near the first order Rayleigh scattering signal than in the Stokes area. This can be useful to get rid of undesired signal where no fluorescence is observed on some spectrometers. The second and third order scattering signal areas are automatically scaled 2 and 3 times, respectively.

method

A string choosing how to handle the scattering signal:

omit

Replace it with NA.

pchip

Interpolate it line-by-line using piecewise cubic Hermitean polynomials (pchip). Pass a by argument to choose the direction of interpolation; see Details.

loess

Interpolate it by fitting a locally weighted polynomial surface (loess). Extra arguments are passed verbatim to loess, which may be used to set parameters such as span.

kriging

Interpolate it by means of ordinary or simple Kriging, as implemented in pracma function kriging. Pass a type argument to choose between the two methods. This method is not recommended due to its high CPU time and memory demands: it has to invert a dense O(N^2) matrix (which easily reaches multiple gigabytes for some EEMs), and compute its product with a vector then take scalar products O(N) times, with N = length(x).

whittaker

Interpolate it by minimising a weighted sum of squared residuals (for known part of the spectrum) and roughness penalty (squared central difference approximations for derivatives by \lambda_\mathrm{em} and \lambda_\mathrm{ex}). See Details for more information and parameters.

add.zeroes

Set intensities at \lambda_\mathrm{em} < \lambda_\mathrm{ex} - \mathtt{add.zeroes}\:\mathrm{nm} to 0 unless they have been measured in order to stabilise the resulting decomposition (Thygesen, Rinnan, Barsberg, and Møller 2004). Set to NA to disable this behaviour.

Raman.shift

Raman shift of the scattering signal of water, \textrm{cm}^{-1}.

...

Passed verbatim from feemscatter generics to feemscatter.feem.

If “pchip” method is selected, the by parameter chooses between interpolating by row, by column, or averaging both, see Details.

If “loess” method is selected, remaining options are passed to loess (the span parameter is of particular interest there).

If “kriging” method is selected, remaining options are passed to kriging.

If “whittaker” method is selected, available parameters include d, lambda, nonneg and logscale, see Details.

cl

If not missing, a parallel cluster object to run the scattering correction code on or NULL for the default cluster object registered via setDefaultCluster.

progress

Set to FALSE to disable the progress bar.

Details

The “pchip” method works by default as described in (Bahram, Bro, Stedmon, and Afkhami 2006): each emission spectrum at different excitation wavelengths is considered one by one. Zeroes are inserted in the corners of the spectrum if they are undefined (NA) to prevent extrapolation from blowing up, then the margins are interpolated using the corner points, then the rest of the spectrum is interpolated line by line. Since pchip requires at least 3 points to interpolate, the function falls back to linear interpolation if it has only two defined points to work with. The by argument controls whether the function proceeds by rows of the matrix (“emission”, default), by columns of the matrix (“excitation”), or does both (“both”) and averages the results to make the resulting artefacts less severe (Pucher, Wünsch, Weigelhofer, Murphy, Hein, and Graeber 2019) (see the staRdom package itself).

The “loess” method feeds the whole FEEM except the area to be interpolated to loess, then asks it to predict the remaining part of the spectrum. Any negative values predicted by loess are replaced by 0.

The “kriging” method (Press, Teukolsky, Vetterling, and Flannery 2007) is much more computationally expensive than the previous two, but, on some spectra, provides best results, not affected by artefacts resulting from line-by-line one-dimensional interpolation (pchip) or varying degrees of smoothness in different areas of the spectrum (loess). Any negative values returned by kriging are replaced by 0.

Whittaker smoothing

The “whittaker” method (Krylov and Labutin 2023) works by minimising a sum of penalties, requiring the interpolated surface to be close to the original points around it and to be smooth in terms of derivatives by \lambda_\mathrm{em} and \lambda_\mathrm{ex}.

The parameters d and lambda should be numeric vectors of the same length, corresponding to the derivative orders (whole numbers \ge 1) and their respective penalty weights (small real numbers; larger is smoother). For interpolation purposes, the default penalty is 10^{-2} \mathbf{D}_1 + 10 \mathbf{D}_2 , which corresponds to d = 1:2 and lambda = c(1e-2, 10).

Any resulting negative values are pulled towards 0 by adding zero-valued points with weight nonneg (default 1) and retrying. Set nonneg to 0 to disable this behaviour. It is also possible to deal with resulting negative values by scaling and shifting the signal between logscale (typically) and 1, interpolating the logarithm of the signal, then undoing the transformation. By default logscale is NA, disabling this behaviour, since it may negatively affect the shape of interpolated signal.

See the internal help page whittaker2 for implementation details.

Value

An object of the same kind (FEEM object / FEEM cube / list of them) with scattering signal handled as requested.

References

Bahram M, Bro R, Stedmon C, Afkhami A (2006). “Handling of Rayleigh and Raman scatter for PARAFAC modeling of fluorescence data using interpolation.” Journal of Chemometrics, 20(3-4), 99-105. doi:10.1002/cem.978.

Krylov IN, Labutin TA (2023). “Recovering fluorescence spectra hidden by scattering signal: in search of the best smoother.” Spectrochimica Acta Part A: Molecular and Biomolecular Spectroscopy, 122441. doi:10.1016/j.saa.2023.122441.

Press WH, Teukolsky SA, Vetterling WT, Flannery BP (2007). “Interpolation by Kriging.” In Numerical recipes: The Art of Scientific Computing (3rd Ed.), chapter 3.7.4, 144-147. Cambridge University Press, New York.

Pucher M, Wünsch U, Weigelhofer G, Murphy K, Hein T, Graeber D (2019). “staRdom: Versatile Software for Analyzing Spectroscopic Data of Dissolved Organic Matter in R.” Water, 11(11), 2366. doi:10.3390/w11112366.

Thygesen LG, Rinnan Å, Barsberg S, Møller JKS (2004). “Stabilizing the PARAFAC decomposition of fluorescence spectra by insertion of zeros outside the data area.” Chemometrics and Intelligent Laboratory Systems, 71(2), 97-106. ISSN 0169-7439, doi:10.1016/j.chemolab.2003.12.012.

See Also

feem, feemcube

Examples

  data(feems)
  plot(x <- feemscatter(
    feems[[1]], widths = c(25, 25, 20, 20),
    method = 'whittaker', Raman.shift = 3500
  ))