Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Question about Normalization of Foward Fast Fourier Transform #465

Open
Ameyanagi opened this issue Sep 27, 2023 · 7 comments
Open

Question about Normalization of Foward Fast Fourier Transform #465

Ameyanagi opened this issue Sep 27, 2023 · 7 comments

Comments

@Ameyanagi
Copy link
Contributor

Ameyanagi commented Sep 27, 2023

@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).

image

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
Copy link
Member

@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

def xftf_fast(chi, nfft=2048, kstep=0.05, _larch=None, **kws):

and
def xftr_fast(chir, nfft=2048, kstep=0.05, _larch=None, **kws):
). 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
Copy link
Contributor Author

@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
Copy link
Member

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

@Ameyanagi
Copy link
Contributor Author

@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')

image

@Ameyanagi
Copy link
Contributor Author

@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
Copy link
Member

@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
Copy link
Contributor Author

@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.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants