.. _xafs-ft_sec:

==============================================
XAFS: Fourier Transforms for XAFS
==============================================

.. module:: _xafs

Fourier transforms are central to understanding and using
XAFS. Consequently, many of the XAFS functions in Larch use XAFS Fourier
transforms as part of their processing, and many of the functions
parameters and arguments described here have names and meanings used
throughout the XAFS functionality of Larch.  For example, both
:func:`autobk` and :func:`feffit` rely on XAFS Fourier transforms, and use
the XAFS Fourier transform function described here.


Overview of XAFS Fourier transforms
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

The standard Fourier transform of a signal :math:`f(t)` can be written as

.. math::
   :nowrap:

   \begin{eqnarray*}
         {\tilde{f}}(\omega) &=& \frac{1}{\sqrt{2\pi}} \int_{-\infty}^{\infty}
       f(t) e^{-i{\omega}t} dt \\
       f(t) &=& \frac{1}{\sqrt{2\pi}} \int_{-\infty}^{\infty}
       {\tilde{f}}(\omega) e^{i{\omega}t} d{\omega} \\
   \end{eqnarray*}

where the symmetric normalization is one of the more common choices of
conventions.  This gives conjugate variables of :math:`\omega` and
:math:`t`. Because XAFS goes as

.. math::

  \chi(k) \sim \sin[2kR + \delta(k)]

the conjugate variables in XAFS are generally taken to be :math:`k` and
:math:`2R`.  The normalization of :math:`\tilde\chi(R)` from a Fourier
transform of :math:`\chi(k)` is a matter of convention, but we follow the
symmetric case above (with :math:`t` replaced by :math:`k` and
:math:`\omega` replaced by :math:`2R`, and of course :math:`f` by
:math:`\chi`).

But there are two more important issues to mention.  First, an XAFS Fourier
transform multiplies :math:`\chi(k)` by a power of :math:`k`, :math:`k^n`
and by a window function :math:`\Omega(k)` before doing the Fourier
transform.  The power-law weighting allow the oscillations in :math:`k` to
emphasize different portions of the spectra, or to give a uniform intensity
to the oscillations.  The window function acts to smooth the resulting
Fourier transform and remove ripple and ringing in it that would result
from a sudden truncation of :math:`\chi(k)` at the end of the data range.

The second important issue is that the continuous Fourier transform
described above is replaced by a discrete transform.  This better matches
the discrete sampling of energy and :math:`k` values of the data, and
allows Fast Fourier Transform techniques to be used.  It does change the
definitions of the transforms used somewhat. First, the :math:`\chi(k)`
data must be on *uniformly spaced* set of :math:`k` values.  The default
:math:`k` spacing used in Larch (including as output from :func:`autobk`)
is :math:`\delta k` = 0.05 :math:`\rm\AA^{-1}`.  Second, the array size for
:math:`\chi(k)` used in the Fourier transform should be a power of 2. The
default used in Larch is :math:`N_{\rm fft}` = 2048.   Together, these
allow :math:`\chi(k)` data to 102.4 :math:`\rm\AA^{-1}`.  Of course, real
data doesn't extend that far, so the array to be transformed is
*zero-padded* to the end of the range.  Since the spacing :math:`\delta R`
of the resulting discrete :math:`\chi(R)` is given as
:math:`\pi/{(N_{\rm fft} \delta k )}`, the extended range and zero-padding
will increase the density of points in :math:`\chi(R)`, smoothly
interpolating the values.   For :math:`N_{\rm fft}` = 2048 and
:math:`\delta k` =  0.05 :math:`\rm\AA^{-1}`, the :math:`R` spacing is
approximately :math:`\delta R` =  0.0307 :math:`\rm\AA`.

For the discrete Fourier transforms with samples of :math:`\chi(k)` at the
points :math:`k_n = n \, \delta k`, and samples of :math:`\chi(R)` at the
points :math:`R_m = m \, \delta R`, the definitions become:

