Real fft of a 3D array with FFTW

General Usage
1

Hi,

I need to take the 3D fft of an array of real data. I’d like to do that using rfft() from the FFTW library rather than the usual fft() to speed up the process and keep memory requirements low.
Now my data is stored in a 3D array, where I have sampled my function f(x,y,z) in a regular
grid, but such that the interval is symmetric around the origin, so x ranges from -x0 to +x0,
y from -y0 to +y0, and z from -z0 to +z0.
I want to get the result also in a symmetric grid of k values with kx in -k0x to +k0x, ky in -k0y tp +k0y, and kz in -k0z to +k0z. All values of x0,y0,z0,kx0,ky0,kz0 are set in advance and all grids have equidistant points.

Now the problem is that I know I must be swapping parts of the input array before sending
that to rfft(), which is done with fftshift(), but I seem to fail to understand the logic behind it
so I get nonsensical results.

So does anybody know how this is done in with rfft()? Please notice that I’m assuming a simple Euler integration to go from the Fourier Transform to the discrete version, and am assuming that f(x,y,z) goes to zero well before reaching the end points of the x,y,z grids.

Best regards and thanks,

Ferran.

2

What is it about fftshift that is confusing? It simply applies a circular shift to each dimension so that whatever was at the edges is now centered, and vice versa.

julia> x = [Int(i∈(1, 4) && j∈(1,4)) for i=1:4, j=1:4]
4×4 Matrix{Int64}:
 1  0  0  1
 0  0  0  0
 0  0  0  0
 1  0  0  1

julia> fftshift(x)
4×4 Matrix{Int64}:
 0  0  0  0
 0  1  1  0
 0  1  1  0
 0  0  0  0

Same could be achieved manually using circshift:

julia> vcat((circshift(hcat((circshift(x[:,j], 2) for j=1:4)...)[i,:], 2)' for i=1:4)...)
4×4 Matrix{Int64}:
 0  0  0  0
 0  1  1  0
 0  1  1  0
 0  0  0  0

julia> x = rand(4, 4);

julia> fftshift(x) == vcat((circshift(hcat((circshift(x[:,j], 2) for j=1:4)...)[i,:], 2)' for i=1:4)...)
true
3

The DFT (what FFTs compute) corresponds to periodic boundary conditions (think “Fourier series” not “Fourier transform”), but this means that the interval is not quite symmetrical (for an even number of points) since you only want the periodic boundary point once, not twice.

In general you have to think carefully about the definitions of the transforms if you want to approximate a continuous Fourier transform by a DFT. There are lots of books on discrete-time signal processing that cover these topics.

What are you trying to accomplish with the Fourier coefficients?