Extended Multiplicative Signal Correction (EMSC)

Multiplicative Scatter Correction (MSC), introduced by Harald Martens, S. A. Jensen and Paul Geladi in 1983 [9] and popularised by Geladi, Douglas MacDougall and Martens in 1985 [1], separates chemical light absorption from physical light scattering in near-infrared reflectance spectra. The method was elegant and effective, but it rests on a strong assumption: each spectrum differs from a reference only by an offset and a stretch. Real spectra often contain more than that. Baselines tilt and bend, and other compounds that we are not interested in add their own bands.

In 1991, Martens and Edward Stark published an extension in the Journal of Pharmaceutical and Biomedical Analysis [2]. They kept the MSC regression and simply added more columns to it: smooth functions of wavelength to describe baseline curvature, and the spectra of known interferents so that their contribution could be estimated and subtracted. The result became known as Extended Multiplicative Signal Correction (EMSC), the name in the title of their paper. (You will also see it called extended multiplicative scatter correction; it is the same method.) In 2003, Martens, Nielsen and Engelsen showed how the extra terms could also describe wavelength-dependent scattering in powder mixtures [3].

The method found some of its most demanding applications in infrared microspectroscopy of biological samples. Achim Kohler, Martens and co-workers used EMSC to separate physical and chemical information in FTIR images of tissue sections [4], and then to model the strong distortions caused by Mie scattering when the size of cells is comparable to the wavelength of infrared light [5]. Building on that work, Bassan, Gardner, Kohler and colleagues developed the resonant Mie scattering correction (RMieS-EMSC) [6], now a standard tool in biomedical infrared spectroscopy. Kohler's group later refined it into Mie extinction EMSC (ME-EMSC) and released it as open-source code [10]. A tutorial by Nils Kristian Afseth and Kohler in 2012 made EMSC accessible to a much wider audience [7], and the same model-based approach has since been carried over to Raman spectra of biological samples [11].

MSC works well when scatter is the main problem and it behaves like a single offset and stretch. But real spectroscopic data often add other complications:

• Sloping or curved baselines, for example because scattering depends on wavelength, or because of instrumental drift.
• Known interfering compounds that vary from sample to sample: water vapour or CO₂ in an FTIR beam path, a packaging film, the paraffin in an embedded tissue section.
• Constituents with known spectra whose variation we want to keep, not to confuse with scatter.

EMSC handles all of these with the same idea MSC uses: describe the spectrum as a sum of known shapes, estimate how much of each shape is present with least squares, and remove the ones we do not want.

MSC models a spectrum ${\mathbf{x}}_i$ as:

EMSC adds a polynomial in wavelength and the spectra of known interferents:

Where:

• $a_i$ and $b_i$ are the same offset and multiplicative terms as in MSC,
• $\mathbf{u}$ is the wavelength axis rescaled to run from −1 to 1, and ${\mathbf{u}}^j$ means each element raised to the power $j$. So $d_{i1}\mathbf{u}$ is a straight tilt and $d_{i2}{\mathbf{u}}^2$ is a parabolic bend,
• ${\mathbf{m}}_k$ is the spectrum of the $k$-th known interferent, and $c_{ik}$ is how much more (or less) of it sample $i$ contains than the reference, after scaling by $b_i$. Because the reference is usually the mean spectrum, it already carries the average amount of every interferent, so $c_{ik}$ can be negative for samples with less than average.

Once all coefficients are estimated, we subtract everything we do not want and undo the multiplicative effect:

Notice what stays: the reference part $b_i{\mathbf{x}}_{\text{ref}}$ (which becomes ${\mathbf{x}}_{\text{ref}}$ after dividing by $b_i$) and the residual, which is where the chemical differences between this sample and the reference live.

The demo below makes this concrete. For a single spectrum, it shows each piece of the EMSC model: the scaled reference, the offset, the polynomial baseline, the interferent, and the residual (magnified five times so you can see it). Add them all up and you get the measured spectrum back.

Look at the residual for samples at the two ends of the colour scale (the darkest blue and the brightest amber). It has peaks exactly where the analyte bands are, positive for high-content samples and negative for low-content ones. That is the chemistry, and it is precisely what EMSC is designed to keep. Smaller bands at 1450 and 1940 nm, with the opposite sign, come from water, which also varies between the simulated samples. Notice also that the interferent coefficient $c$ can be negative: it measures the difference from the amount already present in the reference.

