将 FFTW 中每个频率的所有值相加

Question

我已经使用 FFTW (fftw_plan_dft_r2c_3d) 运行了一个 3D 傅里叶变换，我想总结每个频率的变换值（对数），包括实际未存储在输出数组（我知道大小是Nx x Ny x (Nz/2 + 1)）。如何在不重复计算的情况下做到这一点？

Answer 1

很好的问题。抱歉我的回答有点啰嗦，我想确保我没有犯任何错误。来了——

如果您对所有“切片”进行双重计算，复数到复数 3D FFT 的对数幅度之和将等于实数到复数 3D FFT 的对数幅度之和（最后一个维度的）后者缺少前者。

• 如果Nz 是偶数，则意味着重复计算除第一个和最后一个切片之外的所有切片。
• 如果 Nz 是奇数，则重复计算除第一个以外的所有切片。

（这是因为偶数长度的实数到复数 DFT 包括 -π 弧度角频率（对应于 -1 的相量），而奇数长度的 DFT 则缺少它。我不记得这种模式, 所以我总是在单位圆周围画出 N=4 vs N=3 相量来提醒自己奇数或偶数是否包含-π弧度。）

这是使用 Numpy/Python 对该想法的实验验证，我相信它的实数到复数 FFT 表示法与 FFTW 相匹配：通过Ny = 20 由Nz = 8 实数数组生成一个Nx = 10。计算其复数到复数 3D FFT（产生 Nx 乘 Ny 乘 Nz 复数数组）和其实数到复数 3D FFT（产生 Nx 乘 Ny 乘 (Nz/2+1) 复数）大批）。如果您重复计算除第一个和最后一个切片之外的所有切片，请验证前者的对数大小之和与后者的对数大小之和 相同，因为 @ 987654328@是偶数。

代码：

import numpy as np
import numpy.fft as fft

Nx = 10
Ny = 20
Nz = 8

x = np.random.randn(Nx, Ny, Nz)

Xf = fft.fftn(x)
Xfr = fft.rfftn(x)

energyProduct1 = np.log10(np.abs(Xf)).sum()

lastSlice = -1 if Nz % 2 is 0 else None
energyProduct2 = np.log10(np.abs(np.dstack((Xfr, Xfr[:, :, 1:lastSlice])))).sum()

print('Difference: %g' % (energyProduct1 - energyProduct2))
# Difference: -4.54747e-13

如果你用奇数Nz 重新运行它，你会看到复数到复数和实数到复数之间的差异保持在 0 的机器精度内。

np.dstack((Xfr, Xfr[:, :, 1:lastSlice))（dstack、fft.rfftn 的文档）将rfftn 输出与第三维中的倒数第二个切片堆叠在一起——倒数第二个，因为Nz 是偶数，而您不想这样做重复计算 0 或 -π DFT bin。

当然，另一种方法是计算实数到复数数组的对数幅度之和，加倍，然后减第一个切片和（如果Nz 是偶数）最后一个切片的对数幅度之和。

tl;dr 将实数到复数输出的对数幅度相加。加倍。从此结果中减去第一个切片（在第 3 维中）的总和对数幅度。如果Nz 是奇怪的，你就完成了。如果Nz 是偶数，还减去最后一个切片的 sum-log-magnitudes。