Reusing a complex FFT array’s memory as a real array, “C-style”

General Usage
, ,
1

When doing a complex-to-real FFT (using FFTW.jl), the real input and complex output arrays take roughly the same amount of memory. To quote from the rfft() docstring,

If A has size (n_1, …, n_d), the result has size (n_1 ÷ 2 + 1, …, n_d).

where A is the real input array. To save memory, it’s useful to convert between the two interpretations, i. e. a real array of size (n_1, …, n_d) (leaving some extra unused memory) and a complex array of size (n_1 ÷ 2 + 1, …, n_d). In C, this is simple: A is usually a memory block from malloc(), and the interpretation just depends on how this is indexed.

My question is how to do the same thing in Julia, e. g. take the complex output array from rfft() and re-use it as a real array for some other computation. As an example, let’s take a real array A and a complex array C from rfft():

julia> A = collect(reshape(1.0:5 * 3 * 4, (5, 3, 4)))
5×3×4 Array{Float64, 3}:
[:, :, 1] =
 1.0   6.0  11.0
 2.0   7.0  12.0
 3.0   8.0  13.0
 4.0   9.0  14.0
 5.0  10.0  15.0

[:, :, 2] =
 16.0  21.0  26.0
 17.0  22.0  27.0
 18.0  23.0  28.0
 19.0  24.0  29.0
 20.0  25.0  30.0

[:, :, 3] =
 31.0  36.0  41.0
 32.0  37.0  42.0
 33.0  38.0  43.0
 34.0  39.0  44.0
 35.0  40.0  45.0

[:, :, 4] =
 46.0  51.0  56.0
 47.0  52.0  57.0
 48.0  53.0  58.0
 49.0  54.0  59.0
 50.0  55.0  60.0

julia> C = rfft(A)
       # set elements to consecutive integers
       reshape(C, :) .= collect(1.0:2:2 * length(C)) .+ 1im .* collect(2.0:2:2 * length(C))
       C
3×3×4 Array{ComplexF64, 3}:
[:, :, 1] =
 1.0+2.0im   7.0+8.0im   13.0+14.0im
 3.0+4.0im   9.0+10.0im  15.0+16.0im
 5.0+6.0im  11.0+12.0im  17.0+18.0im

[:, :, 2] =
 19.0+20.0im  25.0+26.0im  31.0+32.0im
 21.0+22.0im  27.0+28.0im  33.0+34.0im
 23.0+24.0im  29.0+30.0im  35.0+36.0im

[:, :, 3] =
 37.0+38.0im  43.0+44.0im  49.0+50.0im
 39.0+40.0im  45.0+46.0im  51.0+52.0im
 41.0+42.0im  47.0+48.0im  53.0+54.0im

[:, :, 4] =
 55.0+56.0im  61.0+62.0im  67.0+68.0im
 57.0+58.0im  63.0+64.0im  69.0+70.0im
 59.0+60.0im  65.0+66.0im  71.0+72.0im

How do I turn C into a real array with the same shape as A? This is what I came up with:

C_as_real = reshape(@view(reshape(reinterpret(eltype(A), C), :)[1:length(A)]), size(A))
5×3×4 reshape(view(reshape(reinterpret(Float64, ::Array{ComplexF64, 3}), 72), 1:60), 5, 3, 4) with eltype Float64:
[:, :, 1] =
 1.0   6.0  11.0
 2.0   7.0  12.0
 3.0   8.0  13.0
 4.0   9.0  14.0
 5.0  10.0  15.0

[:, :, 2] =
 16.0  21.0  26.0
 17.0  22.0  27.0
 18.0  23.0  28.0
 19.0  24.0  29.0
 20.0  25.0  30.0

[:, :, 3] =
 31.0  36.0  41.0
 32.0  37.0  42.0
 33.0  38.0  43.0
 34.0  39.0  44.0
 35.0  40.0  45.0

[:, :, 4] =
 46.0  51.0  56.0
 47.0  52.0  57.0
 48.0  53.0  58.0
 49.0  54.0  59.0
 50.0  55.0  60.0

julia> A == C_as_real
true

Explanation:

  • reinterpret(eltype(A), C): Reinterpret the complex array as real
  • reshape(..., :): Reshape the array as 1D (to be able to get rid of the extra elements)
  • @view(...[1:length(A)]): Discard extra memory
  • reshape(..., size(A)): Reshape back to the original shape

Am I overcomplicating things, or is this how it needs to be done?

2

There have been a few attempts to add in-place real-to-complex transforms to FFTW.jl (for example here) which for some reason have not been merged to this day. There is also the RealFFTs.jl package which seems to serve the same purpose. I haven’t tried it myself but that is probably your best bet if you want an easy solution.

2 Likes
3

Yeah, I initially wanted to do an in-place real-to-complex FFT, and instead ended up pre-allocating the complex output array and using that via mul!, which is OK, except that I now have a complex array where I normally need a real one. So I guess if I can use RealFFTs.jl to do the FFT in-place, I can then allocate the second array as real to begin with and only use it for the following computations.

4

I guess these lines from RealFFTs.jl show how to do what I asked originally in a more C-like (and less safe) way:

1 Like
5

You might be able to do this safely using the Base.reinterpret function, specifically the method

reinterpret(reshape, T, A::AbstractArray{S}) -> B
1 Like
6

I did use reinterpret in the original post. I also saw that specific method, but it does something completely different:

julia> reinterpret(reshape, Float64, C)
2×3×3×4 reinterpret(reshape, Float64, ::Array{ComplexF64, 3}) with eltype Float64:
…
7

Sorry, I didn’t read your original post carefully enough.

8

No worries :slight_smile:

1 Like