FFTW.jl icon indicating copy to clipboard operation
FFTW.jl copied to clipboard

Why is there no plan_rrft! ?

Open fmeirinhos opened this issue 4 years ago • 4 comments

Is there any fundamental reason why the inplace plan_rfft! / plan_irfft! are not implemented?

fmeirinhos avatar Apr 16 '20 10:04 fmeirinhos

Out of my head, rFFTPlans perform real to complex Fourier Transform. So the input is an array of real numbers and the Fourier Transform is stored in an array of complex numbers. Therefore, there is no way ( or no obvious way ) of writing the results onto the input array, as they have different types.

plan_irfft performs the opposite operation, complex to real Fourier Transform. The same logic applies here.

For real inputs you also have plan_r2r (real to real), and this function has its inplace counterpart, plan_r2r!. I don't use these functions though, because there are many types of r2r transform, and I have never spent the time to study which applies where.

Marc-3d avatar Apr 16 '20 13:04 Marc-3d

I was under the impression that inplace plans would compute the transforms in a pre-determined location (since I thought FFT operations could not actually be done in place!).

For repeated rfft / irfft operations into a predetermined array then I guess one just needs to call unsafe_execute! instead of the * operator.

fmeirinhos avatar Apr 18 '20 14:04 fmeirinhos

I will quote the FFTW3 manual about inplace transforms, so I don't make anything up:

in and out point to the input and output arrays of the transform, which may be the same (yielding an in-place transform).

Where in and out refer to the plan creation signature: fftw_plan_dft_2d(int n0, int n1, fftw_complex *in, fftw_complex *out, ... )

Later in the manual it also explains that plans can be reused on different input arrays, as long as certain conditions are met. For instance, they must have the same size as the original array the plan was created on. The manual calls it "New-array Execute Functions".

Regarding the * operator, it is overloaded in fft.jl and eventually calls unsafe_execute!.

function *(p::r2rFFTWPlan{T,K,true}, x::StridedArray{T}) where {T,K}
    assert_applicable(p, x)
    unsafe_execute!(p, x, x)
    return x
end

The code above is for in-place r2r plans, indicated by true in {T,K,true}.

Marc-3d avatar Apr 21 '20 13:04 Marc-3d

There shouldn't be any reason you can't do an in-place (i)rfft. FFTW supports this. The trick is that the real input must be padded by either one or two values depending on whether it is odd or even. This would make the function a little tricky to use in Julia, but I still think it would be very valuable. Regarding data types, the Julia reinterpret function can handle this. Maybe an in-place real p * x would return reinterpret(Complex(typeof(x)), p*x).

A few references in the FFTW docs: http://www.fftw.org/fftw3_doc/Real_002ddata-DFTs.html#Real_002ddata-DFTs http://www.fftw.org/fftw3_doc/Real_002ddata-DFT-Array-Format.html#Real_002ddata-DFT-Array-Format

btmit avatar Sep 19 '20 15:09 btmit