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