.. math::
   :nowrap:

   \begin{eqnarray*}
   \tilde\chi(R_m) &=& \frac{i \delta k}{\sqrt{\pi N_{\rm fft}}} \,
   		       \sum_{n=1}^{N_{\rm fft}} \chi(k_n) \,
                       \Omega(k_n) \, k_n^w e^{2i\pi n m/N_{\rm fft}} \\
   \tilde\chi(k_n) &=& \frac{2 i \delta R}{\sqrt{\pi N_{\rm fft}}} \,
                       \sum_{m=1}^{N_{\rm fft}} \tilde\chi(R_m) \,
                       \Omega(R_m) \, e^{-2i\pi n m/N_{\rm fft}} \\
   \end{eqnarray*}


These normalizations preserve the symmetry properties of the Fourier
Transforms with conjugate variables :math:`k` and :math:`2R`.
Though the reverse transform converts the complex :math:`\chi(R)` to the
complex :math:`\chi(k)` on the same :math:`k` spacing as the starting
data, we often refer to the filtered :math:`\chi(k)` as *q* space.

A final complication in using Fourier transforms for XAFS is that the
measured :math:`\mu(E)` and :math:`\chi(k)` are a strictly real values,
while the Fourier transform inherently treats :math:`\chi(k)` and
:math:`\chi(R)` as complex values. This leads to an ambiguity about how to
construct the complex :math:`\tilde\chi(k)`.  In many formal treatments,
XAFS is written as the imaginary part of a complex function.  This might
lead one to assume that constructing :math:`\tilde\chi(k)` as :math:`0 +
i\chi(k)` would be the natural choice.  For historical reasons, Larch uses
the opposite convention, constructing :math:`\tilde\chi(k)` as
:math:`\chi(k) + 0i`.   As we'll see below, you can easily change this
convention.  The effect of this choice is minor unless one is
concerned about the differences of the real and imaginary parts of
:math:`\chi(R)` or one is intending to filter and back-transform the
:math:`\chi(R)` and compare the filtered and unfiltered data.


Forward XAFS Fourier transforms (:math:`k{\rightarrow}R`)
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

The forward Fourier transform converts :math:`\chi(k)` to :math:`\chi(R)`
and is of primary importance for XAFS analysis.  In Larch, this is
encapsulated in the :func:`xftf` function.

..  function:: xftf(k, chi=None, group=None, ...)

    perform a forward XAFS Fourier transform, from :math:`\chi(k)` to
    :math:`\chi(R)`, using common XAFS conventions.

    :param k:        1-d array of photo-electron wavenumber in :math:`\rm\AA^{-1}`
    :param chi:      1-d array of :math:`\chi`
    :param group:    output Group
    :param rmax_out: highest *R* for output data (10 :math:`\rm\AA`)
    :param kweight:  exponent for weighting spectra by :math:`k^{\rm kweight}`
    :param kmin:     starting *k* for FT Window
    :param kmax:     ending *k* for FT Window
    :param dk:       tapering parameter for FT Window
    :param dk2:      second tapering parameter for FT Window
    :param window:   name of window type
    :param nfft:     value to use for :math:`N_{\rm fft}` (2048).
    :param kstep:    value to use for :math:`\delta{k}` (0.05).

    :returns:  ``None`` -- outputs are written to supplied group.

    Follows the First Argument Group convention, using group members named ``k`` and ``chi``.  
    The following data is put into the output group:

       ================= ===============================================================
        array name         meaning
       ================= ===============================================================
        kwin               window :math:`\Omega(k)` (length of input chi(k)).
	r                  uniform array of :math:`R`, out to ``rmax_out``.
	chir               complex array of :math:`\tilde\chi(R)`.
	chir_mag           magnitude of :math:`\tilde\chi(R)`.
	chir_pha           phase of :math:`\tilde\chi(R)`.
	chir_re            real part of :math:`\tilde\chi(R)`.
	chir_im            imaginary part of :math:`\tilde\chi(R)`.
       ================= ===============================================================

    It is expected that the input ``k`` be a uniformly spaced array of
    values with spacing ``kstep``, starting a 0.  If it is not, the ``k``
    and ``chi`` data will be linearly interpolated onto the proper grid.

    The FT window parameters are explained in more detail in the discussion of
    :func:`ftwindow`.


