I have 3D real data, say, 64 * 64 * 64 real data points, and I want to do a 3D Fourier transform. Specifically, I want to do DFT (exponential transform, containing both sine and cosine terms) on two of the three dimensions and DCT (discrete cosine transform) on the third dimension. The reason why I want to do this is because I'm solving partial differential equations (PDEs) where I want to impose periodic boundary condition in two directions while Neumann boundary condition in the third direction. To achieve this, I want to use fftw_plan_r2r_3d in FFTW, in which the DFT is done by R2HC who outputs data in the halfcomplex format. However, FFTW manual explicitly mentions that multiple R2HC is not recommended:
So now I need to figure out how to interpret the output and reconstruct the correct coefficients after the transform.
Just to make it a bit simpler, I will first ignore the DCT in the third dimension and just consider now 2D R2HC transform. My way to try to understand it is to compare the 2D R2HC transform to a 2D R2C transform. Take 4*4 case as an example, using some synthesized real space data and transform it using R2C without truncation (FFTW truncates it because there is some symmetry if the real space data is purely real and thus some data are redundant) we have:
where the data points related through complex conjugate are boxed by the same colour (except for the first row and column for which I did not draw box). As we can see here, due to the symmetry, you only need one of the two data points with the same colour to get all the information required to reconstruct the whole data. Due to this reason, FFTW truncates the last dimension of the data outputted by multidimensional r2c transform from N to N/2+1, and in our simple example that is:
For R2HC transform, the halfcomplex format means that the truncated coefficients are stored in a double-type array where the real and imaginary parts are separated and arranged in some order (introduced in FFTW manual). The data outputted by 2D R2HC looks like:
And now if I compare 2D R2HC to 2D r2c, problem comes. In 2D R2HC transform, the output coefficients are truncated in both dimensions, and if we try to visualize it in the equivalent way to r2c transform, it's like:
We can see that the two yellow (or orange) boxes are both gone. In this case, when I try to use the output from R2HC to reconstruct the whole coefficients, I can sorta successfully get all the values by combining different pieces except for the one in yellow box (I guess, in my understanding, that makes sense because two yellow boxes are both gone). However, if I call the reverse transform (2D HC2R) and normalize it, the real space data is fully recovered, which means all the information are still preserved. So that means I should be able to reconstruct the yellow one as well it's just I don't know how (due to my lack of understanding of the symmetry in dft).
Now, I wonder how I could reconstruct all the correct coefficients from the raw output of 2D R2HC transform so that I can then use them in solving the PDEs.
