如何从FFTW的2D R2HC变换结果重构所有傅里叶系数？

核心背景 #

我有64×64×64的3D实数据，要做混合傅里叶变换：两个维度用带正弦/余弦项的DFT（对应周期性边界条件），第三个维度用DCT（对应诺依曼边界条件），用来求解偏微分方程（PDE）。计划用FFTW的fftw_plan_r2r_3d实现，其中DFT用R2HC变换，但FFTW手册明确不推荐多次使用R2HC，因此必须搞懂它的半复数输出格式，重构出完整的傅里叶系数。

先从2D R2HC的简化场景入手分析。

2D R2HC输出的系数重构规则 #

2D实数据的傅里叶系数天生满足共轭对称性：对于尺寸为$N \times N$的输入，任意傅里叶系数$F(p,q)$满足：
$$F(p,q) = \overline{F\left((N-p)%N, (N-q)%N\right)}$$
R2HC变换正是利用这个对称性压缩存储量，所有必要信息都在输出数组里，所谓“缺失”的系数只是冗余的对称项，完全可以通过已有数据推导出来。

以4×4数据为例，R2HC输出是和输入同尺寸的double数组，重构完整复数系数的步骤如下：

1. 初始化一个$N \times N$的复数数组F。
2. 填充非冗余区域（R2HC直接输出的部分）：
   • 当$p=0$且$q=0$：$F[0][0] = \text{hc_out}[0][0] + 0i$（直流分量，虚部为0）
   • 当$p=0$且$0 < q < N/2$：$F[0][q]$的实部取hc_out[0][q]，虚部取hc_out[0][N-q]；同时根据对称性，$F[0][N-q] = \overline{F[0][q]}$
   • 当$0 < p < N/2$且$q=0$：$F[p][0]$的实部取hc_out[p][0]，虚部取hc_out[N-p][0]；$F[N-p][0] = \overline{F[p][0]}$
   • 当$0 < p < N/2$且$0 < q < N/2$：$F[p][q]$的实部取hc_out[p][q]，虚部取hc_out[N-p][N-q]；对应的对称项$F[N-p][N-q] = \overline{F[p][q]}$，$F[p][N-q] = \overline{F[N-p][q]}$
   • 当N为偶数时，$p=N/2$或$q=N/2$的位置（比如4×4中的(2,0)、(0,2)、(2,2)），系数虚部为0，直接取hc_out对应位置的值作为实部即可。
3. 所有不在非冗余区域的系数，直接通过共轭对称性推导：$F[p][q] = \overline{F\left((N-p)%N, (N-q)%N\right)}$

你提到的“黄色框缺失系数”，本质就是这些对称推导出来的冗余项，逆变换能完整还原数据也验证了所有信息都被保留在R2HC输出中。

3D混合变换的适配（2D R2HC + 1D DCT） #

回到你的3D PDE场景，用fftw_plan_r2r_3d实现时，需要指定两个维度为FFTW_R2HC（对应DFT/周期性边界），第三个维度为FFTW_REDFT01（对应DCT-II，是诺依曼边界PDE求解的常用选择，其逆变换为DCT-III）。

系数处理步骤：

1. 对两个R2HC维度，按照上面的2D规则重构完整的复数DFT系数。
2. 第三个维度的DCT系数是实值的，直接从输出数组中提取即可，无需共轭处理（DCT本身是实变换，输出全为实数）。
3. 在傅里叶空间中应用PDE的对角化算子（比如拉普拉斯算子对应$-k_x^2 -k_y^2 -k_z^2$，其中$k_z$是DCT的波数），再执行逆变换（FFTW_HC2R + FFTW_REDFT10）回到实空间。

额外建议 #

FFTW手册不推荐多次使用R2HC的核心原因是多维R2HC的存储格式极易混淆，不如拆分处理直观：

• 先用fftw_plan_dft_r2c_3d处理两个周期性维度，得到半复数格式的2D DFT结果。
• 再对第三个维度单独用fftw_plan_r2r_1d执行DCT变换。
  这种拆分方式能降低格式理解的复杂度，减少出错概率。

内容的提问来源于stack exchange，提问作者Dennis