Robust, fast and fully automated peak fitting of Voigt profiles to spectral data

Implements the algorithm as discussed in <paperedoi>. The model fitted is defined as: $$f(x) = a x + b + \sum_{j=1}^K \beta_j V(x \mid \theta_j) + \sum_{i=1} \alpha_i f_i(x)$$ Here $\beta_j > 0$ are the amplitudes of the Voigt profiles $V$ where $\theta_j$ represents the position (pos), Gaussian (gwidth) and Lorentzian (lwidth) width of the Voigt profile. It requires no starting parameters, but optionally some can be passed in start_peaks.

The most impactful hyperparameters are typical_width and add_component_every_iters. The former should be a bit smaller than the typical width of Voigt profiles in the signal. If chosen too big or too small the fit will be usually suboptimal. The hyperparameter add_component_every_iters defines the number of iterations to perform until a new peak is added. Setting it to a bigger value comes at a greater computational cost but usually results in a better fit. The hyperparameter max_iter defines the total number of iterations, impact is similar but usually weaker than add_component_every_iters.

Additionally, it is possible to fit arbitrary positive background functions $f_i(x) > 0$, these are by default zero. Note that $\alpha_i > 0$ must hold. BSplines could be used here.

Implementation uses Voigt profiles as defined by 'RcppFaddeeva::Voigt()'.

Usage

spectralem(
  x,
  y,
  K,
  start_peaks = list(pos = c(), gwidth = c(), lwidth = c()),
  background_model = list(linear = TRUE),
  typical_width = diff(range(x))/100/2,
  add_component_every_iters = 10,
  max_iter = add_component_every_iters * (K + 2 + ceiling(K/2)),
  min_width = 1e-06 * mean(diff(x)),
  possible_peak_positions = range(x),
  print_progress = TRUE
)

Arguments

x: the x coordinates of the signal
y: the function values of the function to be fit evaluated at the x coordinates
K: the number of peaks to fit
start_peaks: optional initial peak parameters; list see examples
background_model: defines the background model
typical_width: typical width of a voigt profile in the signal
add_component_every_iters: number of iterations to perform until a new peak is added
max_iter: the maximal number of iterations
min_width: minimal gaussian and lorentzian width
possible_peak_positions: allowed range of positions for fitted peaks, can be increased to fit cut off peaks
print_progress: whether to print a progress bar

Value

A list with components:

fit_params - A list which contains named components amp, pos, gwidth and lwidth which are vectors of the same length. For each the i-th entry yields the parameters of the i-th Voigt profile. Further named components of the list a and b describe the slope and offset of the linear background respectively. If $f_i(x)$ are specified also contains a vector background_amps representing $\alpha_i$
fit - A vector containing the function values of the fitted function evaluated at x.

Examples

if (FALSE) {
res <- spectralem(
  x, y,
  K = 4,
)

res <- spectralem(
  x, y,
  K = 4,
  start_peaks = list(
    pos = c(332, 577, 697),
    gwidth = c(5, 5, 5),
    lwidth = c(0.1, 0.1, 0.1)
  )
)

res <- spectralem(
  x, y,
  K = 4,
  background_model = list(linear = TRUE, f1, f2, f3)
)
}