..  function:: xftf_fast(chi, nfft=2048, kstep=0.05)

    perform a forward XAFS Fourier transform, from :math:`\chi(k)` to
    :math:`\chi(R)`, using common XAFS conventions.  This version demands
    ``chi`` to include any weighting and windowing, and so to represent
    :math:`\chi(k)k^w\Omega(k)` on a uniform :math:`k` grid. It returns
    the complex array of :math:`\chi(R)`.

    :param chi:      1-d array of :math:`\chi` to be transformed
    :param nfft:     value to use for :math:`N_{\rm fft}` (2048).
    :param kstep:    value to use for :math:`\delta{k}` (0.05).

    :returns:  complex :math:`\chi(R)`.


Reverse XAFS Fourier transforms (:math:`R{\rightarrow}q`)
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Reverse Fourier transforms convert :math:`\chi(R)` back to filtered
:math:`\chi(k)`.  We refer to the filtered :math:`k` space as :math:`q` to
emphasize the distinction between the two.  The filtered :math:`\chi(q)` is
complex.  By convention, the real part of :math:`\chi(q)` corresponds to
the explicitly real :math:`\chi(k)`.

..  function:: xftr(r, chir, group=None, ...)

    perform a reverse XAFS Fourier transform, from :math:`\chi(R)` to
    :math:`\chi(q)`.

    :param r:        1-d array of distance.
    :param chir:      1-d array of :math:`\chi(R)`
    :param group:    output Group
    :param qmax_out: highest *k* for output data (30 :math:`\rm\AA^{-1}`)
    :param rweight:  exponent for weighting spectra by :math:`r^{\rm rweight}` (0)
    :param rmin:     starting *R* for FT Window
    :param rmax:     ending *R* for FT Window
    :param dr:       tapering parameter for FT Window
    :param dr2:      second tapering parameter for FT Window
    :param window:   name of window type
    :param nfft:     value to use for :math:`N_{\rm fft}` (2048).
    :param kstep:    value to use for :math:`\delta{k}` (0.05).

    :returns:  ``None`` -- outputs are written to supplied group.

    Follows the First Argument Group convention, using group members named ``r`` and ``chir``.  
    The following data is put into the output group:

       ================= ===============================================================
        array name         meaning
       ================= ===============================================================
        rwin               window :math:`\Omega(R)` (length of input chi(R)).
	q                  uniform array of :math:`k`, out to ``qmax_out``.
	chiq               complex array of :math:`\tilde\chi(k)`.
	chiq_mag           magnitude of :math:`\tilde\chi(k)`.
	chiq_pha           phase of :math:`\tilde\chi(k)`.
	chiq_re            real part of :math:`\tilde\chi(k)`.
	chiq_im            imaginary part of :math:`\tilde\chi(k)`.
       ================= ===============================================================

    In analogy with :func:`xftf`, it is expected that the input ``r`` be a
    uniformly spaced array of values starting a 0.

    The input ``chir`` array can be either the complex :math:`\chi(R)` array
    as output to ``Group.chir`` from :func:`xftf`, or one of the real or
    imaginary parts of the :math:`\chi(R)` as output to ``Group.chir_re`` or
    ``Group.chir_im``.

    The FT window parameters are explained in more detail in the discussion of
    :func:`ftwindow`.


..  function:: xftr_fast(chir, nfft=2048, kstep=0.05)

    perform a reverse XAFS Fourier transform, from :math:`\chi(R)` to
    :math:`\chi(q)`, using common XAFS conventions.  This version demands
    ``chir`` be  the complex :math:`\chi(R)` as created from :func:`xftf`.
    It returns the complex array of :math:`\chi(q)` without putting any
    values into a group.

    :param chir:     1-d array of :math:`\chi(R)` to be transformed
    :param nfft:     value to use for :math:`N_{\rm fft}` (2048).
    :param kstep:    value to use for :math:`\delta{k}` (0.05).

    :returns:  complex :math:`\chi(q)`.


:func:`ftwindow`: Generating Fourier transform windows
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

As mentioned above, a Fourier transform window will smooth the resulting
Fourier transformed spectrum, removing ripple and ringing in it that would
result from a sudden truncation data at the end of it range.  There is an
extensive literature on such windows, and a lot of choices and parameters
available for constructing windows.  A sampling of windows is shown below.


..  function:: ftwindow(x, xmin=0, xmax=None, dk=1, ...)

    create a Fourier transform window array.

    :param x:        1-d array array to build window on.
    :param xmin:     starting *x* for FT Window
    :param xmax:     ending *x* for FT Window
    :param dx:       tapering parameter for FT Window
    :param dx2:      second tapering parameter for FT Window (=dx)
    :param window:   name of window type
    :returns:  1-d window array.

    Note that if ``dx`` is specified but ``dx2`` is not, ``dx2`` will
    generally take the same value as ``dx``.

    The window type must be one of those listed in the :ref:`Table of
    Fourier Transform Window Types <xafs-ftwin_table>`.

.. index:: Fourier Transform Window types
.. _xafs-ftwin_table:

    Table of Fourier Transform Window Types.

       =================== =========================================================
        window name          description
       =================== =========================================================
        hanning              cosine-squared taper
        parzen               linear taper
        welch                quadratic taper
	gaussian             Gaussian (normal) function window
	sine                 sine function window
	kaiser               Kaiser-Bessel function-derived window
       =================== =========================================================

In general, the window arrays have a value that gradually increases from 0
up to 1 at the low-k end, may stay with a value 1 over some central
portion, and then tapers down to 0 at the high-k end.  The meaning of the
``dx`` and ``dx2``, and even ``xmin``, and ``xmax`` varies a bit for the
different window types.  The Hanning, Parzen, and Welch windows share a
convention that the windows taper up from 0 to 1 between ``xmin-dx/2`` and
``xmin+dx/2``, and then taper down from 1 to 0 between ``xmax-dx2/2`` and
``xmax+dx2/2``.

The conventions for the Kaiser, Gaussian, and Sine window types is a bit
more complicated, and is best given explicitly.  In the formulas below,
``dx`` written as :math:`dx` and ``dx2`` as :math:`dx_2`.  We
define :math:`x_i = x_{\rm min} - dx/2`, :math:`x_f = x_{\rm max} +
dx_2/2`, and :math:`x_0 = (x_f + x_i)/2`, as the beginning, end, and
center of the widows.  For the Gaussian window, the form is:

.. math::

   \Omega(x) = \exp{\bigl[ -\frac{(x - x_0)^2}{2{dx}^2}\bigr]}

The form for the Sine window is

.. math::

   \Omega(x) = \sin{\bigl[ \frac{ \pi(x_f - x)}{x_f - x_i}\bigr]}

between :math:`x_i`  and  :math:`x_f`, and 0 outside this range.
The Kaiser-Bessel window is slightly more complicated:

.. math::
   :nowrap:

   \begin{eqnarray*}
      a       &=& \sqrt{\max{\bigl[0, 1 - \frac{4(x-x_0)^2}{(x_f-x_i)^2} \bigr] }} \\
     \Omega(x) &=& \frac{i_0(a dx) - 1}{i_0(dx) - 1} \\
   \end{eqnarray*}

where :math:`i_0` is the modified Bessel function of order 0.

Catalog of Fourier transform window
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Here, we give a series of example windows, to illustrate the different
window types and the effect of the various parameters.  The meanings of
``xmin``, ``xmax``, ``dx`` and ``dx2`` are identical for the Hanning, Parzen and
Welch windows, and illustrated in the two following figures.

.. subfigstart::

.. _fig-xafsft1a:

.. figure::  ../_images/ftwin_example1.png
    :target: ../_images/ftwin_example1.png
    :width: 100%

.. _fig-xafsft1b:

