Question about Normalization of Foward Fast Fourier Transform

[Ameyanagi]

@newville I have a question about the Normalization of FFT in larch. I was reading the Documentation and saw that normalization is treated as symmetrical, which means it is divided by the 1/sqrt(N).

However, the current implementation seems to be using 1 for forward FFT. (at least for scipy.fftpack.fft)

Is this intentional behavior? I was confused about what is the standard for normalizing the amplitudes in EXAFS analysis.

[newville]

@Ameyanagi Good question, we should check this carefully. In fact, we expect to normally use

from mkl_fft.interfaces.numpy_fft import fft, ifft

and then use these (at https://github.com/xraypy/xraylarch/blob/65b23e7d028d1a6c2315535c32c1bf0f53bb406e/larch/xafs/xafsft.py#L321 and https://github.com/xraypy/xraylarch/blob/65b23e7d028d1a6c2315535c32c1bf0f53bb406e/larch/xafs/xafsft.py#L346). And there we do (or intend to do) the proper scaling for symmetric FFTs.

We're using mkl_fft.interfaces.numpy_fft because it seems to be the fastest. But that might have different scaling from scipy.fftpack (not sure)? It looks like we should probably fall back on the more equivalent numpy.fft routine, and also check that these results are actually symmetric. I have not looked at the details, yet, but if you have any insight please let me know.

[Ameyanagi]

@newville I just had a quick look into the source code. For mkl_fft, it uses backward normalization, which is 1 for forward FFT and 1/N for inverse FFT. https://github.com/IntelPython/mkl_fft/blob/master/mkl_fft/_numpy_fft.py

def fft(a, n=None, axis=-1, norm=None):
    """
    Compute the one-dimensional discrete Fourier Transform.

    This function computes the one-dimensional *n*-point discrete Fourier
    Transform (DFT) with the efficient Fast Fourier Transform (FFT)
    algorithm [CT].

    Parameters
    ----------
    a : array_like
        Input array, can be complex.
    n : int, optional
        Length of the transformed axis of the output.
        If `n` is smaller than the length of the input, the input is cropped.
        If it is larger, the input is padded with zeros.  If `n` is not given,
        the length of the input along the axis specified by `axis` is used.
    axis : int, optional
        Axis over which to compute the FFT.  If not given, the last axis is
        used.
    norm : {"backward", "ortho", "forward"}, optional
        .. versionadded:: 1.10.0

        Normalization mode (see `numpy.fft`). Default is "backward".
        Indicates which direction of the forward/backward pair of transforms
        is scaled and with what normalization factor.

        .. versionadded:: 1.20.0

            The "backward", "forward" values were added.

numpy fft is also the same. https://numpy.org/doc/stable/reference/generated/numpy.fft.fft.html#numpy.fft.fft

[newville]

@Ameyanagi Thanks -- that suggests "no normalization", and then we do the normalization in xafsft. But, we should check that too ;)

[Ameyanagi]

@newville I think we should discuss this very carefully. Because I don't see much difference in R space with the Demeter package, indicating that the implementation of the XAFSFT is the same for both of the programs. If we change one of it, it might affect all of the fitting processes.

By the way,I currently have a m1 Macbook, in which, it falls back to scipy. Do you think this is the reason why the below code does not match? I will look into it, but I will raise up as an issue. The scipy.fftpack seems not to normalize at all. (1 for forward FFT and 1 for inverse FFT.) https://docs.scipy.org/doc/scipy/reference/generated/scipy.fftpack.ifft.html#scipy.fftpack.ifft

from larch import Group 
from larch.xafs import pre_edge, autobk, xftf, xftr
import numpy as np

data = np.loadtxt("tests/testfiles/Ru_QAS.dat")

mu = np.log(data[:, 1] / data[:, 2])
energy = data[:, 0]

group = Group(energy=energy, mu=mu, filename='Ru_QAS.dat')

autobk(group)
xftf(group)
xftr(group)

import matplotlib.pyplot as plt

plt.plot(group.k, group.chi*group.k**2, label='chi')
plt.plot(group.q, group.chiq_re, label='chi_re')
plt.legend()
plt.xlabel('k')
plt.ylabel('chi(k) * k^2')

[Ameyanagi]

@newville The issue above was not related to the FFT, it was because I used default window functions. Sorry for the confusion. xftr(group, dr=0, window='hanning') gave a good match

[newville]

@Ameyanagi Ok, that makes sense. Still, I think it is worth double-checking that we're doing the scaling correctly (or at least consistently). I haven't had a chance to dig into that, but I'll try to get to it soon.

[Ameyanagi]

@newville I summarized the issue of the normalization as the following for clarity.

• Applying Discrete Fourier Transform (DFT) and Inverse Discrete Fourier Transformation (IDFT) will give a scale of N. Which therefore needs a normalization factor.
• The documentation says it is scaled by 1/sqrt(N) for DFT and 1/sqrt(N) for IDFT. But the actual implementation of (mkl_fft.interfaces.numpy_fft) is 1 for DFT and 1/N for IDFT. (This is the reason the q space is correct)
• The fallback of ffts and iffts needs to be double-checked that it is scaling correctly or consistently.

After double-checking the functions, I would suggest changing descriptions of the normalization factors in the documentation.

Do you know any literature in EXAFS that mentions the normalization factors of FFT? I know that on p.200 of "Introduction to XAFS" by Grant Bunker, it is stated that the DFT is usually defined as "1/sqrt(N) * ...", but I don't if many of the software follows this normalization factor. (IFEFFIT, REX2000, XTUNES, and more?)

It might be worth, summarizing how the R space of each software scale within each other, so that we don't get confused about the plots obtained by different software.