An algorithm for the automatic deglitching of X-ray absorption spectroscopy data

Analysis of X-ray absorption spectroscopy data often involves the removal of artifacts or glitches from the acquired signal, a process commonly known as deglitching. Glitches result either from specific orientations of monochromator crystals or from scattering by crystallites in the sample itself. Since the precise energy – or wavelength – location and the intensity of glitches in a spectrum cannot always be predicted, deglitching is often performed on a per spectrum basis by the analyst. Some routines have been proposed, but they are prone to arbitrary selection of spectral artifacts and are often inadequate for processing large data sets. Here, a statistically robust algorithm, implemented as a Python program, for the automatic detection and removal of glitches that can be applied to a large number of spectra, is presented. It uses a Savitzky–Golay filter to smooth spectra and the generalized extreme Studentized deviate test to identify outliers. Robust, repeatable, and selective removal of glitches is achieved using this algorithm.

1. Introduction

In X-ray absorption spectroscopy (XAS), glitches correspond to artifacts in a confined energy range, usually a few data points, manifesting as a sharp variation in measured absorption (Stern & Lu, 1982; Bunker, 2010; Calvin, 2013). Glitches often arise from multiple-diffraction events within the crystal monochromator used to set the energy of the X-rays, resulting in a significant decrease or rise in the intensity of the beam delivered to the sample (Bauchspiess & Crozier, 1984). Alternatively, diffraction from other crystalline phases, for instance in diamond anvil cells, may also create glitches, impacting the measurement of the absorption coefficient [μ(E)] (Sapelkin & Bayliss, 2002).

Analysts can minimize the impact of glitches by using best practices at the time of data collection, which includes ensuring uniform sample preparation, confirming the linearity of equipment at the beamline, and detuning the monochromator to reduce the intensity of the multiple diffraction events; however, these practices do not completely eliminate glitches (Abe et al., 2018). Collecting data free from significant artifacts is particularly challenging in some energy ranges, such as at the Fe K-edge using a Si(111) monochromator (Pickering, 1999). Leaving spurious points in place may interfere with the analysis of XAS data; for instance, glitches may introduce error either in the normalization of the spectrum or the transformation of the extended X-ray absorption fine structure (EXAFS) data into R-space (Abe et al., 2018), and spectra may be improperly clustered. As such, the processing of spectra may require a deglitching step, where data points corresponding to glitches are removed.

Given the erratic nature of glitches, their removal is commonly based on the judgment of the analyst through visual inspection of the data. For example, the program Athena offers a graphical user interface for glitch removal, where one may either conduct point-and-click identification of spurious data points or set a threshold value based on the post-edge normalization curve outside of which data points are removed (Ravel, 2016). Other programs include the capability to remove data points at specified indices (Wellenreuther & Meyer-Klaucke, 2009) or the option to compare μ(E) with respect to the incoming incident beam intensity I₀ (Aberdam, 1998). As a result, these methods require the analyst to manually inspect each XAS spectrum to remove the aberrant values, an approach that is readily applicable to large glitches and small datasets.

However, the growing relevance of time-resolved XAS (Bak et al., 2018) and continuous scanning `quick EXAFS' modes (Prestipino et al., 2011) along with increasing accessibility to laboratory-source XAS instruments (Anklamm et al., 2014; Jahrman et al., 2019) promises large datasets where manual deglitching of spectra becomes impractical. A rapid and robust method for automated deglitching of XAS spectra is needed to address this shortcoming. Previous methods with greater potential for automation used the first derivative of the absorption coefficient to identify glitches (Zhuchkov et al., 2001). While useful, such strategies can result in points adjacent to glitches being erroneously identified, and the applicability to the X-ray absorption near-edge structure (XANES) region is limited. Moreover, a manual threshold must still be set by the analyst to identify glitches based on visual inspection of the data.

Here, we propose a method for the high-throughput deglitching of XAS data. Outliers are identified in the normalized residuals between original and low-pass-filtered data. A second iteration of this process minimizes the occurrence of type I `false positive' errors. With this method, we achieve accurate and repeatable removal of glitches from full-spectrum XAS data. We present deglitching examples on XANES and EXAFS spectra, as well as a negative control. Finally, we discuss potential limitations of our method as well as strategies for improvement.

2. Methods

The deglitching algorithm was implemented as a computer program written in the Python 3.7 programming language. The packages Scipy (Virtanen et al., 2020) and Numpy (Oliphant, 2006) were used for all calculations, and Larch (Newville, 2013) was used for XAS data processing in a Jupyter Notebook environment (Kluyver et al., 2016).

