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

Question

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

Answer 1

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

如果您对所有 '前者缺失的后者的（最后一个维度的）切片。

如果 Nz 是偶数，则表示对除第一个和最后一个切片之外的所有切片进行双重计数。
如果 Nz 是奇数，重复计算除第一个之外的所有切片。

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

这是使用 Numpy/Python 对该想法进行的实验验证，我相信其实数到复数 FFT 的符号与 FFTW 相匹配：通过 Ny = 20 通过 [=15 生成 Nx = 10 =] 真正的数组。计算其复数到复数的 3D FFT（通过 Nz 复数组产生 Nx 乘以 Ny）和其实数到复数的 3D FFT（通过 Ny 通过（Nz/2+1) 复杂数组）。如果您对除第一个和最后一个切片以外的所有部分进行双重计数，请验证前者的对数幅度之和与后者的对数幅度之和相同 , 因为 Nz 是偶数。

代码：

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 区间。

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

tl;dr 对实数到复数输出的对数幅值求和。加倍。从该结果中减去第一个切片（在第 3 维中）的对数和幅值。如果 Nz 是奇数，你就完成了。如果 Nz 是偶数，还减去最后一个切片的对数和幅值。

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

Adding up all the values for each frequency from FFTW

fft

fftw