FFTW / fftw3

DO NOT CHECK OUT THESE FILES FROM GITHUB UNLESS YOU KNOW WHAT YOU ARE DOING. (See below.)
GNU General Public License v2.0
2.72k stars 661 forks source link

Docs for Half-complex should describe connection with sine/cosine expansion #319

Open dlfivefifty opened 1 year ago

dlfivefifty commented 1 year ago

I'm surprised the docs for "Half-complex" don't explain that the transform gives the coefficients in cos/sin series:

julia> m = 5; θ = range(0,2π; length=m+1)[1:m];

julia> FFTW.r2r(cos.(2θ), FFTW.R2HC) # only non-zero is in entry corresponding to cos(2θ)
5-element Vector{Float64}:
  0.0
 -2.220446049250313e-16
  2.5
 -6.709199951388644e-16
 -2.310772834960631e-17

julia> FFTW.r2r(sin.(θ), FFTW.R2HC) # only non-zero is in entry corresponding to sin(θ)
5-element Vector{Float64}:
  1.1102230246251565e-16
 -2.139456371091744e-16
  1.5843448587791657e-16
  0.0
 -2.5

This is particular important in multiple dimensions since it gives precisely the coefficients in the tensor product basis, which is clearer than the explanation in the docs:

julia> FFTW.r2r(sin.(θ) .* cos.(2φ'), FFTW.R2HC) # only non-zero is in entry corresponding to sin(θ)*cos(2φ)
5×6 Matrix{Float64}:
  0.0           0.0           3.33067e-16  0.0          0.0           0.0
  0.0           0.0          -6.41837e-16  0.0          0.0           0.0
  0.0           0.0           4.75303e-16  0.0          0.0           0.0
  1.11022e-16  -5.55112e-17  -5.55112e-17  1.11022e-16  9.61481e-17  -9.61481e-17
 -1.77636e-15  -1.11022e-16  -7.5          2.22045e-16  3.07674e-15  -2.88444e-15
matteo-frigo commented 1 year ago

It depends upon what you mean by "sin/cos series". In a naive interpretation of your comment, the coefficient in the transformed array should be 5 (or perhaps 1), but definitely not 2.5. In addition, if you compute FFTW.r2r(cos(0 * θ), FFTW.R2HC), you get 5 as coefficient, whereas all other frequencies have coefficient 2.5.

The position that FFTW takes is that sines and cosines are a historical relic and only complex exponentials exist. Your cos(2 theta) input is just the sum of two complex exponentials, from which it is "obvious" that the coefficient should be 2.5 and not 5 (where the other 2.5 at negative frequency is redundant and not stored). cos(0 theta) is one complex exponential. This convention at least has the advantage that frequency 0 is not in any way special.

Anyway, that's why things are the way they are, but we are not fascist about this. If you have a suggested text for that section in the manual, we'll be happy to talk about it.

dlfivefifty commented 1 year ago

That’s an odd position since sines/cosines preserve realness of the function, which is very important in many PDE applications. Yes in 1D you can reinterpret it as a symmetry property of the coefficients but as the current docs make clear it doesn’t generalise to higher dimensions in a natural way.

I’ll try to come up with a suitable text.