# Multidimensional Real Fast Fourier Transform

**URL:** <https://discourse.julialang.org/t/multidimensional-real-fast-fourier-transform/56498>\
**Category:** General Usage\
**Tags:** fftw, dsp\
**Created:** [March 4, 2021, 8:09pm UTC](https://discourse.julialang.org/t/multidimensional-real-fast-fourier-transform/56498 "2021-03-04T20:09:46Z")\
**Posts on this page:** 9\
**Page:** 1

<div class="post-metadata">

**Author:** ![tiago](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tiago/32/22567_2.png) [@tiago](https://discourse.julialang.org/u/tiago)\
**Post date:** [March 4, 2021, 8:09pm UTC](https://discourse.julialang.org/t/multidimensional-real-fast-fourier-transform/56498/1 "2021-03-04T20:09:46Z")

</div>

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:

```julia
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
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

```julia
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
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

---

<div class="post-metadata">

**Author:** ![jgidi](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jgidi/32/26206_2.png) [@jgidi](https://discourse.julialang.org/u/jgidi)\
**Post date:** [March 4, 2021, 11:16pm UTC](https://discourse.julialang.org/t/multidimensional-real-fast-fourier-transform/56498/2 "2021-03-04T23:16:53Z")

</div>

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!

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [March 4, 2021, 11:50pm UTC](https://discourse.julialang.org/t/multidimensional-real-fast-fourier-transform/56498/3 "2021-03-04T23:50:40Z")

</div>

> [@tiago](#):
>
> 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` .

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](http://fftw.org/fftw3_doc/Multi_002ddimensional-Transforms.html#Multi_002ddimensional-Transforms):

 ![image](https://global.discourse-cdn.com/julialang/original/3X/7/3/73f9509e7245b40957f79013b03d43c96e476524.png)

> [@tiago](#):
>
> _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?

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

---

<div class="post-metadata">

**Author:** ![rafael.guerra](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rafael.guerra/32/216610_2.png) [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)\
**Post date:** [March 5, 2021, 7:50am UTC](https://discourse.julialang.org/t/multidimensional-real-fast-fourier-transform/56498/4 "2021-03-05T07:50:43Z")

</div>

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](https://www.researchgate.net/publication/3432376_2-D_symmetry_Theory_and_filter_design_applications)_

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.

---

<div class="post-metadata">

**Author:** ![tiago](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tiago/32/22567_2.png) [@tiago](https://discourse.julialang.org/u/tiago)\
**Post date:** [March 5, 2021, 10:36am UTC](https://discourse.julialang.org/t/multidimensional-real-fast-fourier-transform/56498/7 "2021-03-05T10:36:09Z")

</div>

@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 `rfft`the 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!

---

<div class="post-metadata">

**Author:** ![rafael.guerra](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rafael.guerra/32/216610_2.png) [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)\
**Post date:** [March 5, 2021, 10:47am UTC](https://discourse.julialang.org/t/multidimensional-real-fast-fourier-transform/56498/8 "2021-03-05T10:47:52Z")

</div>

@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)](http://www.fftw.org/fftw-paper-ieee.pdf)_

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

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [March 5, 2021, 1:28pm UTC](https://discourse.julialang.org/t/multidimensional-real-fast-fourier-transform/56498/9 "2021-03-05T13:28:56Z")

</div>

> [@tiago](#):
>
> in the majority of the cases I have symmetrical real data, so this seems promising!

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

---

<div class="post-metadata">

**Author:** ![fgerick](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fgerick/32/13228_2.png) [@fgerick](https://discourse.julialang.org/u/fgerick)\
**Post date:** [March 5, 2021, 1:46pm UTC](https://discourse.julialang.org/t/multidimensional-real-fast-fourier-transform/56498/10 "2021-03-05T13:46:12Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![tiago](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tiago/32/22567_2.png) [@tiago](https://discourse.julialang.org/u/tiago)\
**Post date:** [March 5, 2021, 4:14pm UTC](https://discourse.julialang.org/t/multidimensional-real-fast-fourier-transform/56498/11 "2021-03-05T16:14:52Z")

</div>

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