Multidimensional Real Fast Fourier Transform

General Usage
,
1

Hello,
I’m using fourier transformations to solve a partial differential equation in two dimensions. Since performance is super important in my case and I only deal with real data, so i’m using the pre-computed plan of the rfft, plan_rfft and the respective inverse, plan_irfft.

The big advantage of using a rfft instead of the normal fft, it’s the fact that we only need to compute half of the spectrum, N/2 + 1, due to the hermitian symmetry. Which is big a plus, since I need to perform computations in the Fourier domain, this implies lighter computations.

But, if we go to the fftw.jl package documentation, we can see that the half of the dimension is only valid for the first dimension of the array:

help?> rfft
search: rfft rfftfreq irfft brfft plan_rfft plan_irfft

  rfft(A [, dims])

  Multidimensional FFT of a real array A, exploiting the      
  fact that the transform has conjugate symmetry in order to  
  save roughly half the computational time and storage costs  
  compared with fft. If A has size (n_1, ..., n_d), the       
  result has size (div(n_1,2)+1, ..., n_d).

As example:

julia> using FFTW

julia> A=rand(Int64,9,9);

julia> rfft(A)
5×9 Array{Complex{Float64},2}:
  3.62383e18+0.0im         …  -8.17985e18+5.85874e18im
 -2.70232e18+5.41946e19im     -2.63532e19+1.40708e19im
  6.06796e19-1.54719e19im     -1.78341e19+6.2371e18im 
  2.68829e19+8.20984e18im      -7.1641e19+6.47834e19im
  4.83938e19+7.09874e19im      1.23248e19-1.17264e19im

And if we check what the documentation says about multidimensional fft, we see that a multidimensional fft is just a ´fft´applied to each dimension

help?> fft
...
  A multidimensional FFT simply performs this operation       
  along each transformed dimension of A.
...

But this is not what is implemented when dealing with the real transformation, correct? Only computes half of the spectrum in the first dimension, as we can see below:

julia> A=rand(Int64,4,4);

julia> rfft(A)
3×4 Array{Complex{Float64},2}:
 2.49889e19+0.0im          1.49133e19-8.09715e18im  2.81949e19+0.0im          1.49133e19+8.09715e18im
 1.76655e19+7.75576e18im  -5.13541e18-2.90954e19im  1.62245e19+1.41345e19im  -2.73817e18+1.48001e19im
 -1.5568e19+0.0im         -3.64835e18+1.08543e19im   2.4693e19+0.0im         -3.64835e18-1.08543e19im

julia> fft(A)
4×4 Array{Complex{Float64},2}:
 2.49889e19+0.0im          1.49133e19-8.09715e18im  2.81949e19+0.0im          1.49133e19+8.09715e18im
 1.76655e19+7.75576e18im  -5.13541e18-2.90954e19im  1.62245e19+1.41345e19im  -2.73817e18+1.48001e19im
 -1.5568e19+0.0im         -3.64835e18+1.08543e19im   2.4693e19+0.0im         -3.64835e18-1.08543e19im
 1.76655e19-7.75576e18im  -2.73817e18-1.48001e19im  1.62245e19-1.41345e19im  -5.13541e18+2.90954e19im

My question is, how can I implement the rfft in such way that I deal with matrices with dimension (N/2 + 1) x (N/2 + 1) and not the (N/2 + 1) x N returned by rfft. Will I have to add a method to plan_rfft? If so, I’m afraid that will be too complex to deal with.

Thanks,
Tiago

2

Hello!

I am sorry to tell you this, but what you propose can’t be done.

As you said, a multidimensional fft consists in the successive execution of the Fourier transform along each dimension of the array. The problem here is that, after the first Fourier transform is applied to the real data, the result is now complex, and the following Fourier transforms don’t have hermitian symmetry. This is the reason why the multidimensional rfft may take advantage of the hermitian symmetry and compute a half of the spectrum, but only along the first dimension transformed.

I hope this explanation helps you!

4 Likes
3

You can’t, because that’s not how a multidimensional FFT of real inputs works. You can only save ≈half the data and computation regardless of the dimensionality. That is, you can only cut one dimension in half, and FFTW for convenience chose the first dimension (in Julia column-major order, corresponding to the last dimension on C).

See this explanation from the FFTW manual:

image
image2170×558 166 KB

Conceptually it is transforming along each dimension, and then discarding a redundant half of the output.

In practice, it first performs the rfft along the first dimension, which produces complex data of half the size, and then transforms along the remaining dimensions — but since the result of the first transform was complex, the remaining computations are ordinary complex FFTs.

5 Likes
4

What do you think about the following reference:
Reddy, H. C., Khoo, I. H. and Rajan, P. K. [2003] 2-D symmetry: theory and filter design applications, IEEE Circuits and Systems Magazine, vol. 3, no. 3, pp. 4-33

They decompose the input signal into symmetrical and anti-symmetrical components and process each of them using appropriate symmetry properties.
For 2D FFT computation, the authors claim that while for general non-symmetric input signals the method may not reduce the overall complexity, it should still facilitate parallel processing.

1 Like
7

@jgidi Yes, it helped a lot your explanation. I had already realized that the second rfft couldn’t be performed after the first because of that, if we by hand do rfft(â,2) where A is the output of rfft(A) gives an expected error.

I thought adding a method to the rfft in such way that I could compute only the first half of the spectrum in the second dimension, by exploiting the hermitian symmetry. But I totally forgot that the fact you don’t have real input anymore after applying the rfftthe first time.

Thank you, it helped a lot!

@rafael.guerra Thanks for the reference, I didn’t knew this one. I’ll read this carefully in order to understand if I can take advantage, in the majority of the cases I have symmetrical real data, so this seems promising!
Once again, thanks!

@stevengj Thank you for that posting that manual, I didn’t know about it I’ll read more carefully! It’s curious that it’s consider the possibility to discard the second half of the spectrum in the other dimensions, but concluded that it was not computationally efficient.
Really elucidated me!

1 Like
8

@tiago, you might have overlooked it but in that site they recommend for more information regarding FFTW, to see the paper:
“The Design and Implementation of FFTW3,” by M. Frigo and S. G. Johnson, which was an invited paper in Proc. IEEE 93 (2), p. 216 (2005)

As one can see, @stevengj is the co-author and an authority in this topic.

1 Like
9

(An FFT of real-symmetric data is also known as a discrete cosine transform. In this case you can in principle save another factor of two for each symmetric direction. Similarly for anti-symmetric real data with discrete sine transforms.)

1 Like
10

You might be interested in GitHub - FourierFlows/FourierFlows.jl: Tools for building fast, hackable, pseudospectral partial differential equation solvers on periodic domains, if that fits your problem equation.

2 Likes
11

@rafael.guerra I just read the paper and thank you for the reference, it really helped me with some questions that I had.

@stevengj Thanks! I will implement in my code the plan_dct and plan_idct to see if I gain more computational speed, of what I’ve been reading, it seems to me that it is the case, thanks once more for the help, it’s really helping me!!

@fgerick After reading the documentation example and I can say that is an interesting package, with a cool idea. But it doesn’t fit in my case, actually I’m implementing a specific numerical to algorithm to solve a quite specific equation.