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

Question

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

Answer 1

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

如果对所有“切片”进行双重计数，则复数到复数 3D FFT 的对数幅值总和将等于实数到复数 3D FFT 的对数幅值总和(最后一个维度)前者缺少后者。

• 如果 Nz 为偶数，则意味着对除第一个和最后一个切片之外的所有切片进行双重计数。
• 如果 Nz 为奇数，则对除第一个之外的所有切片进行双倍计数。

(这是因为偶数长度的实数到复数 DFT 包括 -π 弧度角频率(对应于 -1 的相量)，而奇数长度的 DFT 则没有达到它。我从来不记得这种模式，所以我总是在单位圆周围绘制 N=4 与 N=3 相量，以提醒自己奇数或偶数是否包含 -π rad。)

这是使用 Numpy/Python 对这个想法进行的实验验证，我相信其实数到复数 FFT 的表示法与 FFTW 的表示法相匹配:通过 Ny = 20 生成 Nx = 10 code> by Nz = 8 实数数组。计算其复数到复数 3D FFT(产生 Nx × Ny × Nz 复数数组)和其实数到复数 3D FFT(产生 Nx × 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 bin。

当然，另一种方法是计算实数到复数数组的对数幅值总和，将其加倍，然后减去第一个切片和最后一个切片(如果 Nz 为偶数)的对数幅值之和。

tl;dr 对实数到复数输出的对数幅度求和。加倍。从此结果中减去第一个切片(在第三维中)的和对数幅度。如果 Nz 是奇数，那么你就完成了。如果 Nz 是偶数，还要减去最后一个切片的和对数幅值。