Noise reduction and peak detection in x-ray diffraction data by linear and nonlinear methods

Significant progress has been made in the development of low-pass filters for removing noise in visible–near-ultraviolet (UV/VIS) spectra.^1–8 These filters trade on the fact that information appears as point-to-point correlations, whereas noise is characterized by point-to-point fluctuations. Thus, in general, the low-order coefficients obtained in a Fourier transform contain the information and the high-order coefficients of the noise. Successful filters capitalize on this separation in reciprocal (Fourier) space (RS) to minimize or even eliminate noise while preserving information. In principle, the Fourier reconstructions of the results back to direct (spectral) space (DS) are ideally faithful, noise-free versions of the originals.

Filtering can be either linear or nonlinear. A linear filter is a reversible filter that operates by attenuating coefficients, ideally mainly in the white-noise region. In contrast, a nonlinear filter operates by replacing white-noise coefficients with “most-probable” values. Attenuation of high-order coefficients, termed apodization, requires compromises. If the attenuation is too abrupt, the reconstruction exhibits Gibbs oscillations (ringing.) If the attenuation is too gradual, a mix of noise leakage (type-I errors) and information loss (type-II errors) occurs. Of the many transfer (window) functions that have been proposed, the recent quantification of type-II errors⁷ in RS reduces the number to two: the binomial filter introduced by Marchant and Marmet³ in 1982 and the Gauss–Hermite filter introduced by Hoffman et al.⁴ in 2002.

The 1967 “maximum-entropy” spectral-deconvolution (sharpening) approach of Burg,⁹ designed to extract harmonics buried in noise in stationary time series, has recently been shown to contain a true maximum-entropy solution, where coefficients in the information-containing region remain unchanged, and coefficients in the white-noise region are replaced with most-probable, model-independent extrapolations of the low-order coefficients.^{10, 11} While extrapolation solves the apodization problem, the corrected maximum-entropy (CME) approach goes further: it also provides model-independent estimates of the locations (energies or in the present case angles) of the singularities that give rise to features in the data.¹² The question addressed here is whether the same noise-reduction, feature-detection, and feature-location advantages seen in UV/VIS applications can be realized in the analysis of x-ray diffraction data. The results indicate that they can.

The diffraction data that we analyze are a set of 3214 points from 15.0014° to 59.9834° equally spaced by 0.014°, as shown in Fig. 1. These data consist of localized structures separated by large regions containing no signal, similar to those seen in Raman scattering.¹³ This provides a complete selection of examples and, hence, serves as a useful guide for reducing noise and detecting weak features in other applications as well. We approach noise reduction in two ways: linearly, using Gauss–Hermite (GH) DS convolution, and nonlinearly, through maximum-entropy calculations. Both are described in detail below. Both are necessary: the CME method operates only on features, so it cannot be used in the empty stretches between them. Here, noise can be reduced only by linear methods. In many cases, linear filtering may be good enough. We show examples of both below.

The purpose of this section is to define procedures. Comments and explanations are kept to a minimum.

We assume that the data consist of an odd number N of real values $f j$⁠, $1≤j≤N$⁠, equally spaced in the relevant variable. The CME formalism requires N to be odd because the elements $R μ ν= R ν μ ∗$ of its two-dimensional characteristic matrix, defined below, must be complex conjugates when reflected about the main diagonal. Reflection about a main diagonal necessarily requires an odd number of points. This CME restriction also simplifies mathematics in general.

The reality condition $f j= f j ∗$ requires $F − n= F n ∗.$

From now on, the Fourier coefficients of the data are assumed to be $A ′ n$ and $B ′ n$⁠, so the prime notation is dropped.

Thus, increasing N decreases the uncertainties in the Fourier coefficients. This counterintuitive result is easily demonstrated given the large dataset here.

Because $| F ¯ n | 2$ typically decrease exponentially, the least disruptive filters are maximally flat, that is, $H n=1$ to as large a value of n as practical. This is the Butterworth criterion,¹⁵ which states that as many low-order derivatives as possible should vanish in a Taylor-series expansion of $H n→H(n)$ about $n=0$⁠.