The deglitching program was designed to be compatible with the data structure used by Larch for XAS analysis. The only required input for the deglitching program are data channels corresponding to the energy (E) and the absorption coefficient [μ(E)] within a Larch group. Absorption data may or may not be normalized prior to deglitching. Additional parameters may be specified here to tune the deglitching program, though default parameters should work in most circumstances.

All spectra presented here were collected at the bending magnet beamline of the Dow–Northwestern–Dupont Collaborative Access Team (DND-CAT), Sector 5 at the Advanced Photon Source at Argonne National Laboratory, Lemont, IL, USA. Fe K-edge data are presented for drinking water treatment residual samples collected from drinking water treatment plants in the USA. The other samples include a denture adhesive cream and cobalt (II, III) oxide (Aldrich). Energy was set using a Si(111) double-crystal monochromator. Intensity values for the incident beam (I₀), the transmitted beam (I_T1), and the secondary transmitted beam (I_T2) were collected using Oxford ionization chambers with 29.6 cm path lengths. Fluorescence data were captured using Vortex ME4 silicon drift detectors. XAS data were collected at an energy interval of 10 eV from 150 to 20 eV below the edge energy, 0.5 eV in the XANES region (up to k = 3 Å⁻¹), and 0.05 Å⁻¹ in the EXAFS region.

3. Algorithm

In the following a description of the deglitching process, also depicted in Fig. 1, is given. After data channels for energy and μ(E) are provided to the algorithm (Fig. 1, box 1), the μ(E) channel is fit with a Savitzky–Golay filter to provide a smoothed representation of the data, μ_SG(E) (Fig. 1, box 2 and subplot A). The Savitzky–Golay filter uses a least-squares fitting procedure to fit each point over a rolling window of odd length (w_SG), effectively acting as a low-pass filter (Savitzky & Golay, 1964) (Fig. 2). Savitzky–Golay window length and the polynomial order of the rolling fit are adjustable parameters set at 9 and 3, respectively, by default. The results of the filter are then subtracted from the normalized absorption to compute the residuals between the original and filtered data, δμ(E)_A (Fig. 1, box 3),

Given that the Savitzky–Golay filter acts as a low-pass filter, δμ(E)_A will reflect the contribution of high frequencies to μ(E). One expects a normal distribution of residuals should Gaussian noise be the only source of misfit between the data and a fit aiming to represent the true values of the data (Trutna et al., 2003). Glitches, being aberrant high-frequency spectral features, will result in residuals that are outliers compared with the rest of the residuals. In simple cases, only glitches will be outliers in these residuals. However, in the absence of a glitch, these residuals will vary slightly based on regional features of the XAS spectrum. One can generally expect larger values for |δμ(E)_A| near the absorption threshold, owing to the wider range of absorbances contained in a fitting window within this region (Fig. 3) and, potentially, the attenuation of some high-frequency features. Identifying outliers on untreated residuals may result in the false positive identification of glitches in the XANES region or failure to identify subtle glitches in the EXAFS region. In addition, low-pass filters cannot completely attenuate high-frequency signals. This is apparent in Fig. 1, box 2 and subplot C: the initial Savitzky–Golay filter is not smooth in the region surrounding aberrant points. Glitches will impact δμ(E)_A at any point within (w_SG/2) − 1 points of a glitch, increasing the variance of these residuals relative to those outside the Savitzky–Golay filter window.

Figure 1 The deglitching procedure, outlined textually (left) and graphically (right). Numbers and colors in the text outline correspond with colors and labels in the graphical outline. Subplot A shows the μ(E) data, the filtered data, and initial and normalized residuals. Subplot B shows the histogram for the first set of residuals normalized to the rolling median absolute deviation (MAD) (Step 4). The red lines show the threshold for outlier values, and the red dots are used to draw attention to the outlying points. Subplot C shows the interpolation of candidate points and compares the initial and final filters with their respective non-normalized residuals. Finally, subplot D shows the resulting EXAFS from the deglitching algorithm.

Figure 2 An example of the rolling window used for the Savitzky–Golay filter. Data shown are the same as the data in Fig. 1. A point centered in a rolling window of odd length (here, nine points) undergoes a polynomial fit using a least-squares approach.

Figure 3 Comparison between regionally normalized and non-regionally normalized residuals between the μ(E) data and Savitzky–Golay filter. Residuals shown the same example as Steps 3 and 4 in Fig. 1. Regionally normalized residuals are divided by a rolling-window median absolute deviation (MAD), and non-regionally normalized residuals are divided by the MAD for the full spectrum. Residuals are plotted on a point-index basis, and the top axis shows how this translates into energy. Minor tick marks are placed every 50 eV, and gray vertical lines indicate changes in the sampling interval (10 eV in the pre-edge region, 0.5 eV in the XANES region, and 0.05 Å^−1 in the EXAFS region). XANES residuals from a point index of approximately 50 to 200 are larger than non-glitch EXAFS residuals (such as those above a point index of 280, or 7600 eV) in non-normalized residuals. In contrast, normalized residuals are of a consistent magnitude in the absence of glitches.

