PLS algorithm

The concept of PLS originated from the Whittaker smoother [10], which was initially designed to smooth one-dimensional equally spaced data. Consider a data sequence y of length N, sampled at equal time intervals. The algorithm smoothes a noisy sequence y into a smooth sequence z by minimizing the following objective function: $S(\text{z})=(\text{y}-\text{z}{)}^{\text{T}}(\text{y}-\text{z})+\lambda {\text{z}}^{\text{T}}{\text{D}}^{\text{T}}\text{DZ}\mathrm{,}$ (1) where D is an m-th-order difference matrix with N - m rows and N columns. In this study, the discussion is limited to the 2nd-order difference (m = 2). Thus, D is expressed as $\text{D}=\left(\begin{array}{cccccccc}1& -2& 1& 0& \cdots & 0& 0& 0\\ 0& 1& -2& 1& \cdots & 0& 0& 0\\ ⋮& ⋮& ⋮& ⋮& \ddots & ⋮& ⋮& ⋮\\ 0& 0& 0& 0& \cdots & 1& -2& 1\end{array}\right)\mathrm{.}$ (2) The first term in Eq. (1) represents the fidelity of z to data y, whereas the second term represents the smoothness of z. The balance between fidelity and smoothness is controlled by the parameter λ. A larger λ leads to a smoother z.

In a typical spectrum, some data points have only baseline and noise components, whereas the rest of the points contain additional signal components. For baseline estimation, it is optimal to obtain a slowly varying baseline if there is a way to distinguish between both cases by considering only the former and ignoring the latter. Therefore, combining the Whittaker smoother with the concept of asymmetric least squares [17], Eilers and Boelens [14] first introduced a weight for each point to obtain a smooth baseline. Let W be a diagonal matrix with a weight array w of data points on its diagonal. Thus, the objective function in Eq. (1) becomes $S(\text{z})=(\text{y}-\text{z}{)}^{\text{T}}\text{W}(\text{y}-\text{z})+\lambda {\text{z}}^{\text{T}}{\text{D}}^{\text{T}}\text{DZ}\mathrm{.}$ (3) Setting the partial derivatives $\frac{\partial S}{\partial {\text{z}}^{\text{T}}}=-2\text{W}(\text{y}-\text{z})+2\lambda {\text{D}}^{\text{T}}\text{DZ}$ (4) to zero, we obtain the solution for minimizing Eq. (3): $\text{z}=(\text{W}+\lambda {\text{D}}^{\text{T}}\text{D}{)}^{-1}\text{W}\text{y}\mathrm{.}$ (5) Intuitively, the weight should be 0 if the data point is definitely in the peak region and 1 if the data point does not contain any peak component of the signal. However, the regions containing no signal are difficult to determine because of the presence of noise and peaks. Thus, several methods have been proposed using iterative algorithms to circumvent peak finding.