Use EMSC rather than plain MSC when you see:

• Baseline tilt or curvature that varies between samples and is still visible after MSC,
• wavelength-dependent scattering, typical of powders with different particle sizes and of biological cells in infrared microscopy,
• a known interferent whose amount varies: water vapour, CO₂, a solvent, a container, an embedding medium,
• samples whose chemical variation is large enough to disturb MSC, but whose main constituent spectra are known (see "good" spectra below).

Stick with MSC (or SNV) when the scatter is a simple offset and stretch, there is no visible residual baseline, and no known interferent varies. Simpler is better when it works.

EMSC is a multiple linear regression, done once per spectrum, with the same design matrix for every spectrum.

Step 1: choose the reference. Usually the mean spectrum of the training set, exactly as in MSC:

Step 2: build the design matrix. Put every shape the model is allowed to use in one column:

Each row of $\mathbf{Z}$ corresponds to one wavelength.

Step 3: fit every spectrum. Solve the least squares problem with the normal equations (or, better, a QR-based solver):

The vector ${\hat{\mathbf{\beta }}}_i$ contains $a_i,b_i,d_{i1},\dots ,d_{iJ},c_{i1},\dots ,c_{iK}$. Because $\mathbf{Z}$ is the same for every spectrum, the matrix $({\mathbf{Z}}^T\mathbf{Z}{)}^{-1}{\mathbf{Z}}^T$ can be computed once and reused.

Step 4: correct. Subtract the unwanted parts and divide by $b_i$, as in the formula above.

Step 5: store the model. The reference, the polynomial order, the interferent spectra and the minimum and maximum wavelengths used to build $\mathbf{u}$ are all part of the preprocessing model. New spectra are corrected with exactly the same ones. If the new spectra are later cut to a shorter wavelength range, reuse the stored minimum and maximum for $\mathbf{u}$; rescaling the shorter axis to run from −1 to 1 again would change the meaning of every polynomial coefficient.

Here the spectra suffer from three problems at once: an offset and stretch (which MSC can fix), a curved baseline that differs from sample to sample, and a varying amount of an interfering compound. The quality of each correction is again judged by how well the height of the analyte band follows the true content (R²).

Things to try:

1. Select MSC. The spectra are pulled together, but the curved baselines are left behind and the R² stays well below 1.
2. Select EMSC with order 0 and the interferent switched off. That is exactly MSC again.
3. Increase the order to 1. A straight tilt cannot describe a curved baseline, and here it even makes things worse (R² drops from about 0.81 to 0.71). Order 2 removes the bends and R² jumps to about 0.97.
4. Switch the interferent on. Its band at 1695 nm sits right next to the analyte band, so removing it explicitly pushes R² close to 1.
5. Set the baseline curvature to zero and the interferent to zero. Now MSC is already enough, and the extra terms bring nothing.

The polynomial order controls how flexible the baseline part is:

• Order 0: no polynomial. This is plain MSC.
• Order 1: adds a straight tilt.
• Order 2: adds a parabolic bend. A common default.
• Order 3 and above: rarely needed. A very flexible polynomial can start to imitate broad chemical bands and remove them.

A good rule is to start with order 2, look at the corrected spectra and at the residuals, and increase the order only if a clear baseline structure remains. Then check the choice with cross-validation of the model you actually care about.

If you know that a compound varies across samples and you do not want it in your model, include its pure spectrum as a column of $\mathbf{Z}$. EMSC will estimate how much of it each sample contains and subtract it. Martens and Stark called this spectral interference subtraction [2].

Include an interferent when:

• you have a reasonably pure spectrum of it, measured on the same wavelength axis,
• its amount varies between samples,
• it is not what you are trying to measure.

Do not include it when:

• it is constant across all samples (then it is already part of the reference),
• it is (or is strongly correlated with) the property you want to predict, because EMSC would remove the very signal you need.

The same trick works the other way around. If you know the spectra of the main chemical constituents, you can add them as columns too, but keep their contribution in the corrected spectrum [2],[7]. With "good" spectra ${\mathbf{g}}_l$ (coefficients $h_{il}$) in the model, the fit becomes