Residuals must be normalized regionally to account for differences in residual values between XANES and EXAFS regions. Here, a rolling window approach, similar to the Savitzky–Golay filter, is adopted. Average-based values, like standard deviation, are strongly impacted by outliers, while median values are comparatively more resistant to outliers. For this reason, a rolling median is found for |δμ(E)_A| to acquire a rolling median absolute deviation. Much like the Savitzky–Golay filter, these median values are calculated from a rolling window, w_m. The window size for the median calculation is selected such that, should a glitch be contained within w_m, more than half the points in w_m will not have included the glitch in their Savitzky–Golay filter fitting window w_SG. This may be defined as

where w_m is the window length for calculating the median absolute deviation, w_SG is the Savitzky–Golay filter window length, and L_g is the maximum number of points corresponding to a single glitch. Normalized residuals r are computed by dividing δμ(E)_A by the local median absolute deviation of residuals to account for regional variability in the fit,

where r_i corresponds to the normalized residual value at E_i. At the beginning and end of the data, calculations are performed using truncated windows. Provided Savitzky–Golay parameters that do not filter true signal, these normalized residuals will be normally distributed with a standard deviation of approximately 1.48, corresponding to the conversion from median absolute deviation to standard deviation based on the cumulative distribution function of the standard normal distribution (Filliben & Heckert, 2003) (Fig. 1, box 4 and subplot B). This may be confirmed by plotting the normalized residuals and a histogram of their distribution; well defined oscillations in the residuals are indicative of an overlong Savitzky–Golay filter window length; and non-normal distributions, such as bimodal distributions, are an indication of poor normalization due to the selection of the filter window length (w_SG) and/or the maximum glitch length (L_g).

From the normalized residuals, outliers are mathematically identified using a generalized extreme Studentized deviate test (generalized e.s.d.) (Fig. 1, box 5 and subplot B). The generalized e.s.d. identifies outliers in normally distributed data when provided with a maximum number of outliers and a significance value for the outlier identification (Rosner, 1983; Czesla et al., 2019); these values are set at 10% of the data points contained in the spectrum and 0.025 by default, respectively. The first step of the generalized e.s.d. is to calculate the mean of the dataset. Next, the data point furthest from the mean value is identified. The test statistic is calculated by finding the distance of this point from the mean value in terms of number of standard deviations. Using a t distribution at the provided significance level, the maximum acceptable distance from the mean is calculated for a dataset of a given length. This provides a critical value which is, again, in units of number of standard deviations from the mean. The most distant point is removed and the calculation is repeated until the maximum number of outliers are reached. If the critical value of a given point is greater than the test statistic, that point and all previously analyzed points are identified as outliers (Rosner, 1983).

Outliers identified in this first stage are taken as candidate points for glitches. A sufficiently large glitch may result in a poor fit from the Savitzky–Golay filter on points within the same filter window as the glitch (Fig. 1, box 2 and subplot C), so outliers in the first pass (candidate points) may not all be glitches. In the data in Fig. 1, for instance, although there are only four obvious glitch points, seven points are identified as outliers (Fig. 1, box 5 and subplot B). An interpolation step is included in the deglitching procedure to account for the limitations of low-pass filtering (Fig. 1, subplot C); because high frequencies are dampened and not eliminated by low-pass filters (including the Savitzky–Golay filter), points near the glitch are poorly fit by the filter. This leads to the identification of three candidate points that, while near glitches, are not part of any glitches. Candidate points are removed from a copy of the data, and μ(E) at these points is interpolated using a cubic spline (Fig. 1, box 6). The resulting μ(E)_B only differs from μ(E) at the candidate points. A second set of Savitzky–Golay filtered data, μ_SG(E)_B, is generated based on this copy of the data using the same parameters as before (Fig. 1, box 7A). From here, the initial process is repeated, finding the residuals [δμ(E)_B] between the original data [μ(E)] and the second Savitzky–Golay filter [μ_SG(E)_B], defined as

As before, these residuals are normalized by the regional median absolute deviation, this time calculated on δμ(E)_B. The window for this second normalization shrinks to (2L_g) + 1 points, effectively limiting the maximum number of points in a glitch to L_g. This adjustment is made possible by the interpolation of candidate points, which minimizes the impact of glitches on the filter's fit for other points within w_SG. Outliers are identified in δμ(E)_B using the generalized e.s.d. The outliers identified in this second pass that are within one half of the Savitzky–Golay window length, w_SG, of a candidate point are taken to be glitches and removed. In our example, the interpolated points for the non-glitch candidate points essentially overlap with the original data (Fig. 1, box 6). As a result, the residuals between the second filter and these three candidate points are small, and these points are not identified as glitches.

4. Results and discussion

The deglitching algorithm was tested on a range of elements and data collection modes, as shown in Fig. 2. For all data, the deglitching procedure was applied across the full spectrum and used a Savitzky–Golay window length of 9, a significance level of 0.025, and a maximum glitch length of 4 (the default values for the program). In each of these circumstances, the algorithm successfully identified and removed glitches across the full spectrum of data without removing non-glitch points. Given the diverse datasets included in these examples – which include XANES glitches, minor EXAFS glitches, and data for various elements collected in both absorbance and fluorescence mode – the default parameters should be adequate for most applications wherein the data are sampled on a similar energy grid (0.5 eV in the XANES, 0.05 Å⁻¹ in the EXAFS).

Using the sampling regime of the example data, the default parameters would result in a Savitzky–Golay filtering window that spans 4.0 eV in the XANES and 0.4 Å⁻¹ in the EXAFS. Initial normalization of the residuals would occur over a 12.0 eV window in the XANES and a 1.2 Å⁻¹ window in the EXAFS, and normalization of the final residuals would cover 4.0 eV in the XANES and 0.4 Å⁻¹ in the EXAFS. For full-spectrum deglitching, these ranges should not be larger than the range of energies with high-derivative features (i.e. near the rising edge and the white line). This is especially important for the normalization of the final residuals, which is based on the maximum glitch length. In most cases, in the interest of maintaining more consistent rolling window lengths, L_g may be kept equal to (w_SG − 1)/2.

Figs. 4, 5, and 6 show the results of the deglitching algorithm as applied to distinct datasets. Fig. 4 shows an example of the algorithm simultaneously removing glitches in the XANES [Fig. 4(B)] and EXAFS [Fig. 4(C)] regions for data collected in fluorescence mode at the Fe K-edge. Fig. 5 shows data for a sample where the algorithm removed a subtle EXAFS glitch at a high K value in a sample of denture adhesive cream collected in transmission mode at the Zn K-edge. Finally, Fig. 6 presents a negative control for analysis, a cobalt (II–III) oxide sample collected in transmission mode with no apparent glitches, where the deglitching algorithm identified no points as glitches. All glitches identified and removed in these examples correspond to monochromator glitches, but not all monochromator glitches in these examples resulted in aberrant data points.

Figure 4 Results of the deglitching algorithm on Fe K-edge fluorescence data from a drinking water treatment residual sample. Subplot A shows the full spectrum of flattened, normalized XAS data; I0; and the energies of removed glitches. Subplot B shows the deglitched XANES data along with the removed glitch. Subplot C compares the original and deglitched EXAFS data. In total, five points, including one point in the XANES and four points in the EXAFS, are identified as glitches and removed.

Figure 5 Results of the deglitching algorithm on Zn K-edge absorbance data from a denture adhesive cream sample. Subplot A shows the full spectrum of flattened, normalized XAS data; I0; and the energies of removed glitches. Subplot B compares the original and deglitched EXAFS data. One point near 13 Å in K-space, corresponding with a subtle monochromator glitch, is removed.

Figure 6 Results of the deglitching algorithm on Co K-edge absorbance data from a cobalt (II–III) oxide sample. Subplot A shows the full spectrum of flattened, normalized XAS data and I0; subplot B shows the EXAFS data. No points are removed by the deglitching algorithm.

The deglitching algorithm is unlikely to require major adjustments when applied to data with sporadic glitches sampled on a similar energy grid as the provided examples. Extended glitches with regular, absorption-like features present a greater challenge for this algorithm, given the intended use case of full-spectrum deglitching consistent across the XANES and EXAFS regions. For more challenging data, several parameters of the algorithm may be adjusted: lengthening the Savitzky–Golay filter window will result in a more stringent low-pass filter, which may be useful either for intense and/or extended glitches or for data sampled on a finer energy grid; the significance level α for outlier identification may also be changed, where a higher value will result in more aggressive outlier identification, which may be useful for low-intensity glitches; and the algorithm may be limited to a specific region of interest, either specified in terms of energy (e.g. 7300–7800 eV) or as XANES or EXAFS regions (defined in the program as being divided at 150 eV above the absorption threshold).

This algorithm provides repeatable and robust glitch removal in typical XAS data across the full spectrum. An analyst may confirm that the algorithm is performing according to expectations on one data scan, then apply those settings to all similar data, permitting the rapid processing of large XAS datasets.

The deglitching algorithm is available for download at https://github.com/wallacesam/deglitching.