Equation (17) suggests that the optimum filter sets $H n=1$ to a cutoff index $n= n c$ and zero beyond. However, the DS coefficients of this “brick-wall” filter are given by the sinc function. Upon DS convolution, phase reversals and slow decay generate large oscillations when operating on any localized structure.

The associated $h j$ can be represented analytically, but a more practical approach is simply to obtain them using Eq. (14).

The results are approximations because the main goal of the procedure is to establish the most-probable line shapes. As a result, some roots perform only a “cosmetic” function and do not represent true singularities. These typically have a relatively large values $Γ ν$⁠. Examples are given below.

This approach also offers a path⁸ to deconvolution (sharpening) spectra that is simpler than the relatively complicated formalism developed by Burg.⁹ It is only necessary to replace $θ→θ+iα$ in Eq. (19), in which case, $Γ ν→ Γ ν−α$ for each $ν$⁠. With this substitution, every broadening parameter is reduced by $α$⁠. By increasing $α$⁠, we can observe the $Γ ν$ in sequence, “sweeping” through the singularities and identifying those that may otherwise be hidden.

The material for which XRD data were collected is a two-phase ZrO₂−33  mol. % YTaO₄ compound consisting of a monoclinic (I2/a, space group 15) YTaO₄-based phase and tetragonal (p42/nmc, space group 137) ZrO₂-based phase.¹⁷ This material is being investigated as a future thermal-protection material for components in the hottest regions of gas-turbine engines. It is synthesized from 325-mesh (<42 μm diameter) ZrO₂, Y₂O₃, and Ta₂O₅ oxide powders (Elemental Metals, Randolph, NJ) using a ball milling approach with steel grinding media followed by annealing at 1500 °C in air for 100 h. The powder XRD patterns were obtained using a Malvern Panalytical Empyrean x-ray diffractometer (Almelo, the Netherlands), set up with a Cu-K_α source, a Bragg-Brentano HD PreFIX module, combined with a reflection spinner stage, and a GaliPIX 3D with FASS detector in scanning line (1D) mode with an active length of 7.1763°. Data were acquired using a 1/8° divergence slit, a 1/2° antiscatter slit, a 10 mm mask, and a 0.04 rad Soller slit, scanned between diffraction angles 2θ of 15° –100° at a rate of 0.12°/s, resulting in a total scan time of 768.13 s, while samples were spun at 0.25 revolutions/s. No postprocessing was conducted after data collection, and thus, both the Cu-Kα₁ and Cu-Kα₂ x-ray emission peaks of the copper anode contributed to this unprocessed data.

As noted earlier, linear filters involve trade-offs among information loss, noise leakage, and Gibbs oscillations. To establish the transfer-function coefficients $H n$ relevant to the current data, we consider the largest window centered on the dominant peak that essentially does not overlap adjacent features. This is the 501-point window centered on 30.1074° and extending from 26.6074° to 33.6074° in Fig. 1. By Eq. (15), a larger N yields smaller uncertainties in the Fourier coefficients and, hence, gives better white-noise suppression. Figure 2 shows the resulting Fourier coefficients plotted as $ln⁡( C n)$⁠, along with the transfer function for order $M=4$ and a Gauss–Hermite cutoff $n G H=40$⁠. The presence of a second main peak in this cluster, which appears as a shoulder in Fig. 1, is evident in the slow oscillation seen as $ln⁡( C n)$ descends from its initial value. The approximately linear decrease in the average value of ln(C_n) with n shows that the lineshapes are essentially Lorentzian. For Gaussian lineshapes, ln(C_n) decreases quadratically with n. For Voight lineshapes, which are a convolution of Lorentzian and Gaussian functions, ln(C_n) decreases as the product of the Fourier transforms of the constituents, following the rule giving the Fourier coefficients of a convolution. A Gaussian cutoff is not seen in the data of Fig. 2, so the lineshapes here are essentially Lorentzian.

