Skip to contents

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)
)
}