We now implemented FourierTools.jl a first version of which will appear soon. I want to use this discussion to outline how a proper padding and cutting in Fourier space (also mixed versions) needs to be done.
- We want: Reversibility for cases where we do not reduce sizes.
- We need to not change the integrated power spectrum
The latter part is crucial when it comes to dealing with even sized dimensions in FFT and RFFT of real valued data alike. The lowest (biggest negative) frequency stores in these cases the sum of amplitude of the negative highest frequency PLUS the corresponding value at the aliased positive highest frequency (e.g. [-kxmax, ky, kz] + [+kxmax, ky, kz]). Note that this value is generally NOT real except for special choices of ky and kz.
When changing the region of interest in FFT- (or RFFT-) space the algorithm consists of two steps:
fix_dim_before(Array, dim)-
fix_dim_after(Array, dim)
Both of these routines need to be run for each dimension consecutively. The first before the size change and the second after it.
fix_dim_before(Array, dim) is only active for a dimension if it is even-sized and enlarged
In this case the highest frequency subslice(dim) needs to be halfed in value and this half also copied (without conjugation) to the equivalent positive frequency (e.g. [+kxmax, ky, kz]).
fix_dim_after(Array, dim) is only active for a dimension if it is even-sized and reduced
Here we need to copy the lowest frequency which is just about dropped (e.g. [+kxmax, ky, kz]) and sum it to the other side (here [-kxmax, ky, kz]))
We implemented this for array views (overloading getindex, such that no extra memory allocation is needed.
A remark on just taking the real part and ignoring such (rather complicated) actions:
An even lowest frequency in FFTs (and RFFTs) store both, the positive and negative highest frequencies. Going back to real-space after padding in Fourier-space without copying will then drop some of the energy at the moment when the real part is taken. When reverting this situation (cutting away the previously introduced zeros) it is clear to see that this real-space energy is still lost.
An interesting consequence is that also some other functions need modifications on which we are currently working: For example: multiplication with a complex exponential needs a modification to the lowest frequency along an even dim.
It seems to be that also calculating the Energy or Power spectrum needs a modification to such frequencies: sum((a+a)^2) !== 2*sum(a^2)
The should be fairly easy to check via a comparison to the real-space integral of the square of the signal.
If array sizes are fairly large, all of these effects are small, but wouldn’t it be nice to get this correct?
This is what FourierTools.jl will be trying to achieve.