Figure 2 shows that the onset of white noise occurs at approximately $n=100$⁠. Accordingly, we design the transfer function to be flat for eight derivatives $(M=4)$ with a value of 0.5 at slightly less than $n=100$ $( n c=40)$⁠. This is somewhat conservative with respect to noise while yielding an effective balance among the conflicting apodization errors. Here, IRED preprocessing was done with a cutoff $n c=92$⁠.

The DS convolution coefficients $h j$ corresponding to $H n$ of Fig. 2 are shown on an expanded scale in Fig. 3. Only 27 of the 501 coefficients have values greater than 0.01% of the peak value and need to be retained. The useful coefficients take the form of a broadened delta function with overshoots at points $251±4$ of about 13% of the peak value. For sharper DS features, these phase reversals would generate Gibbs oscillations, but the effect here is essentially negligible. With over an order of magnitude fewer coefficients in this example, the sliding DS convolution is much more efficient than doing filtering in RS, even though the results are the same. Despite their “Mexican hat” appearance, these $h j$ are not a wavelet, since wavelets must integrate to zero.¹⁸ 

The effectiveness of this approach can be appreciated from Fig. 4, which shows segments of the original (black) data and that resulting from the 27-point GH convolution of the entire spectrum using the convolution coefficients shown in Fig. 3.²¹ Noise is eliminated on the scale of the figure, while the amplitudes of the extrema and their diffraction angles are retained.

The above parameters can be adjusted to optimize other criteria as desired. The most usual adjustment is the size of the window. If this is changed but the filtering is intended to remain the same, then the GH cutoff $n c$ would be scaled proportionally.

As our first example of nonlinear filtering, consider the features for scattering angles in a 101-point segment from 18. 4874° to 19.8874°, shown on an expanded scale in Fig. 5. Its small amplitude relative to the 30° maximum feature provides a better illustration of the capability of CME to eliminate noise and return data about the singularities present in a segment. The upper noisy (black) trace is the segment taken directly from Fig. 1. The segment contains three obvious structures, a dominant feature near 18.9° and two smaller features near 19.2° and 19.4°. The lower noisy (blue) trace shows the data after IRED preprocessing. The lower smooth (cyan) trace is the result of CME processing the IRED-preprocessed data. The upper smooth (red) trace completes the analysis, with the endpoint-discontinuity segment being reinserted, yielding the filtered version of the original. The result exhibits essentially no noise.

The processing parameters of Fig. 5 were obtained with the aid of Fig. 6. This figure shows $C n$ of the IRED-preprocessed data of Fig. 5 on a logarithmic scale. Both upper traces exhibit a white-noise level of $ln⁡( C n)∼(−6)$ with an onset at $n c∼8$⁠. Accordingly, $M=8$ was the order selected for CME processing. The IRED cutoff was selected to be 20 to remain well away from the information-containing range. The CME extrapolation is excellent, not only accurately describing the slope but also continuing the oscillations apparent in the nine available low-order coefficients, as expected from the theory. In the calculation, the dc value of the data was lowered by 0.2 so $ln⁡( C n)$ at $n=0$ is an accurate extrapolation of the other values.

Many approaches can reduce noise, but CME is unique in being able to return model-independent estimates of parameters describing singularities. Figure 7 shows the roots of the characteristic polynomial of order $M=8$ for the data of Fig. 5. The roots essentially divide into two classes: three “true” singularities characterized by widths $Γ ν$ of the order of $0.11 °$⁠, four “cosmetic” singularities with broadening parameters $Γ ν> 0.15 °$⁠, and one “intermediate” singularity with $Γ ν= 0.135 °$⁠. By inspection, the “intermediate” singularity can be attributed to the shoulder at approximately $18.8 °$⁠, although this feature is by no means obvious in Fig. 5. In this example, the “true” roots are separated by at least twice their broadening parameters and so are relatively independent. Consequently, their values can be expected to be relatively accurate representations of the parameters of the corresponding singularities of the response function describing diffraction.

