Python Calculates Fourier Transforms

Technical background

Fourier transforms can be used in almost any computationally relevant field, such as filtering in communications, lattice inverse space calculations in materials, and inverse force field energy terms in molecular dynamics, to name a few. The simplest example is the calculation of the potential energy of a periodic box$k\sum_i\frac{q_i}{r_i}$In itself, it is a form similar to a harmonic series, and it is difficult to find an exact solution. But the use of Fourier transform in the Edward summation method can be done in the inverse easy space for energy calculation, can approach the exact solution. This article mainly introduces the principle of Fourier transform and the corresponding Python code implementation.

DFT Principle

The DFT calculation is essentially a matrix operation:

\[y_k=\sum_{n=0}^{N-1}x_ne^{-j\frac{2\pi nk}{N}},0\leq k\leq N-1 \]

\[x_n=\frac{1}{N}\sum_{k=0}^{N-1}y_ke^{j\frac{2\pi nk}{N}},0\leq n\leq N-1 \]

If written in the form of a matrix, that is:

\[\left[\begin{matrix} y_1\\ y_2\\ y_3\\ ...\\ y_{N-1} \end{matrix}\right]=\left[\begin{matrix} 1&&1&&1&&...&&1\\ 1&&e^{-j\frac{2\pi}{N}\cdot1}&&e^{-j\frac{2\pi}{N}\cdot2}&&...&&e^{-j\frac{2\pi}{N}\cdot(N-1)}\\ 1&&e^{-j\frac{2\pi}{N}\cdot2}&&e^{-j\frac{2\pi}{N}\cdot4}&&...&&e^{-j\frac{2\pi}{N}\cdot2(N-1)}\\ ...&&...&&...&&...&&...\\ 1&&e^{-j\frac{2\pi}{N}\cdot(N-1)}&&e^{-j\frac{2\pi}{N}\cdot2(N-1)}&&...&&e^{-j\frac{2\pi}{N}\cdot(N-1)(N-1)} \end{matrix}\right]\left[ \begin{matrix} x_1\\ x_2\\ x_3\\ ...\\ x_{N-1} \end{matrix}\right] \]

Similarly, the matrix form of the inverse Fourier transform is:

\[\left[ \begin{matrix} x_1\\ x_2\\ x_3\\ ...\\ x_{N-1} \end{matrix}\right] =\left[\begin{matrix} 1&&1&&1&&...&&1\\ 1&&e^{j\frac{2\pi}{N}\cdot1}&&e^{j\frac{2\pi}{N}\cdot2}&&...&&e^{j\frac{2\pi}{N}\cdot(N-1)}\\ 1&&e^{j\frac{2\pi}{N}\cdot2}&&e^{j\frac{2\pi}{N}\cdot4}&&...&&e^{j\frac{2\pi}{N}\cdot2(N-1)}\\ ...&&...&&...&&...&&...\\ 1&&e^{j\frac{2\pi}{N}\cdot(N-1)}&&e^{j\frac{2\pi}{N}\cdot2(N-1)}&&...&&e^{j\frac{2\pi}{N}\cdot(N-1)(N-1)} \end{matrix}\right] \left[\begin{matrix} y_1\\ y_2\\ y_3\\ ...\\ y_{N-1} \end{matrix}\right] \]

If the notation parameter$W_{N,n,k}=e^{-j\frac{2\pi}{N}nk}$, then its conjugate$W_{N,n,k}^*=e^{j\frac{2\pi}{N}nk}$is the parameter of the inverse Fourier transform. And by the nature of the complex function, the parameter is periodic:$W_{N,n+N,k}=W_{N,n,k+N}=W_{N,n,k}$, the conjugate parameter is the same. Finally there is another very important property:$W_{N/m,n/m,k}=W_{N/m,n,k/m}=W_{N,n,k}$, based on this property, it is possible to turn large-scale operations into small-scale computations. Without considering these parametric properties, we can do a preliminary simple implementation of DFT using Python.