A pioneering method is the asymmetric least squares (AsLS) method developed by Eilers and Boelens [14]. In this method, an asymmetry parameter p, recommended to be set between 0.001 and 0.01, is introduced, and the weights are assigned as $w_i=\{\begin{array}{ll}p\hfill & y_i>z_i\mathrm{,}\hfill \\ 1-p\hfill & y_i⩽z_i\mathrm{.}\hfill \end{array}$ (6) The algorithm first sets w to an array of all 1s and solves Eq. (5) to obtain the initial z. Then, the data points are divided into two zones: positive deviations yi > zi and negative residuals $y_i⩽z_i$. By applying the asymmetric weight function given by Eq. (6), the positive deviations, which are supposed to contain signals, gain fewer weights than the negative residuals in the next iteration. The baseline, i.e., z, is updated by solving Eqs. (5) using the new w. By applying this new z, we can again obtain the updated weight array w. The iteration ends when the weight array w remains unchanged or the number of iterations reaches the maximum pre-set value. The final updated baseline z is the result of the estimation method.

Given that the AsLS method only assigns two values to the weights of each point, this assignment is unreasonable for data containing different signal peak heights. To improve the weight assignment, the adaptive iteratively reweighted penalized least squares (airPLS) [15] and asymmetrically reweighted penalized least squares (arPLS) [16] methods were successively developed.

The airPLS algorithm considers the effects of signal peak heights. It uses an exponential function to ensure that the weights are small for large peak heights and vice versa. This method also adapts to two-segment weight assignment. The subtraction d = y - z, d_- is part of d in the region where $\left\{d_i\right|d_i=y_i-z_i<\mathrm{0,}{d}_i\in d\}$, and d₊ is part of the region where $\left\{d_i\right|d_i=y_i-z_i>\mathrm{0,}{d}_i\in d\}$. The weight array w for each iteration step t in the airPLS method is constructed as follows: $w_i=\{\begin{array}{ll}0\hfill & y_i⩾z_i\hfill \\ \text{exp}\left[\frac{t(y_i-z_i)}{{\displaystyle {\sum }_{d_j\in d_{-}}\left|d_j\right|}}\right]\hfill & y_i<z_i\mathrm{.}\hfill \end{array}$ (7) The step number t in Eq. (7) is introduced to speed up the convergence of the iterations. The iteration of this method stops when the maximum number of iterations or the following termination criterion is fulfilled: ${\displaystyle \sum _{d_j\in d_{-}}}\left|d_j\right|<0.001\times \left|y\right|\mathrm{.}$ (8) Given that both the airPLS and AsLS methods attempt to set the weights of the data points to zero, or very close to zero, they tend to result in a baseline lower than the accurate baseline. To overcome this problem, the arPLS method is proposed, in which an asymmetric weight array is constructed as follows: $w_i=\{\begin{array}{cc}\frac{1}{1+\text{exp}\left[\frac{2}{{\sigma }_{-}}((y_i-z_i)-(-m_{-}+2{\sigma }_{-}))\right]}& y_i⩾z_i\\ 1& y_i⩽z_i\mathrm{,}\end{array}\mathrm{,}$ (9) where m_- and σ_- are respectively the mean and standard deviation of d_-.

The arPLS method performs baseline correction as successfully as the AsLS and airPLS methods without a priori peak detection. Indeed, it obtains a more satisfactory baseline in most cases. Similar to the previous method, the weight assignment is still constructed in different forms for both zones (d_- and d₊). In addition, it does not consider a specific form of the weights at the points corresponding to different signal peak heights. When dealing with higher and wider peaks in a rapidly changing baseline, the results are no longer reliable. In the next section, we propose a unified method to construct weights for the PLS method by invoking the Bayesian rule. The new algorithm outperforms arPLS when it comes to handling rapidly changing baselines and varying peak heights under different noise levels.

Bayesian approach of weight assignment

Generally, experimental data, denoted by yi, consist of three components: baseline zi, noise ϵi, and possible signal si: $y_i=s_i+z_i+{ϵ}_i\mathrm{.}$ (10) The baseline and noise are commonly present throughout the spectrum, whereas signals are only present in certain parts. Therefore, some data points in different spectral regions have signal contributions, whereas others do not. As an analogy to the good-and-bad data model [18], the likelihood $\text{prob}(y_i|\beta \mathrm{,}\mathcal{I})$ of the data point i can be considered as a weighted combination of case Bi, which contains no signal but only baseline and noise, and case ${\overline{\text{B}}}_i$, which contains also a signal, as follows: $\text{prob}\left(y_i|\beta \mathrm{,}\mathcal{I}\right)=\beta \text{ prob}\left(y_i|{\overline{\text{B}}}_i\mathrm{,}\mathcal{I}\right)+(1-\beta )\text{ prob}\left(y_i|{\text{B}}_i\mathrm{,}\mathcal{I}\right)\mathrm{,}$ (11) where 0 < β < 1 and the nuisance parameter β represents the probability that a datum contains a signal component. All relevant information $\mathcal{I}$ includes the nature of the physical situation and knowledge of the experiment. The prior probability of case Bi is $\text{prob}\left({\text{B}}_i|\beta \mathrm{,}\mathcal{I}\right)=1-\beta \mathrm{.}$ (12) The task is to obtain the appropriate weight using Eq. (5) for each point to achieve a smooth baseline zi. According to the Bayesian theorem [19], the posterior probability for case Bi containing no signal is $\text{prob}\left({\text{B}}_i|y_i\mathrm{,}\mathcal{I}\right)=\frac{\text{prob}\left(y_i|{\text{B}}_i\mathrm{,}\beta \mathrm{,}\mathcal{I}\right)\text{ prob}\left({\text{B}}_i|\beta \mathrm{,}\mathcal{I}\right)}{\text{prob}\left(y_i|\beta \mathrm{,}\mathcal{I}\right)}\mathrm{,}$ (13) where $\text{prob}\left(y_i|{\text{B}}_i\mathrm{,}\beta \mathrm{,}\mathcal{I}\right)$ is the likelihood function of the data point i in the absence of a signal.

Starting with the common assumption of Gaussian noise with zero mean and standard deviation σ, the likelihood function of case Bi containing no signal can be approximated by a Gaussian probability density function: $\text{prob}\left(y_i|{\text{B}}_i\mathrm{,}{z}_i\mathrm{,}\sigma \mathrm{,}\mathcal{I}\right)=\frac{1}{\sqrt{2\pi }\sigma }\text{exp}\left[-\frac{{(y_i-z_i)}^2}{2{\sigma }^2}\right]\mathrm{.}$ (14) When the datum contains a signal contribution si, the likelihood of ${\overline{\text{B}}}_i$ is given by $\text{prob}\left(y_i|{\overline{\text{B}}}_i\mathrm{,}{z}_i\mathrm{,}{s}_i\mathrm{,}\sigma \mathrm{,}\mathcal{I}\right)=\frac{1}{\sqrt{2\pi }\sigma }\text{exp}\left[-\frac{{(y_i-z_i-s_i)}^2}{2{\sigma }^2}\right]\mathrm{.}$ (15) Because si is unknown before the analysis, it must be marginalized by $\begin{array}{c}\text{prob}\left(y_i|{\overline{\text{B}}}_i\mathrm{,}{z}_i\mathrm{,}\sigma \mathrm{,}\mathcal{I}\right)={\displaystyle {\int }_0^{\infty }}d{s}_i\text{ prob}\left(y_i|{\overline{\text{B}}}_i\mathrm{,}{z}_i\mathrm{,}{s}_i\mathrm{,}\sigma \mathrm{,}\mathcal{I}\right)\\ \text{ prob}\left(s_i|\mathcal{I}\right)\mathrm{.}\end{array}$ (16) As a convention, si is always positive in the task of baseline correction. Here, a common scale factor μ is introduced to describe the signal; let μ be the mean value of all signal heights in the spectrum. The maximum entropy principle leads to a prior probability of si [19] under a Poisson distribution $\text{prob}\left(s_i|\mu \mathrm{,}\mathcal{I}\right)=\frac{1}{\mu }\text{exp}\left[-\frac{s_i}{\mu }\right]\mathrm{.}$ (17) Substituting Eqs. (15) and (17) into Eqs. (16) for si, we obtain the likelihood function of a data point containing signal as follows: $\begin{array}{l}\text{prob}\left(y_i|{\overline{\text{B}}}_i\mathrm{,}{z}_i\mathrm{,}\sigma \mathrm{,}\mu \mathrm{,}\mathcal{I}\right)\\ =\frac{1}{2\mu }\left(1+\text{erf}\left[\frac{y_i-z_i}{\sqrt{2}\sigma }-\frac{\sigma }{\sqrt{2}\mu }\right]\right)\\ \text{exp}\left[-\frac{y_i-z_i}{\mu }+\frac{{\sigma }^2}{2{\mu }^2}\right]\mathrm{.}\end{array}$ (18) Therefore, the likelihood function in Eq. (11) can be rewritten as follows: $\begin{array}{c}\text{prob}\left(y_i|z_i\mathrm{,}\beta \mathrm{,}\sigma \mathrm{,}\mu \mathrm{,}\mathcal{I}\right)\beta \text{ prob}\left(y_i|{\overline{\text{B}}}_i\mathrm{,}{z}_i\mathrm{,}\sigma \mathrm{,}\mu \mathrm{,}\mathcal{I}\right)\\ +\left(1-\beta \right)\text{}\text{prob}\left(y_i|{\text{B}}_i\mathrm{,}{z}_i\mathrm{,}\sigma \mathrm{,}\mu \mathrm{,}\mathcal{I}\right)\\ =\beta \frac{\text{exp}\left[-\frac{y_i-z_i}{\mu }+\frac{{\sigma }^2}{2{\mu }^2}\right]}{2\mu }\left(1+\text{erf}\left[\frac{y_i-z_i}{\sqrt{2}\sigma }-\frac{\sigma }{\sqrt{2}\mu }\right]\right)\\ +(1-\beta )\frac{\text{exp}\left[-\frac{{(y_i-z_i)}^2}{2{\sigma }^2}\right]}{\sqrt{2\pi }\sigma }.\end{array}$ (19) Finally, by substituting the likelihood expressed by Eqs. (19) and the prior probability given by Eq. (12) into Eqs. (13), we obtain the posterior probability for a data point i containing no signal: $\begin{array}{l}\text{prob}\left({\text{B}}_i|y_i\mathrm{,}{z}_i\mathrm{,}\beta \mathrm{,}\sigma \mathrm{,}\mu \mathrm{,}\mathcal{I}\right)\\ =\frac{\text{prob}\left(y_i|{\text{B}}_i\mathrm{,}{z}_i\mathrm{,}\beta \mathrm{,}\sigma \mathrm{,}\mu \mathrm{,}\mathcal{I}\right)(1-\beta )}{p\left(y_i|z_i\mathrm{,}\beta \mathrm{,}\sigma \mathrm{,}\mu \mathrm{,}\mathcal{I}\right)}\\ \frac{(1-\beta )\text{ prob}\left(y_i|{\text{B}}_i\mathrm{,}{z}_i\mathrm{,}\sigma \mathrm{,}\mathcal{I}\right)}{\beta \text{}\text{prob}\left(y_i|{\overline{\text{B}}}_i\mathrm{,}{z}_i\mathrm{,}\sigma \mathrm{,}\mu \mathrm{,}\mathcal{I}\right)+(1-\beta )\text{ prob}\left(y_i|{\text{B}}_i\mathrm{,}{z}_i\mathrm{,}\sigma \mathrm{,}\mathcal{I}\right)}\\ \frac{1}{1+f(y_i\mathrm{,}{z}_i\mathrm{,}\beta \mathrm{,}\sigma \mathrm{,}\mu )},\end{array}$ where $\begin{array}{l}f(y_i\mathrm{,}{z}_i\mathrm{,}\beta \mathrm{,}\sigma \mathrm{,}\mu )=\frac{\beta }{1-\beta }\sqrt{\frac{\pi }{2}}\frac{\sigma }{\mu }\\ \left(1+\text{erf}\left[\frac{y_i-z_i}{\sqrt{2}\sigma }-\frac{\sigma }{\sqrt{2}\mu }\right]\right)\text{exp}\left[{\left(\frac{y_i-z_i}{\sqrt{2}\sigma }-\frac{\sigma }{\sqrt{2}{\mu }_i}\right)}^2\right]\mathrm{.}\end{array}$ (20)

Bayesian reweighted penalized least squares algorithm

Thus far, we have obtained the posterior probability of a data point containing no signal, expressed by Eq. (20). In this study, this probability value is used as the weight for each point in Eq. (5). Thus, the weight of each datum i can be expressed in a unified format: $\begin{array}{l}{w}_i=\frac{1}{1+\frac{\beta }{1-\beta }\sqrt{\frac{\pi }{2}}F(y_i-z_i\mathrm{,}\sigma \mathrm{,}\mu )}\\ \text{where }F(d\mathrm{,}\sigma \mathrm{,}\mu )=\frac{\sigma }{\mu }\left(1+\text{erf}\left[\frac{d}{\sqrt{2}\sigma }-\frac{\sigma }{\sqrt{2}\mu }\right]\right)\\ \text{exp}\left[{\left(\frac{d}{\sqrt{2}\sigma }-\frac{\sigma }{\sqrt{2}\mu }\right)}^2\right]\end{array}$ (21) In practice, nuisance parameters such as the signal chance β, signal scale μ, and Gaussian noise standard deviation σ for each data point must be properly assigned to use the algorithm. Similar to existing algorithms for PLS-based baseline estimation, we also divide the values of subtraction d = y - z into a positive part d₊ and negative part d_- in the iterations of our algorithm. Then, an approximation of σ and μ can be obtained using different parts of the subtraction: the standard deviation σ can be estimated from the quadratic mean of d_-, whereas the signal scale μ can be estimated from the mean value of d₊.

As mentioned in the previous section, β reflects the probability of signal contribution in the data. According to Eq. (21), if the signal-to-noise ratio (SNR) μ² / σ² [20] and signal scale μ are constant, the calculated weight wi can be expressed as a function of the subtraction di for different values of β, as shown in Fig. 1. The weight functions wi corresponding to different β values are clearly different.

Fig. 1Full-screen View

(Color online) Weight function wi, defined in Eq. (21), as a function of the subtraction di = yi - zi for different values of β. The SNR parameter ${\mu }^2/{\sigma }^2=100$ and signal scale μ=0.1 are fixed. The region between grey vertical lines highlights a notable change in wi when β approaches one.

Thus, a rational value for β is required. Given that the proportion and distribution of signals in the data are unknown, it is convenient to use the average probability of no signal appearance in the total data set as a uniform β factor assigned to each data point, as suggested by the concept of good-and-bad model [18]. In this study, the average probability of no signals in the total data set can be defined as follows: $\langle \text{prob}\left({\text{B}}_i|y_i\mathrm{,}\sigma \mathrm{,}\mathcal{I}\right)\rangle \equiv \frac{1}{N}{\displaystyle \sum _{i=1}^N}\text{prob}\left({\text{B}}_i|y_i\mathrm{,}\sigma \mathrm{,}\mathcal{I}\right)\mathrm{.}$ (22) Note that β is updated in each iteration and the iteration terminates when β approaches $1-\langle \text{prob}(B_i|y_i\mathrm{,}\sigma \mathrm{,}\mathcal{I})\rangle$. All in all, the algorithm flow of the Bayesian reweighted penalized least squares (BrPLS) method is shown in Algorithm 1.

Reliability tests

We used synthetic spectra to test the reliability of different PLS-based baseline-estimation algorithms. Synthetic spectra were generated by adding signal components and random noise to a pre-defined baseline as follows: $d(x)=s(x)+b(x)+n(x)\mathrm{,}$ (23) where s(x), b(x), n(x), and d(x) denote the signal, baseline, noise, and resulting synthetic data, respectively. The synthetic signals were eight Gaussian peaks whose intensities, positions, and widths are listed in Table 1. Typical baseline patterns in PIXE [21-23] and Schottky spectra [24] contain linear, concave, or convex regions. Therefore, we chose two types of baseline for the reliability tests: linear and sinusoidal. Random noise was generated by using a random number from a Gaussian distribution. The SNR in terms of energy was defined in the absence of a baseline, according to the following equation [20]: $SNR=10{\text{log}}_{10}(E_{\text{s}}/E_{\text{n}})\mathrm{,}$ (24) where E_s refers to the expected energy value of the signal and E_n refers to the expected noise energy value in the spectrum. The SNR of the low-noise spectrum was set to 40 dB and that of the high-noise spectrum was set to 20 dB.

The maximum number of iterations was set to 50 for the AsSL, airPLS, arPLS, and BrPLS methods. The early termination criterion for AsSL, arPLS, and BrPLS was a standard deviation of the difference between two estimated baselines in nearby iterations less than 10^-6. For airPLS, Eq. (8) was used as the early termination criterion.

To start the tests of these four methods, an appropriate smoothness parameter λ in Eq. 3 should be determined. Recall that a larger λ leads to a smoother estimated baseline; if λ is too large or too small, the estimated baseline will not match the curvature of the baseline or will be highly distorted by high peaks. In addition, the sensitivities of the four methods to the same λ are different. To obtain the best λ for each method, the parameter space of λ was scanned from 10² to 10¹² in 10-base logarithmic steps of 0.1. By the end of the iterations, the algorithm reaches a final set of zi; then, the root mean square error (RMSE) can be calculated as follows: $\text{RMSE}(\lambda )=\sqrt{\frac{{{\displaystyle {\sum }_{i=1}^N\left(y_i-z_i\right)}}^2}{N}}\mathrm{.}$ (25) This RMSE expression reflects the quadratic mean of the estimated signal. If the baseline shape is linear, as λ increases, the RMSE curve first increases rapidly and then converges, as shown in Figs. 2a, c. In this case, λ values can be chosen at any location where the RMSE is stable. For example, the stable RMSE(λ) region shown in Fig. 2a corresponds to values of λ ranging from 10^8.1 to 10¹² and that shown in Fig. 2c corresponds to values of λ ranging from 10⁸ to 10¹². If the baseline contains depression or protrusion regions, the RMSE as a function of λ can become more complicated, as shown in Figs. 2b, d. A suitable λ value is always chosen at the stable plateau of RMSE(λ), where the estimated baseline is close to the actual baseline. For the specific case of a sinusoidal baseline, as in Fig. 2b, the optimal values of λ that can be selected range from 10^6.5 to 10⁸ for the AsLS and airPLS methods, from 10^6.5 to 10^10.2 for the arPLS method, and from 10^6.5 to 10^9.9 for the BrPLS method. Similarly, for the case shown in Fig. 2d, the optimal λ values that can be selected range from 10⁷ to 10^8.6 for the AsLS and airPLS methods, and from 10^4.5 to 10^8.6 for the arPLS and BrPLS methods. Generally, the expected baseline must be as smooth as possible. Therefore, the maximum optional λ was used in the tests for each method, as listed in Tab. 2.

Fig. 2Full-screen View

(Color online) RMSE as a function of the λ parameter calculated from baseline estimation results. Figs. a and c depict the case of linear synthetic baseline, while Figs. b and d depict the case of sinusoidal synthetic baseline. The SNRs used for generating the synthetic spectra were 20 dB in Figs. a and b and 40 dB in Figs. c and d. The green, red, blue, and orange dash-dot lines represent the results of the AsLS, airPLS, arPLS, and BrPLS methods, respectively. The grey lines indicate the location of the RMSE plateau.

To quantitatively assess the performance of each method, we introduced the goodness of the baseline estimation factor R² [25]: $R^2=1-\frac{{\displaystyle {\sum }_{i=1}^N{\left(b_i-z_i\right)}^2}}{{\displaystyle {\sum }_{i=1}^N{\left(\langle b\rangle -b_i\right)}^2}}\mathrm{,}$ (26) where bi and $\langle b\rangle$ denote the synthetic baseline and mean of all the baseline values, respectively. Ideally, the estimated baseline values match the synthetic values, i.e., R² = 1.

The results are presented in Fig. 3 and Table 2. The BrPLS and arPLS methods have R² values close to 1, but the R² values of the other two methods deviate from 1. The results suggest that the AsSL and airPLS methods are inadequate for spectral-baseline estimation. The differences in R² between the BrPLS and arPLS methods are less than 0.5%, indicating no significant differences between both methods in these analysis tests. In Figs. 3(c) and 3(d), the estimated baseline near the data boundary is slightly distorted within the 1-σ error interval of the noise owing to the limitations of the PLS-based methods, where only the 2nd-order difference calculation (m = 2) of the baseline was used in Eqs. (1) and (2). In any case, the differences between the estimated and synthetic baselines of both the arPLS and BrPLS methods are within the 1-σ interval of the synthetic Gaussian noise. This means that the estimated baselines of both methods are reliable within the 1-σ error interval. Note also that the BrPLS method has a relatively better ability to handle the signal peaks in the convex regions of the baseline, as indicated by the arrows in Fig. 3(d). This will be further evaluated in the next section.

Fig. 3Full-screen View

(Color online) Results of baseline estimation of synthetic spectra. The solid grey lines present the synthetic spectra. The top graph shows the estimated baselines for the four methods and all the original data. The graph below shows the difference between the estimated and synthetic baselines. The interval from -σ to σ of the Gaussian noise in the y-axis is highlighted in grey. The green, red, blue, and orange dash-dot lines are estimated baselines using the AsLS, airPLS, arPLS, and BrPLS methods, respectively. The synthetic spectra depicted in Figs. (a, c) were generated with SNR = 20 dB, whereas those shown in Figs. (b, d) were generated with SNR = 40 dB. The parts marked with arrows show that the BrPLS method outperforms the arPLS method when dealing with the condition in which signal peaks are located on the baseline protrusion regions.

Comparison between the arPLS and BrPLS methods

According to the results presented in the previous section, the BrPLS and arPLS methods outperform the other two methods regardless of the baseline shape. Moreover, the BrPLS method has a slight advantage over the arPLS method when dealing with cases in which the signal peaks are located in the convex regions of the baseline. In this section, we discuss and compare the advantages and limitations of both algorithms by applying them to the analysis of other synthetic spectra.

The difference between the BrPLS and arPLS methods increases in regions where the baseline has protruding parts. To estimate such a baseline, the most challenging locations are typically those at which the signal peaks lie near the highest point of the convex part and at the inflection points of the rising and falling edges. Thus, a typical convoluted Landau-Gaussian function was used to generate the synthetic baseline in the synthetic spectrum, and three Gaussian signal peaks were placed in the aforementioned parts of the baselines. The intensities and widths of the signal peaks are listed in Table 3. The SNR was set to 40 dB. The maximum number of iterations and early termination criteria for both methods were the same as those in the previous tests. Similarly, the λ values were selected in the same manner as in the previous tests.

The results are presented in Fig. 4 and Table 3. Both methods produce reasonable baseline results. However, when the signal peaks are located at the rising edge of the baseline, the two methods yield unsatisfactory estimations. This is another inherent limitation of PLS-based methods as they strive to balance the smoothness of the baseline and fidelity to the original spectrum. Signal peaks located at the rapidly changing parts of the baseline represent the case in which the two factors are conflicting. As the peak width and height increase, the baseline obtained by the arPLS method further deviates from the actual value when the signal peak is located at the baseline protrusion regions. In contrast, the BrPLS method yields significantly better results, as shown in Fig. 4(i). This is because the BrPLS method considers the heights of different signal peaks and assigns different weight values to the parameter μ in Eq. (21). Conversely, the weight function in Eq. (9) for the arPLS method is a symmetric sigmoid curve that does not consider signal peak characteristics.

Fig. 4Full-screen View

(Color online) Baseline estimation results using the arPLS (blue dash-dot lines) and BrPLS (dash orange lines) methods. Nine simulation cases, from (a) to (i), share the same convoluted Landau-Gaussian baseline under the same SNR of 40 dB. Additional information for the nine simulations is listed in Table 3. The solid grey lines represent the synthetic spectra. The top graph in each case shows the estimated baselines and original data. The middle graph shows the difference between the estimated and synthetic baselines. The intervals in the y-axis, highlighted in light grey, represent the -3σ-to-3σ intervals of the Gaussian noise. The bottom graph shows the estimated signal magnitudes after baseline subtraction.

Application to real-world Schottky spectra

Schottky spectroscopy [26] is an important method for beam diagnostics in heavy-ion accelerator facilities. It is routinely used in precision mass and lifetime measurement experiments in nuclear physics studies [27]. The corresponding Schottky spectrum reflects the power density of signals induced by circulating ions in a storage ring. The periodic motion of an ion is represented by the resonance peaks at each harmonic of the revolution frequency in the measured spectrum [28]. Accurate determination of the peak position and area in the power density spectrum is essential for data analysis. The revolution frequency, which is proportional to the mass-to-charge ratio of the ion, can be derived from the peak position. The peak area can be derived from the baseline-corrected spectrum. It is proportional to the number of corresponding ions and can be used to determine the lifetime of the ions in the storage ring.

In the Experimental Storage Ring at Darmstadt [29] and experimental Cooler Storage Ring (CSRe) at Lanzhou [30], ions circulate with frequencies ranging from several hundred kHz to a few MHz. The bandwidth of the measured Schottky spectrum is typically less than 2 MHz [31-33]. The baseline of this type of spectrum can be approximated as a straight line and, in some cases, can be corrected using simple 1st-order polynomial fittings [34, 35]. However, the baseline becomes complicated for the full spectrum, which contains several harmonics of the full momentum acceptance of the storage ring. Owing to complex signal processing in hardware [24], the baseline appearing in this spectrum does not always have a theoretical Lorentzian shape [36]. Different types of baseline distortions may occur during the actual signal processing. Thus, a method that can effectively estimate and correct Lorentzian-type baselines is required. The newly developed BrPLS method offers a new option in this direction.

The power spectral density of Schottky noise is multiplicative [37] and our Bayesian assumption is based on the additive components in the data, as expressed in Eq. (10). Therefore, a logarithmic transformation of the power densities in the original Schottky spectra was applied first. The processed spectra were then adjusted to non-negative values by subtracting the minimum values. Experimental data were obtained from the measurements of ⁸⁶Kr projectile fragments [38-40] at CSRe in 2018. Each frame in the Schottky spectrum had a data acquisition time of 2.73 s. The statistical uncertainty of the amplitude at a given frequency point was approximately 1.4% [41]. Particle identification in the spectrum was performed according to [26]. During baseline estimation, steps similar to those of the previously described synthetic data treatment were applied. The parameters and other settings of the arPLS and BrPLS methods were the same as those described in the previous section. The results of a simple baseline fit with a local 1st-order polynomial (fitting frequency range between ± 30 kHz of the signal peak center) are also presented for comparison. The peak position and area under the frequency peak were determined by fitting Gaussian functions to the peaks after baseline corrections.

The baseline estimation results for a single frame are shown in Fig. 5. Consistent with the previous section, the BrPLS method could effectively handle the presence of peaks in the baseline protrusion regions. In contrast, the arPLS method underestimated the baseline to the extent of introducing more spurious signal peaks, as indicated by the blue line in Fig. 6. In Fig. 5, five signal peaks are marked for further discussion. Peaks 1 and 4 belong to the same ion species ⁸²Se³⁴⁺ but differ in the harmonic number of the ion revolution frequency. Peak 1 is located in the flat baseline part whereas peak 4 is located in the middle of the falling edge of the baseline. Peaks 2 and 5 belong to the same ion species ⁸⁴Br³⁵⁺ and are located in the relatively flat part of the baseline. The small ion peak 3 belongs to ⁷⁴Ga³¹⁺ ions and is located near the top of the baseline-protrusion region. The estimated revolution frequencies were consistent among the three methods within the 1σ error range, as shown in Fig. 7. The central values vary slightly as a function of the peak position and method, particularly when the peak is located in the convex part of the baseline. The peak area as a function of the elapsed time was used to determine the ion storage lifetime after ion injection into the ring. More spectra shown in the top graph of Fig. 5 were analyzed for this task, and the differences in the performance of the three methods became more apparent. Ideally, the ion lifetimes estimated from the ion peaks at different harmonic numbers were expected to be the same. As shown in Fig. 8, the three methods achieved this goal. However, arPLS incurred a relatively larger error for ion peak 4, which is located in the baseline protrusion region. In this region, the ion peak areas obtained by the arPLS method were somewhat overfitted, resulting in an underestimation of the baseline. Some blue points obtained using this method and depicted in Fig. 8(b) deviate from the overall trend, leading to larger uncertainty in the corresponding lifetime, as shown in Fig. 8(c). The overall difference between the results of the polynomial fitting and BrPLS methods is negligible. It is expedient to use a local linear fit to the baseline to estimate the area of individual peaks when the baseline shape of the full spectrum is very complex. However, this method relies heavily on the degree of linearity of the baseline within the selected local range. Additionally, to estimate the neutral atomic lifetime of the corresponding nuclide of an ion from its storage lifetime in the ring, it is necessary to use other stored ions in the same spectrum as a reference. A local linear fit to individual peaks is challenging to guarantee the consistency of the baseline estimate for all the peaks located in different parts of the spectrum. In contrast, the results obtained by the newly developed BrPLS method are generally more stable and reliable than those of the other two methods.

Fig. 5Full-screen View

(Color online) Schottky spectra experimentally measured and corresponding baseline estimation results. The top graph shows the 2D power spectral density in the frequency domain as a function of time after ion injection. The time resolution is 2.73 s. The middle graph shows one frame of spectrum measured at t = 25 s. The power spectrum after baseline correction is presented in the bottom graph. Five ion peaks are marked for further discussion. Ion peaks belonging to the same ion species but with different harmonic numbers of the revolution frequency are identified by the same color box: peaks 1 and 4 belong to ⁸²Se³⁴⁺, peaks 2 and 5 belong to ⁸⁴Br³⁵⁺, and peak 3 belongs to ⁷⁴Ga³¹⁺.

Fig. 6Full-screen View

(Color online) Baseline-corrected spectrum, as shown in the bottom graph of Fig. 5, zoomed into the frequency range of (-500, 500) kHz.

Fig. 7Full-screen View

(Color online) Revolution frequency determined from the baseline-corrected spectrum. The peak positions were extracted from Gaussian peak fitting. The revolution frequencies were normalized by the corresponding harmonic numbers. The labels of the ion peaks are the same as in Fig. 5.

Fig. 8Full-screen View

(Color online) Lifetime estimation results of ⁸²Se³⁴⁺(neutral atom half-life T_1/2 = 87.6 Ey [42]) using data (peaks 1 and 4) in Fig. 5. The left graphs (a) and (b) show peak areas, normalized to the initial area after baseline correction, as a function of the time. The solid lines are exponential-function fitting results. The labels of the ion peaks are the same as those in Fig. 5. The error bars represent statistical uncertainties of the signal areas extracted from the Gaussian fits. The deduced storage half-lives of the corresponding ions in CSRe are shown in graph (c).