As a method of assessing peaks, Fig. 8 shows the result of evaluating Eq. (19) with $α$⁠, defined in the discussion following Eq. (24), increasing from zero in 0.05° steps. The features turn sequentially into delta functions as $α$ increases, then fade out as $α$ increases further. The series highlights the shoulder at $18.78 °$⁠, which is not immediately obvious in Fig. 5. It can be noted that the signal-to-noise ratio of the two highest-angle peaks in this example is about 3, and that of the shoulder is about 1.

As a second example, Fig. 9 shows the main peak at 30.107° together with its shoulder on the higher-angle side. These data, consisting of 51 points from 29.757° to 30.457°, are shown on an expanded scale as the black trace in Fig. 9. The window width of 51 points was selected to center the structure and to avoid including the feature near 30.52°. As a nonlinear process, the CME functions most effectively with the fewest included structures.

The red trace is the IRED-preprocessed version. The Fourier coefficients $C n$ of this segment are shown on a logarithmic scale in Fig. 10. IRED preprocessing is seen to yield about half again as many information-containing coefficients for analysis. The break near $n c=12$ between the information-containing and noise-dominated coefficients is not as distinct as usual because the structure is a composite. The white-noise level of $ln⁡( C n)$ is approximately $(−5)$⁠, in contrast to the value $(−8)$ in Fig. 2, consistent with Eq. (3) given the factor-of-10 difference in the number of points in the window.

Figure 11 shows the locations of the roots of the characteristic polynomial in the complex plane. These are nominally the poles of the response function. The eight roots for $M=8$ clearly fall into two groups, two with noticeably smaller widths and six with larger widths scattered across the plot. Therefore, the diffraction maxima occur at 30.102° and 30.182° with half-width-half-maximum values of 0.060° and 0.062°, respectively. These values are consistent with the overall line shape, although only the width of the dominant feature can be assessed visually. It is not possible to estimate the amplitude of the feature, giving the shoulder structure by this procedure, but from Fig. 9, the amplitude is somewhat less than half that of the main peak.

An analysis of the peak at the side of the main structure yields a location of 30.536° and a broadening parameter $Γ= 0.063 o$⁠. The calculation was done for $M=1$⁠. Repeating the calculation for $M=2$ yielded a split peak, in agreement with the prediction of analytic calculations.²⁰ The split for $M=2$ occurs here because two Fourier coefficients are insufficient for accurate DS inversion.

The methods developed for UV/VIS spectroscopy are found to work well for diffraction data. The diffraction data consist of a relatively large number of values compared to UV/VIS applications, requiring somewhat different approaches to analysis. Aside from having to define the width of the scanning window, procedures follow those developed for use in UV/VIS spectroscopy. The goal of optimal enhancement of x-ray diffraction data by linear and nonlinear methods can be reached.

L.V.L. acknowledges the support of this work by the Vingroup Innovation Foundation (VinIF) (Grant No. VINIF.2022.STS.46) and the Vietnam Academy of Science and Technology (Grant No. THTEXS.01/21-24). Y.D.K. acknowledges the support of this work by the National Research Foundation of Korea (NRF) (Grant No. NRF-2020R1A2C1009041.) The x-ray diffraction instrument within the University of Virginia Nanoscale Materials Characterization Facility (NMcf) was fundamental to this work. J.A.D. and H.N.G.W. are grateful to the Office of Naval Research for support of this work through Grant No. N00014-21-1-2460 managed by David Shifler. The authors also acknowledge the anonymous reviewer who suggested that the shoulder in Fig. 9 may be the Cu- $K α 2$ replica of the main peak, which was found to be correct.

The authors have no conflicts to disclose.

Long V. Le: Writing – original draft (equal). Jeroen A. Deijkers: Data curation (equal). Young D. Kim: Writing – review & editing (equal). Haydn N. G. Wadley: Writing – review & editing (equal). David E. Aspnes: Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.