and the correction subtracts only the nuisance terms:

The good columns are fitted but not subtracted. Their job is to give chemical variation a column of its own, so that it is not mistaken for a change in $b_i$, which is the main weakness of MSC when samples differ a lot in composition. The good argument of the code below does exactly this.

There is a limit. If the reference, once the offset, polynomial and interferents are taken out, looks just like the analyte spectrum, EMSC cannot tell "more analyte" from "a larger $b_i$", and the correction removes the very signal you want. That happens, for example, when the only thing that varies in the chemistry is one constituent on a featureless background. Adding the analyte as a good spectrum helps, but it cannot fully cure a reference that is itself almost a combination of the other columns: that is the collinearity problem described in the limitations below.

import numpy as np
import matplotlib.pyplot as plt

def emsc(spectra, wavelengths, reference=None, poly_order=2,
         interferents=None, good=None):
    """
    Extended Multiplicative Signal Correction.

    spectra      : array (n_samples, n_wavelengths)
    wavelengths  : array (n_wavelengths,)
    reference    : array (n_wavelengths,), default: mean spectrum
    poly_order   : order of the polynomial baseline (0 = plain MSC)
    interferents : list of spectra to estimate and subtract
    good         : list of constituent spectra to estimate but keep
    Returns (corrected spectra, coefficients, reference).
    """
    spectra = np.asarray(spectra, dtype=float)
    if reference is None:
        reference = spectra.mean(axis=0)
    interferents = [] if interferents is None else list(interferents)
    good = [] if good is None else list(good)

    # Wavelength axis rescaled to [-1, 1] for numerical stability
    u = 2 * (wavelengths - wavelengths.min()) / np.ptp(wavelengths) - 1

    cols = [np.ones_like(u), reference]
    cols += [u ** j for j in range(1, poly_order + 1)]
    cols += [np.asarray(k, dtype=float) for k in interferents]
    n_remove = len(cols)            # columns 2 .. n_remove-1 are subtracted
    cols += [np.asarray(g, dtype=float) for g in good]   # fitted, but kept
    Z = np.column_stack(cols)

    coefs = np.linalg.lstsq(Z, spectra.T, rcond=None)[0].T    # one row per sample
    a, b = coefs[:, 0], coefs[:, 1]
    unwanted = coefs[:, 2:n_remove] @ Z[:, 2:n_remove].T       # polynomial + interferents
    corrected = (spectra - a[:, None] - unwanted) / b[:, None]
    return corrected, coefs, reference

# Example: scatter + curved baseline + a varying interferent
rng = np.random.default_rng(42)
wl = np.linspace(1100, 2500, 300)
u = 2 * (wl - wl.min()) / np.ptp(wl) - 1
analyte = np.exp(-((wl - 1730) ** 2) / (2 * 25 ** 2))
interferent = np.exp(-((wl - 1660) ** 2) / (2 * 18 ** 2))
# The sample matrix has bands of its own (it is not just a flat background)
matrix = (0.3 + 0.4 * np.exp(-((wl - 1450) ** 2) / (2 * 60 ** 2))
              + 0.5 * np.exp(-((wl - 1940) ** 2) / (2 * 70 ** 2)))

conc = rng.uniform(0, 1, 20)
raw = np.array([
    rng.uniform(-0.1, 0.1)                                    # offset
    + rng.uniform(0.8, 1.2) * (matrix + 0.2 * c * analyte)    # stretch x chemistry
    + rng.uniform(-0.05, 0.05) * u
    + rng.uniform(-0.05, 0.05) * u ** 2                       # curved baseline
    + rng.uniform(0, 0.2) * interferent                       # interferent
    + rng.normal(0, 0.002, wl.size)
    for c in conc
])

corrected, coefs, ref = emsc(raw, wl, poly_order=2, interferents=[interferent])
msc_only, _, _ = emsc(raw, wl, poly_order=0)                  # plain MSC, for comparison

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 4))
ax1.plot(wl, raw.T, lw=1); ax1.set_title('Raw')
ax2.plot(wl, corrected.T, lw=1); ax2.set_title('After EMSC')
plt.tight_layout()
plt.show()

