pythonnumpyfftnumerical-methodspyfftw

Calculation of Laplacian in real pyFFTW

For the forward (multidimensional) FFTW algorithm you can specify that the input numpy.ndarray is real, and the output should be complex. This is done when creating the byte-aligned arrays that go in the arguments of the fft_object:

import numpy as np
import pyfftw

N = 256  # Input array size (preferrably 2^{a}*3^{b}*5^{c}*7^{d}*11^{e}*13^{f}, (e+f = 0,1))
dx = 0.1  # Spacing between mesh points
a = pyfftw.empty_aligned((N, N), dtype='float64')
b = pyfftw.empty_aligned((N, N//2+1), dtype='complex128')
fft_object = pyfftw.FFTW(a, b, axes=(0, 1), direction='FFTW_FORWARD')

The output array is not symmetric and the second axis is truncated up to the positive frequencies. For the complex FFT you can compute the laplacian with the following np.ndarray

kx, ky = np.meshgrid(np.fft.fftfreq(N, dx), np.fft.fftfreq(N, dx))  # Wave vector components
k2 = -4*np.pi**2*(kx*kx+ky*ky)  # np.ndarray for the Laplacian operator in "frequency space"

How would it be done in the truncated case? I thought about using:

kx, ky = np.meshgrid(np.fft.fftfreq(N//2+1, dx), np.fft.fftfreq(N, dx))  # The axes conven-
#                                                                        tions are different

But, would this really work? It seems like it is neglecting the negative frequencies in the "y" direction.

Solution

  • I'm not familiar with pyfftw, but with the numpy.fft module it would work just fine (assuming you use rfftfreq as mentioned in the comments).

    To recap: for a real array, a, the fourier transform, b, has a Hermtian-like property: b(-kx,-ky) is the complex conjugate of b(kx,ky). The real version of the forward fft discards (most of) the redundant information by omitting the negative kys. The real version of the backward fft assumes that the values at the missing frequencies can be found by complex conjugating the appropriate elements.

    If you had used the complex fft and kept all frequencies, -k2 * b would still have the Hermitian-like property. So the assumption made by the real backward fft still holds and would give the correct answer.

    I guess with pyfftw it will work just fine provided that you specify a float64 array of the correct size for the output for the direction=FFT_BACKWARD case.