Initial Python implementation

No optimizations have been made here, just an example:

import numpy as np

def dft(x):
    y = np.zeros_like(x, dtype=np.complex64)
    N = [0]
    for k in range(N):
        y[k] = (x * (-1j*2**k*(N)/N))
    return y

def idft(y):
    x = np.zeros_like(y, dtype=np.float32)
    N = [0]
    for n in range(N):
        x[n] = ((y * (1j*2**n*(N)/N)) / N)
    return x

N = 128
x = (N).astype(np.float32)
y0 = dft(x)
y1 = (x)
print ((y0, y1))

yr = (N).astype(np.float32)
yi = (N).astype(np.float32)
y = yr + 1j*yi
x0 = idft(y)
x1 = (y).real
print ((x0, x1))
# True
# True

Both of the outputs areTrue, which also means that this calculation is fine.

FFT Fast Fourier Transform

First we organize all the parameter-related optimization points:

\[\left\{ \begin{matrix} W_{N,n+N,k}=W_{N,n,k+N}=W_{N,n,k}\\ W_{N/m,n/m,k}=W_{N/m,n,k/m}=W_{N,n,k}\\ W_{N,\beta,\frac{N}{2\beta}}=W^*_{N,\beta,\frac{N}{2\beta}}=-1\Rightarrow W_{N,n,k}\cdot W_{N,\beta,\frac{N}{2\beta}}=-W_{N,n,k}\\ W_{N,N-n,k}=W_{N,N,k}\cdot W_{N,-n,k}=W_{N,-n,k}=W_{N,n,N-k} \end{matrix} \right. \]

At this point if we take the original input$x_n$Split into odd and even groups (if the total number N is not even, you can generally padding the input array):

\[\left\{ \begin{matrix} x_{2r}\\ x_{2r+1} \end{matrix},0\leq r\leq \frac{N}{2}-1 \right. \]

Then there is:

\[y_k=\sum_{n=0}^{N-1}x_ne^{-j\frac{2\pi nk}{N}}=\sum_{r=0}^{\frac{N}{2}-1}x_{2r}W_{N,2r,k}+\sum_{r=0}^{\frac{N}{2}-1}x_{2r+1}W_{N,2r+1,k} \]

If we take the$x_{2r}$cap (a poem)$x_{2r+1}$viewed as two separate input data, then the above decomposition can be further optimized:

\[\begin{align*} y_k&=\sum_{r=0}^{\frac{N}{2}-1}x_{2r}W_{N,2r,k}+\sum_{r=0}^{\frac{N}{2}-1}x_{2r+1}W_{N,2r+1,k}\\ &=\sum_{r=0}^{\frac{N}{2}-1}x^{(odd)}_{r'}W_{\frac{N}{2},r',k}+W_{N,1,k}\sum_{r=0}^{\frac{N}{2}-1}x^{(even)}_{r'}W_{\frac{N}{2},r',k}\\ &=y^{(odd)}_k+W_{N,1,k}y^{(even)}_k \end{align*} \]

The same reasoning can be obtained:

\[\begin{align*} y_{k+\frac{N}{2}}&=y^{(odd)}_{k+\frac{N}{2}}+W_{N,1,k+\frac{N}{2}}y^{(even)}_{k+\frac{N}{2}}\\ &=y^{(odd)}_{k+\frac{N}{2}}-W_{N,1,k}y^{(even)}_{k+\frac{N}{2}} \end{align*} \]

This is known as the butterfly operation (image from reference link):

The significance of this operator is that if we originally did a Fourier transform of $2^N$ point data, using the original DFT operation we would need to do $2^{2N}$ multiplications and $2^N(2^N-1)$ additions, but this method reduces the computation to $2\cdot2^N+2^{\frac{N}{2}}$ multiplications and $2^{\ frac{N}{2}}(2^\frac{N}{2}-1)$ additions. Do a decomposition, the complexity from $ O (2 ^ {2N}) $ down to $ O (2 ^ N) $ (Note: here $ N $ with the previous use of the total number of data points is not a meaning, where the $ N $ refers to the total number of data points is an integer multiple of 2, just the two expression habits are commonly used $ N $). Relevant code to achieve the following:

import numpy as np

def dft(x):
    y = np.zeros_like(x, dtype=np.complex64)
    N = [0]
    for k in range(N):
        y[k] = (x * (-1j*2**k*(N)/N))
    return y

def dft2(x):
    y = np.zeros_like(x, dtype=np.complex64)
    N = [0]
    for k in range(N//2):
        c1 = (-1j*2**k*(N//2)/(N//2))
        c2 = (-1j*2**k/N)
        y1 = (x[::2] * c1)
        y2 = (x[1::2] * c1)
        y[k] = y1 + c2 * y2
        y[k+N//2] = y1 - c2 * y2
    return y

N = 128
x = (N).astype(np.float32)
y0 = dft2(x)
y1 = (x)
print ((y0, y1))
# True

The output of the run isTrue, indicating that the results are consistent. Note that the code here does not take padding into account, and cannot be used as a formal code implementation, but merely an algorithm demonstration. Since it can be divided once, it can be divided many times, until it can not be divided, or divided to a specified parameter. This is the principle of the multiple butterfly operation:

For simplicity you can use recursion for the calculation:

import numpy as np

def dft(x):
    y = np.zeros_like(x, dtype=np.complex64)
    N = [0]
    for k in range(N):
        y[k] = (x * (-1j*2**k*(N)/N))
    return y

def dftn(x, N_cut=2):
    y = np.zeros_like(x, dtype=np.complex64)
    N = [0]
    if N > N_cut:
        y1 = dftn(x[::2])
        y2 = dftn(x[1::2])
    else:
        return dft(x)
    for k in range(N//2):
        c2 = (-1j*2**k/N)
        y[k] = y1[k] + c2 * y2[k]
        y[k+N//2] = y1[k] - c2 * y2[k]
    return y

N = 1024
x = (N).astype(np.float32)
y0 = dftn(x)
y1 = (x)
print ((y0, y1))
# True

The implementation here uses a recursive approach that combines the DFT algorithm and the butterfly operation method implemented earlier, and the results obtained are also correct. The optimization method of butterfly operation used here is the basic idea of FFT Fast Fourier Transform.

N-point fast Fourier transform

The so-called N-point FFT is actually to take only N data points to perform Fourier transform. Then there are many ways to take the data points, such as taking only the first N data points, or downsampling and then take the first N data points, and then add the window, after the operation of the window function, to do the Fourier transform of the data points in each window. The simplest way is the rectangular window, and the common ones are Hanning window and Hamming window, which are not analyzed in detail here. It is worth noting that if the downsampling method is used, the sampling rate needs to follow the Nyquist sampling theorem and be greater than twice the target frequency. especially for the periodic boundary conditions and remote interaction scenarios, the contribution of the high frequency region cannot be ignored.

As for why not use the Fourier transform of the full-domain data points, even though we can use the fast Fourier transform to reduce the computational complexity to$O(N\log N)$(Here's$N$(which is the number of data points) level is also not applicable for those scenarios with large-scale data transfer and computation, so the idea of using a reduced number of Fourier Transform points can be a good balance between performance and accuracy for most scenarios. The appearance of the form function further optimizes the problem of data leakage at the truncation.

Summary outline

This paper introduces the basic principles of discrete Fourier transform and fast Fourier transform and its corresponding Python code implementation, and compares the calculation results with the fft function integrated by numpy. In fact, there are now many mature tools for FFT computation, whether it is the fft module of scipy on CPU or the cufft dynamic link library on GPU, both of which have very good performance. But you still have to really understand the principle behind the calculation, and the related physical images, in order to use this powerful tool more appropriately.

This article was first linked to:/dechinphy/p/

Author ID: DechinPhy