band = np.argmin(abs(wl - 1730))
r2 = lambda X: np.corrcoef(X[:, band], conc)[0, 1] ** 2
print(f'R2 raw : {r2(raw):.3f}')        # about 0.27
print(f'R2 MSC : {r2(msc_only):.3f}')   # about 0.95
print(f'R2 EMSC: {r2(corrected):.3f}')  # about 0.999

% Example: scatter + curved baseline + a varying interferent
rng(42);
wl = linspace(1100, 2500, 300);
u = 2 * (wl - min(wl)) / (max(wl) - min(wl)) - 1;
analyte = exp(-((wl - 1730).^2) / (2 * 25^2));
interferent = exp(-((wl - 1660).^2) / (2 * 18^2));
% The sample matrix has bands of its own (it is not just a flat background)
matrix = 0.3 + 0.4 * exp(-((wl - 1450).^2) / (2 * 60^2)) ...
             + 0.5 * exp(-((wl - 1940).^2) / (2 * 70^2));

conc = rand(20, 1);
raw = zeros(20, numel(wl));
for i = 1:20
    raw(i, :) = (-0.1 + 0.2 * rand) ...
        + (0.8 + 0.4 * rand) * (matrix + 0.2 * conc(i) * analyte) ...
        + (-0.05 + 0.1 * rand) * u + (-0.05 + 0.1 * rand) * u.^2 ...
        + 0.2 * rand * interferent + 0.002 * randn(size(wl));
end

[corrected, coefs] = emsc(raw, wl, [], 2, {interferent});
mscOnly = emsc(raw, wl, [], 0, {});            % plain MSC, for comparison

subplot(1, 2, 1); plot(wl, raw'); title('Raw');
subplot(1, 2, 2); plot(wl, corrected'); title('After EMSC');

[~, band] = min(abs(wl - 1730));
R = corrcoef([raw(:, band), mscOnly(:, band), corrected(:, band), conc]);
fprintf('R2 raw : %.3f\nR2 MSC : %.3f\nR2 EMSC: %.3f\n', R(1, 4)^2, R(2, 4)^2, R(3, 4)^2);

function [corrected, coefs, reference] = emsc(spectra, wl, reference, polyOrder, interferents, good)
    % Extended Multiplicative Signal Correction. Rows are samples.
    % interferents: cell array of spectra to estimate and subtract
    % good:         cell array of constituent spectra to estimate but keep
    if nargin < 6, good = {}; end
    if isempty(reference), reference = mean(spectra, 1); end
    u = 2 * (wl(:) - min(wl)) / (max(wl) - min(wl)) - 1;

    Z = [ones(numel(u), 1), reference(:)];
    for j = 1:polyOrder
        Z = [Z, u.^j];
    end
    for k = 1:numel(interferents)
        Z = [Z, interferents{k}(:)];
    end
    nRemove = size(Z, 2);                         % columns 3..nRemove are subtracted
    for g = 1:numel(good)
        Z = [Z, good{g}(:)];                      % fitted, but kept
    end

    coefs = (Z \ spectra')';                      % one row per sample
    unwanted = coefs(:, 3:nRemove) * Z(:, 3:nRemove)';   % polynomial + interferents
    corrected = (spectra - coefs(:, 1) - unwanted) ./ coefs(:, 2);
end

emsc <- function(spectra, wl, reference = NULL, poly_order = 2,
                 interferents = list(), good = list()) {
  # Extended Multiplicative Signal Correction. Rows are samples.
  # interferents: spectra to estimate and subtract
  # good:         constituent spectra to estimate but keep
  if (is.null(reference)) reference <- colMeans(spectra)
  u <- 2 * (wl - min(wl)) / diff(range(wl)) - 1

  Z <- cbind(1, reference)
  if (poly_order > 0) for (j in 1:poly_order) Z <- cbind(Z, u^j)
  for (k in interferents) Z <- cbind(Z, k)
  n_remove <- ncol(Z)                          # columns 3..n_remove are subtracted
  for (g in good) Z <- cbind(Z, g)             # fitted, but kept

  coefs <- t(qr.solve(Z, t(spectra)))          # one row per sample
  unwanted <- matrix(0, nrow(spectra), ncol(spectra))
  if (n_remove > 2) {
    idx <- 3:n_remove
    unwanted <- coefs[, idx, drop = FALSE] %*% t(Z[, idx, drop = FALSE])
  }
  corrected <- (spectra - coefs[, 1] - unwanted) / coefs[, 2]
  list(corrected = corrected, coefs = coefs, reference = reference)
}

# Example: scatter + curved baseline + a varying interferent
set.seed(42)
wl <- seq(1100, 2500, length.out = 300)
u <- 2 * (wl - min(wl)) / diff(range(wl)) - 1
analyte <- exp(-((wl - 1730)^2) / (2 * 25^2))
interferent <- exp(-((wl - 1660)^2) / (2 * 18^2))
# The sample matrix has bands of its own (it is not just a flat background)
matrix_spec <- 0.3 + 0.4 * exp(-((wl - 1450)^2) / (2 * 60^2)) +
  0.5 * exp(-((wl - 1940)^2) / (2 * 70^2))

conc <- runif(20)
raw <- t(sapply(conc, function(c) {
  runif(1, -0.1, 0.1) + runif(1, 0.8, 1.2) * (matrix_spec + 0.2 * c * analyte) +
    runif(1, -0.05, 0.05) * u + runif(1, -0.05, 0.05) * u^2 +
    runif(1, 0, 0.2) * interferent + rnorm(length(wl), 0, 0.002)
}))

fit <- emsc(raw, wl, poly_order = 2, interferents = list(interferent))
msc_only <- emsc(raw, wl, poly_order = 0)$corrected   # plain MSC, for comparison

par(mfrow = c(1, 2))
matplot(wl, t(raw), type = "l", lty = 1, main = "Raw")
matplot(wl, t(fit$corrected), type = "l", lty = 1, main = "After EMSC")

band <- which.min(abs(wl - 1730))
r2 <- function(X) cor(X[, band], conc)^2
cat(sprintf("R2 raw : %.3f\nR2 MSC : %.3f\nR2 EMSC: %.3f\n",
            r2(raw), r2(msc_only), r2(fit$corrected)))

Strengths

• Handles more than offset and stretch. Tilts, curvature and wavelength-dependent scatter are described by the polynomial terms.
• Removes known interferents explicitly, instead of hoping a later model will learn to ignore them.
• Stays interpretable. Every coefficient has a name and a meaning, and the coefficients themselves (for example $b_i$ or the interferent amount $c_{ik}$) are often worth plotting.
• Can protect chemistry when the main constituent spectra are known.

Limitations

• More terms, more risk. A high polynomial order, or many interferents, gives the regression enough freedom to absorb real chemical signal.
• Needs a good reference, exactly like MSC.
• Needs good interferent spectra. A poorly measured interferent spectrum adds its own errors to every corrected spectrum.
• Collinearity. If an interferent spectrum looks very similar to the analyte or to a polynomial shape, their coefficients become unstable. The same happens when the reference itself, minus the offset, polynomial and interferents, is almost a copy of the analyte spectrum: then the analyte variation is absorbed into $b_i$ and removed. This is the same sensitivity problem discussed in the normal equations article.

• Try MSC first. Move to EMSC when MSC leaves visible baseline structure, or when you know an interferent varies [7],[8].
• Start with order 2. Increase only if a baseline structure clearly remains, and confirm with cross-validation.
• Check the coefficients. Plot $b_i$, $d_{ij}$ and $c_{ik}$ against sample number or known properties. Unexpected trends often reveal problems with sampling or measurement.
• Keep everything fixed for new data. Use the same reference, the same polynomial order and the same interferent spectra for the test set and for future samples. Never recompute the reference from test data.
• Do not stack EMSC with SNV. They correct similar effects; choose one.

EMSC is MSC with more columns. The measured spectrum is explained as an offset, a scaled reference, a smooth polynomial baseline and the spectra of known interferents, all estimated together by least squares. The unwanted pieces are subtracted, the multiplicative effect is divided out, and what remains is the reference plus the chemical differences.

It is powerful because it is explicit: you decide which shapes are "nuisance" and which are "chemistry". That same freedom is its main risk, so keep the model as simple as the data allow and validate the choice.