.. figure::  ../_images/ftwin_example2.png
    :target: ../_images/ftwin_example2.png
    :width: 100%

.. subfigend::
    :width: 0.45
    :label: fig-xafsft1

    Fourier Transform window examples and illustration of parameter meaning
    for the Hanning, Parzen, and Welch windows.  Note that
    :math:`\Omega(x=x_{\rm min}) = \Omega(x=x_{\rm max}) = 0.5`, and that
    the meaning of ``dx`` is to control the taper over which the window
    changes from 0 to 1.  Here, ``xmin=5`` and ``xmax=15``.



Some more window functions:

.. subfigstart::

.. _fig-xafsft2a:

.. figure::  ../_images/ftwin_example3.png
    :target: ../_images/ftwin_example3.png
    :width: 100%

.. _fig-xafsft2b:

.. figure::  ../_images/ftwin_example4.png
    :target: ../_images/ftwin_example4.png
    :width: 100%

.. subfigend::
    :width: 0.45
    :label: fig-xafsft2

    Fourier Transform window examples and illustration of parameter
    meaning.  On the left, a comparison of Welch, Parzen, and Hanning with
    the same parameters is shown.  On the right, the effect of ``dx2`` is
    shown as a different amount of taper on the high- and low-``x`` end of
    the window.  As before, ``xmin=5`` and ``xmax=15``.

The Gaussian, Sine, and Kaiser-Bessel windows are illustrated next. These
go to 1 at the average of ``xmin`` and ``xmax``, but do not stay at 1 over
a central portion of the window -- they taper continuously.  The Gaussian
window is a simple Gaussian function, and is not truncated according to
``xmin`` and ``xmax``, and the ``dx`` parameter sets the width.  The Sine
and Kaiser-Bessel windows both go to zero at  ``xmin-dx/2`` and ``xmax +
dx/2``.  For very large values of ``dx``, the Kaiser-Bessel window
approaches a nearly Gaussian lineshape.

.. subfigstart::

.. _fig-xafsft3a:

.. figure::  ../_images/ftwin_example5.png
    :target: ../_images/ftwin_example5.png
    :width: 100%

.. _fig-xafsft3b:

.. figure::  ../_images/ftwin_example6.png
    :target: ../_images/ftwin_example6.png
    :width: 100%

.. subfigend::
    :width: 0.45
    :label: fig-xafsft3

    Fourier Transform windows. On the left, a comparison of Kaiser-Bessel,
    Sine, and Gaussian windows with the same parameters is shown.  On the
    right, the effect of ``dx`` is shown for the Kaiser-Bessel window, and
    a closer comparison to a Gaussian window is made.


