Chapter 5.3. Fourier transforms in EXAFS

In their seminal series of papers, Stern, Sayers and Lytle first showed that Fourier transformation (FT) of extended X-ray absorption fine-structure (EXAFS) spectra can be a useful tool in data analysis, both as a graphical representation of the distribution of atoms and as an intermediate step in data processing (Sayers et al., 1971; Stern, 1974; Lytle et al., 1975; Stern et al., 1975). The utility of FTs is well established, but some understanding is necessary if they are to be used correctly. Some useful review articles and books include Lee et al. (1981), Stern & Heald (1983), Koningsberger & Prins (1988) and Bunker (2010).

The applicability of FT is naturally suggested by the mathematical form of the single-scattering EXAFS equation as a sum of damped sine waves. The atoms that are farther away from the central absorbing atom give rise to more rapid oscillations in the EXAFS signal when plotted versus the electron wavenumber k. Within some limitations, FT of the EXAFS spectra allows the contributions from atoms at different distances to be separated and characterized.

In its simplest approximation, the EXAFS is a simple sum over the contributions from atoms from which the outgoing photoelectron wave is scattered. In the single-scattering formulation, the sum over atoms is recast as a sum over groups (coordination `shells') of atoms, which are at a similar distance and are taken to have the same atomic number; atoms with different atomic numbers are conceptually assigned to distinct shells (or subshells).

The variation in distances to the atoms within a shell is accounted for by the use of an EXAFS Debye–Waller factor (DWF). If the distribution of atoms is broad and highly skewed, it may also be necessary to include additional terms in the phase shift. The DWF approximation is well justified if the higher order (n > 2) cumulants C_n of the distributions of atoms in a shell are sufficiently small compared with , where k_max is the upper limit of the data-analysis range. If they are not, the atoms must be reassigned to additional subshells or higher order terms in the cumulant expansion should be included. This criterion directly follows from the relative importance of terms in the expansion = (Bunker, 1983, 2010) and the fact that the higher order cumulants have the greatest weight at k_max. Here, A_s(k) and φ_s(k) are the structural amplitude and phase, i.e. the amplitude and phase with the scattering amplitude divided out and the scattering phase subtracted out.

In the more general multiple-scattering theory (Ashley & Doniach, 1975; Lee & Pendry, 1975; Rehr & Albers, 2000), the sum over atoms becomes a sum over scattering paths; nevertheless, the basic form of the EXAFS equation as a sum of damped sine waves is preserved.

For simplicity, we use the conventional EXAFS expression for χ = (μ − μ₀)/μ₀ as a sum over coordination shells, where the amplitude A_j of the jth shell is given by and the phase is given by Here, N_j, R_j and are the average coordination number, the average distance and the mean-square variation in distance, respectively, from the central atom to the atoms in shell j, which as described above are taken to be of the same atomic number, |F_j(k)| is the scattering amplitude [the modulus of the complex scattering factor F_j(k)] and is its phase, δ_l(k) is the phase shift of the electron of final-state angular momentum l owing to the attractive central atom potential, which contributes twice, and and the electron mean free path λ(k) account for multi-electron interactions. The mean free path limits the effective range of the outgoing photoelectron and is an essential ingredient for understanding EXAFS.

The form of the EXAFS equation suggests that the spectral contribution from a specific coordination shell will be localized in peaks at specific distances in the Fourier transform. The peaks are systematically shifted to shorter apparent distances by about 0.4 Å owing to the negative slope of the phase-shift functions, especially the term 2δ_l(k) associated with the central atom potential, which is the same for all shells and the slope of which is similar for most central absorbing atoms.

In the next few sections we will discuss FTs in general and then later specialize to EXAFS. The Fourier transform (FT) of a function f(t) (here considered as a signal) is an operation that, given f(t), yields another function that represents a decomposition of the signal f(t) into a sum of distinct sinusoidal components dependent on a variable angular frequency ω. The FT can be considered to be a continuous-frequency analogue of a list of Fourier series coefficients, which represent the magnitudes of discrete frequency components in a signal. The resulting complex function has a modulus (amplitude) that shows peaks with heights and locations that indicate the strengths of the various different frequency components of which the signal f(t) is composed. Since is a complex function it carries both amplitude (modulus) and phase information.

The FT of a function f(t) is defined as where i = (−1)^1/2 and, by Euler's relation, exp(iωt) = cosωt + isinωt. The FT preserves the information content of f(t), and the original signal can be recovered exactly by performing an inverse FT (here abbreviated `IFT'). This is performed simply by changing i → −i in the exponential, or alternatively by complex conjugating, forward transforming and complex conjugating again. Various conventions exist as to which sign for i is considered to be the forward direction and how to distribute the prefactors in front of the integrals, the necessary requirement being that a cumulative prefactor of 1/(2π) must be included for a full transform/inverse transform cycle.

In actual data analysis signals only extend over a finite range and the integral is accordingly performed over a finite range. This is equivalent to multiplying the function f(t) by a window function that extends over a finite range.

Consider the integral , where z = ωτ/2. The sin(z)/z function has a maximum value of 1 at z = 0, with decreasing, oscillating tails; the width of its central maximum as a function of ω is inversely proportional to τ (see Fig. 1). The total area of the right side (in ω space) is 2π irrespective of τ. As τ → ∞ (as in the limits of the FT above) the height of the peak becomes infinite and the width goes to zero, but the area stays constant. These properties (with the area normalized to 1) are characteristic of a Dirac delta `function' δ(ω) (which is actually not a function, but a generalized function or `distribution'). Its characteristic feature is that for any variable x, , i.e. integrating the product of any function f times a delta function at some location x₀ plucks out the value of the function there, f(x₀). Evidently δ(x) = δ(−x). We conclude that = = . This expresses the property of orthogonality between functions of the form that are of different frequencies.

Figure 1 Sin(z)/z for z = ωτ/2 and τ = 1, 2, 4, 8.

The concept of convolution is important for understanding FTs in EXAFS analysis. The convolution of two functions f(t) and g(t) is a binary operation (represented by ) that produces a new function by planting a shifted replica of the function f at every point, weighted by the value of the function g at that point, and summing (integrating) them; the result is a function of the shift variable. Changing the dummy integration variable in the integral above, t′ → t − t′′, shows that the integral is symmetrical between f and g, so = . Convolution is an essential operation for describing spectral broadening, probabilities, diffraction and many other concepts.

The convolution theorem (which can be proved using the considerations above) says that the FT of a convolution of two functions is the product of the FT of each of the separate functions: This confirms that indeed = . Another corollary is that the FT of a product of two functions is the convolution of the FTs of the separate functions: Similar equations apply for inverse FTs.

In passing, we point out that an efficient way to perform convolutions using fast Fourier transforms (see below) is to separately transform the two functions, multiply the transforms together and then inverse-transform.

The EXAFS data χ(k) are functions of the wavenumber k, which is the analogue of t above. The analogue of ω is 2R, the effective length of the scattering path, i.e. twice the inter­atomic distance for single scattering, plus the (negative) slope of the other phase-shift terms. The data range for χ(k) may be limited at the low end (typically to k ≥ 2–3 Å⁻¹) because of the complications (such as high-order multiple scattering) that may occur in the near-edge region and that are not modelled by the EXAFS equation. They are restricted on the high end (typically k ≤ 12–16 Å⁻¹) by a limited signal-to-noise ratio, because the EXAFS signals decay owing to the XAFS Debye–Waller factor and the intrinsic decrease of atomic scattering amplitudes at high k, which is more pronounced for low atomic number elements. The most favourable case is a highly regular structure containing high atomic number atoms measured at temperatures of only few kelvin.

Fourier transforming a sine wave over an infinite range gives two infinitely sharp peaks [Dirac delta functions δ(ω − ω₀)] of opposite sign, one corresponding to the signal frequency and the other corresponding to the negative of that frequency. This is because (by Euler's relation) sin(ω₀t) = [exp(iω₀t) − exp(−iω₀t)]/(2i) and FT[exp(−iω₀t)] = (2π)^1/2δ(ω − ω₀).

Limiting the range of the FT to a finite interval between k₁ and k₂ is equivalent to multiplying the data by a rectangular box function (`window function') of unit height that extends from k₁ to k₂. The FT of a product of two functions (for example a window function and a sine wave) is equivalent to convoluting the FTs of the two separate functions. Convoluting with the FT of the window function broadens and distorts the ideal FT of the data alone. This is an inevitable fact that the data analyst must cope with.

The FT of a rectangular window function extending from k₁ to k₂ is where and Δk = k₂ − k₁. The exponential term simply gives an overall phase factor (unit modulus), which is not interesting, but the sin(Δkr)/(Δkr) dependence of the window FT produces a central peak with oscillating tails. If one Fourier transforms an ideal sine wave over a finite range in k-space the result is not perfectly sharp peaks, but broader peaks with a width that is inversely proportional to the FT k-range and with oscillating tails, which appear as side lobes when the modulus of the FT is plotted in r-space. This phenomenon is sometimes described as `ringing' or `spectral leakage'. This effect can be reduced by using a window function that rises and falls less sharply than a rectangular box. The goal is to minimally alter the k-space data over most of its range, while tapering the ends. A veritable zoo of such window or `apodization' functions (for example Hanning, Bessel) exist, each with their own advantages and disadvantages and their own adherents and detractors. The main goal is to limit the spread of the peaks in r-space while preserving the information content of the data, a procedure that may require some experimentation for different data sets. The key point is to be consistent in the choice of windows so the data and reference (or theoretical) spectra are processed in the same way.

The amplitudes of experimental EXAFS data decrease fairly rapidly as a function of k because of the intrinsic behaviour of scattering amplitudes and also the EXAFS Debye–Waller factors. The FT of a decaying sine wave (in k-space) is a broadened and asymmetric peak; the decay in k-space can be viewed as being caused by different r-space components destructively interfering with each other. We wish to minimize the width and tails of the r-space peaks. A simple way to do this is to make the k-space data appear approximately like a uniform amplitude sine wave, which is normally accomplished by multiplying χ(k) by a positive power of k^w, where typically w = 1–3. This, in combination with a tapered window function, gives reasonably well localized peaks in r-space. Choosing a k-space range, k-weight and window function is generally an act of compromise because the optimization that gives the best localization for one shell may be different to that for another shell because of the differing k-dependencies of their corresponding χ(k) signals.

These concerns are relevant whether the data are to be fitted in r-space or inverse-transformed into k-space (see below). The key goal is to minimize spectral leakage between peaks in r-space. Even if one limits the fitting range to a certain region in r-space, if the contributions from other peaks are leaking into that region, and such contributions are not considered to be part of the fitting model, systematic errors will occur.

The XAFS oscillations can be expressed as The two terms correspond to positive-frequency and negative-frequency components, which produce distinct r-space peaks. The complex function A(k)exp[iφ(k)] can be determined by zeroing out everything except the desired positive-frequency components, inverse-transforming to k-space and multiplying by 2i. This is possible if the positive- and negative-frequency components (and the contributions from other shells) are well separated from each other, which implies that their separation must be much larger than the broadening due to the k-window. This combination of criteria is equivalent to requiring that a sufficient number of oscillations of the lowest frequency component be included within the k-space transform window.

The resulting complex function Aexp(iφ) has real (α) and imaginary (β) parts, , from which the amplitude and phase can be determined as and . The phase constructed in this way is only determined modulo 2π and will have discontinuities, but these can be automatically removed in a postprocessing step to produce a continuous φ(k) which can then be analysed and fitted [along with A(k)] as desired in k-space.

In this way, the single-shell contributions can be isolated from each other provided that the transform peaks of the various shells are well isolated in r-space. This is not always the case. Firstly, the distances between the shells of atoms need to be large enough that their signals destructively interfere within the data range. This occurs when the signals `beat' with each other. For a beat to occur the difference in phase must be an odd multiple of π, giving a minimum condition 2k_beatΔR > π. To even begin to resolve the peaks the beat should occur in the middle of the k-space window. To fully resolve the peaks the difference in distance should be greater yet. This shows the advantages of a high k-window cutoff, but other practical considerations limit this.

Often adjacent peaks that are actually strongly interfering with each other will have destructive interference between the peaks (because of the nonlinearity of the modulus operation), and because of the small amplitude there this could naively be viewed as a good location to choose an r-space window boundary. Doing so can result in serious systematic errors whether one is fitting the data in k-space or in r-space.

Excluding the contributions from negative frequencies as well as the positive-frequency contributions from other shells requires a suitable choice of r-space window. This involves similar considerations as in the choice of a k-space window. We know that inverse-transforming the product of the r-space window (which could be rectangular or tapered) will result in the inverse-transformed data being convoluted with the inverse FT of the r-space window. This implies that the inverse-transformed data will be affected near the k-space window edges over a width that is inversely proportional to the r-space window. Typically, this results in the inverse-transformed data being distorted within about ∼1–2 Å⁻¹ of the k-space window boundary. Fortunately, this effect can be made to cancel out to a good approximation or deconvoluted (Bunker, 2010).

Fourier transforms of EXAFS data are not normally calculated as integrals over continuous variables as described in the previous sections; instead, a discrete Fourier transform (DFT) is performed, normally via a computationally efficient class of algorithms called fast Fourier transforms (FFTs).

To accomplish this, the sampled EXAFS data in energy space are interpolated onto a uniform grid in k-space that is fine enough to adequately sample the highest frequency in the data. The minimum sampling frequency to use is at least twice that of the highest frequency oscillation that might be contained in the data. This minimum sampling frequency (i.e. sampling twice per oscillation of the fastest component) is called the Nyquist frequency. This criterion is explained further below. Because of the limitations of sampling, any high-frequency signal that might be present in the data and that is higher than the Nyquist frequency will appear to be of (`alias to') a lower frequency. If the high-frequency components are due to noise, they too will alias down into the pass-band of the data and increase the errors. Smoothing with an analog filter before sampling reduces this problem, but this is not always practical to do. For this reason, it is good practice to oversample the signal at a rate several times greater than the Nyquist frequency in order to prevent noise or spurious high-frequency signals aliasing into the pass-band.

The DFT may be defined as The orthogonality relation for this discrete case (which corresponds to the continuous case considered earlier) is , where the Kronecker delta δ_lm is 1 for m = l and 0 otherwise, where l and m are integers. This implies that performing DFT on a pure signal of frequency l (in units of a base frequency), i.e. signal exp(−2πinl/N), will produce a single entry of magnitude N in the DFT array position corresponding to this frequency, with zeroes in all other locations. However, this is only true for signal frequencies between zero and frequency N/2 (in units of the base frequency; here, we assume N is even for simplicity), because the entries for N/2 ≤ m ≤ N correspond to negative-frequency components that essentially mirror the elements for 0 ≤ m ≤ N/2 and they provide no additional information for a real-valued input signal [the basis functions exp(2πinm/N) for the second half of the array are just complex conjugates of the first half, so the highest frequency signal that is resolvable occurs at N/2]. Another way to say this is that the highest frequency component of a signal that is resolvable by a DFT is sampled twice in a cycle, which is Nyquist's sampling theorem. Higher frequency signals present in the data alias to lower frequencies.

Experimental XAFS data normally come in digital form, sampled on some grid in energy space, which is then converted to k-space, where k = [0.26247(E − E₀)]^1/2 in eV Å units. The EXAFS data are interpolated onto a uniform grid extending from k_min to k_max with grid spacing δk: χ(k_n) = χ[k_min + (n − 1)δk] = χ_n and δk = (k_max − k_min)/(N_pts − 1). We assume that the digital values have sufficient bit depth (typically 32 bits) that they can be approximated as continuous variables.

For EXAFS, the frequency that we are referring to is the rapidity of the oscillations in k-space, which is related to the interatomic distances: the more distant the coordination shells of atoms, the more rapidly the EXAFS signal oscillates in k-space. The electron mean free path and other effects normally restrict the relevant distances to less than ∼5 Å. Considering Nyquist frequency issues, and an appropriate level of oversampling, in practice a k-space grid of δk = 0.05 Å⁻¹ has generally proved sufficient.

Comparing the forms of the continuous and discrete transforms tells us that the grid spacings δk and δr are related to the number of points in the sum N by Nδkδr = π. If we were to take N to be the number of k-space data points between k₁ and k₂, then δr = π/(k₂ − k₁ + δk), which numerically is much too coarse to adequately sample the peaks in r-space. For example, if k₁ = 2 Å⁻¹, k₂ = 12 Å⁻¹ and δk = 0.05 Å⁻¹ there are 201 points and the r-space grid would be 0.31 Å. The solution is to artificially extend the array size to a larger number (effectively increasing the resolution) and fill in the non­existent high-k data with zeroes (`zero-padding'). Interestingly, this is equivalent to interpolating between the data points that one would obtain using the coarser mesh, and it agrees precisely with the results one obtains from sampling the continuum integral (including the window function). In practice using N = 2048 has worked well (this is a power of two, which some FFT algorithms require) and gives δr ≃ 0.03Å for δk = 0.05 Å⁻¹.

Another aspect of DFTs (and FFTs) is that they are periodic: if you change n → n + N or m → m + N you obtain the same value. This can be pictured by thinking of the array of N points as a ring of array locations. In this case the positive R values appear in the low-index part of the array and the negative R values are mirrored near the end of the array. The maximum value of r that is representable is therefore (N/2)δr = π/δk.

Suppose one has N_true ≃ Δk/δk independent data points in k-space, where Δk is the width of the range in k-space. If one then transforms the data to r-space, and then inverse-transforms the data over the entire r-space range, the original function is recovered exactly. If one instead selects a slice of width Δr (by zeroing the contributions outside this range) and then inverse-transforms, the total possible information content of these `Fourier-filtered' data, which in principle could be spread out over the complete r-space range Nδr/2, is reduced by the fraction of the total r-space that is contained in the pass-band. Therefore, the fraction of the total information that is contained in the selected r-space slice is F = (ΔR/δr)(N/2), which corresponds to a total number of independent data points N_true(ΔR/δr)(N/2) = (Δk/δk)(Δr/δr)(N/2) = 2ΔkΔr/π. This is a useful expression that provides an approximate upper limit on the number of fitting parameters when fitting data that are only appreciable over limited data ranges of width Δk and Δr.

As stated above, the FT of the EXAFS [here called ] is a complex-valued function with a real and an imaginary part. The FT itself is a linear operation: the FT of a sum of signals (here called `shells') is the sum of the FTs of these signals. If one extracts the real or imaginary parts of the complex FT, then linearity (superposition) still holds for various shells, because taking the real or imaginary part is also a linear operation. However, the most common representation of the FT is the modulus , which is the square root of the sums of the squares of the real and imaginary parts. The modulus is not a linear function of , and superposition does not hold for those modulus plots. This makes the interpretation of the modulus plots a little more subtle because of interference.

As described in previous sections, the observed FT is the FT of the raw k-weighted data convoluted with the FT of the k-space window function. Fitting in r-space does not change this fact, and concerns about spectral leakage also apply here. If the data are sharply truncated in k-space, the `ringing' will introduce spurious side-lobes that definitely cannot be interpreted as real shells or peaks. Tapering the edges of the k-window reduces this effect, but not entirely. When comparing different data sets, it is essential to use the same transform ranges and k-weightings. Well designed analysis programs transform the theoretically modelled data in the same manner as the experimental data so as to compensate for these effects. The width of FT peaks is controlled to a large degree by the FT range, k-weighting and k-window shape. The width cannot easily be used as a measure of the disorder of atoms in a shell.

It is important to understand what peaks in the modulus of a FT represent. Essentially, the position of a peak represents the average slope of the phase of χ(k) for that shell over the transform range. This can be affected by any number of things. It is primarily affected by the average distance to atoms in the shell, but is also affected to a lesser extent by the transform range and k-weighting, because phases generally have some curvature as a function of k. Changing the energy reference E₀ by a few electronvolts effectively adds a term proportional to 1/k to the phase, which shifts the average phase slope. Highly skewed disordered distributions of atoms in the shell add a contribution (third cumulant) to the phase that is proportional to k³, which shifts the average phase slope and peak positions, sometimes to a very significant degree. Interference from other shells can apparently shift peaks. In fact, an asymmetric distribution of atoms that is skewed positively (to higher r) would appear to be a shorter distance. Therefore, it is unreliable to directly determine distances from FT peak positions with any degree of accuracy. Detailed analysis should be performed through fitting in k-space or r-space.

For atoms with sufficiently high atomic number (for example silver), scattering amplitudes and phases can exhibit significant structure as a function of k over the transform range, owing to Ramsauer–Townsend-type resonances in the backscattering atom potential (Lee et al., 1981), and this amplitude and phase modulation can strongly affect the peak shapes. For example, if the scattering amplitude has a strong minimum in the middle of the transform range, the FT peak can split into two peaks that naively might be interpreted as two shells of atoms. It is possible to compensate for this amplitude and phase modulation before performing the FT (this variant is called an `optical transform') if only one backscattering element is involved. If this is performed correctly, the amplitude and imaginary part of the FT will peak at the same location.

The FT peak height of a shell is a measure of the average amplitude of the signal of this shell over the transform range in k-space. This is affected by the coordination number: if the coordination number is doubled, and everything else is kept the same, then the peak height will double. However, the peak heights are also strongly affected by the transform window and the k-weighting used, and the Debye–Waller factor (which is typically different for various shells and is temperature dependent). If there is compositional disorder in a shell (small variations in atomic number), differences in the scattering phase shifts of the various atoms of different atomic number can cause some reduction in amplitudes and an apparent reduction in coordination number.

If one compares a material at a series of temperatures, different shells will normally have different variations in the Debye–Waller factor with temperature. This differential temperature dependence of peak heights versus temperature can be used to distinguish atoms that are tightly bound from those that are more weakly bound to the central atom.

The modulus alone tells one nothing about the phase of the signal: if you were to change the phase by π radians then χ(k) and would invert, but the modulus would be unchanged. It is often useful to plot either the real or imaginary part of the signal along with the modulus. Comparing the shapes of the imaginary (or real) parts with respect to the modulus gives an estimate of the constant part of the scattering phase, which depends on the atomic number. These plots can be useful for approximately identifying the type of atoms that are present in a shell. Such a plot is shown in Fig. 2. Detailed analysis, however, should be performed by curve-fitting or other methods.

Figure 2 FT modulus (solid line) and imaginary part (dashed line) of k3-weighted MnO data measured at 80 K, transformed from k = 2 to 12 Å−1.

In summary, despite many advances in other analysis methods, FT methods continue to find broad application in EXAFS data analysis. They must, however, be used with care and an understanding of their limitations.