Examples: Forward XAFS Fourier transforms
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Now we show some example Fourier transforms, illustrating the real and
imaginary parts of the :math:`\chi(R)` as well as the magnitude, the effect
of different windows types, and Fourier filtering to :math:`\chi(q)`.  We
use a single XAFS dataset from FeO for all these examples, with a
well-separated first and second shell.  The full scripts to generate the
figures shown here are included in the *examples/xafs/* folder.


We start with a comparison of a small value of ``dk`` and a larger value.
A script that runs :func:`xftf`, changing on ``dk`` would look like::

    xftf(dat1.k, dat1.chi, kmin=3, kmax=13, dk=1, window='hanning',
         kweight=kweight, group=dat1)

    dat2 = group(k=dat1.k, chi=dat1.chi) # make a copy of the group
    xftf(dat2.k, dat2.chi, kmin=3, kmax=13, dk=5, window='hanning',
         kweight=kweight, group=dat2)

would result in the following results:

.. subfigstart::

.. _fig-xafsft4a:

.. figure::  ../_images/xft_example1.png
    :target: ../_images/xft_example1.png
    :width: 100%

.. _fig-xafsft4b:

.. figure::  ../_images/xft_example2.png
    :target: ../_images/xft_example2.png
    :width: 100%

.. subfigend::
    :width: 0.45
    :label: fig-xafsft4

    Comparison of the effect of different values of ``dk`` on real XAFS
    Fourier transforms.  Increasing ``dk`` reduces peak heights and tends
    to broaden peaks, but the effects are rather small.


A script that runs :func:`xftf` with consistent parameters, but different
window types::

    xftf(dat1.k, dat1.chi, kmin=3, kmax=13, dk=4, window='hanning',
         kweight=kweight, group=dat1)

    dat2 = group(k=dat1.k, chi=dat1.chi) # make a copy of the group
    xftf(dat2.k, dat2.chi, kmin=3, kmax=13, dk=4, window='parzen',
         kweight=kweight, group=dat2)

    dat3 = group(k=dat1.k, chi=dat1.chi) #
    xftf(dat3.k, dat3.chi, kmin=3, kmax=13, dk=4, window='welch',
         kweight=kweight, group=dat3)

    dat4 = group(k=dat1.k, chi=dat1.chi) #
    xftf(dat4.k, dat4.chi, kmin=3, kmax=13, dk=4, window='kaiser',
         kweight=kweight, group=dat4)

    dat5 = group(k=dat1.k, chi=dat1.chi) #
    xftf(dat5.k, dat5.chi, kmin=3, kmax=13, dk=4, window='gaussian',
         kweight=kweight, group=dat5)

would result in the following results:

.. subfigstart::

.. _fig-xafsft_winsa:

.. figure::  ../_images/xft_example3.png
    :target: ../_images/xft_example3.png
    :width: 100%

.. _fig-xafsft_winsb:

.. figure::  ../_images/xft_example4.png
    :target: ../_images/xft_example4.psng
    :width: 100%

.. subfigend::
    :width: 0.45
    :label: fig-xafsft_wins

    Comparison of the effect of different window types on real XAFS Fourier
    transforms.

We now turn our attention to the different components of the Fourier
transform.  As above, it is most common to plot the magnitude of the
Fourier transform.  But, as the transformed :math:`\chi(R)` is complex, it
can be instructive to plot the real and imaginary components, as shown
below::

    newplot(dat1.r, dat1.chir_mag, xmax=8, label='chir_mag',
            show_legend=True, legend_loc='ur', color='black',
            xlabel=r'$R \rm\, (\AA)$', ylabel=r'$\chi(R)\rm\,(\AA^{-3})$' )

    plot(dat1.r, dat1.chir_re, color='red', label='chir_re')
    plot(dat1.r, dat1.chir_im, color='blue', label='chir_im')

which results in

.. _fig-xafsft_ri:

.. figure:: ../_images/xft_example5.png
   :target: ../_images/xft_example5.png
   :width: 65%
   :align: center

   The real and imaginary components of the XAFS Fourier transform.

In fact, in the analysis discussed with :func:`feffit`, the real and
imaginary components are used, not simply the magnitude.

Examples: Reverse XAFS Fourier transforms, Fourier Filtering
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

A reverse Fourier transform will convert data from :math:`\chi(R)` to
:math:`\chi(q)`.  This allows a limited range of frequencies (distances) to
be isolated and turned back into a :math:`\chi(k)` spectrum.  Here, we show
two different :math:`R` windows to filter either just the first shell of
the spectra, or the first two shells, and compare the resulting filtered
:math:`\chi(q)`.


.. subfigstart::

.. _fig-xafsft_filtera:

.. figure::  ../_images/xft_example6.png
    :target: ../_images/xft_example6.png
    :width: 100%

.. _fig-xafsft_filterb:

.. figure::  ../_images/xft_example7.png
    :target: ../_images/xft_example7.png
    :width: 100%


.. subfigend::
    :width: 0.45
    :label: fig_xafsft_filter

    Reverse XAFS Fourier transform, or Fourier filtering.  Here, one can
    see the effect of different window sizes on the Fourier filtered
    spectrum.  Including the first two peaks or shells reproduces most of
    the original spectrum, with only high-frequency components removed.

Note that it is ``chiq_re`` that is compared to the k-weighted ``chi``
array.


