带有 fftw 库的 3d c2c fft

Question

我正在尝试使用 FFTW 库进行 3D FFT，但我在逆变换方面遇到了一些困难。

首先我通过以下方式进行前言转换：

fftwf_plan_dft_3d(_dimensions[0], _dimensions[1], _dimensions[2], (fftwf_complex*)_inputBuffer, (fftwf_complex*)_outputBuffer, FFTW_FORWARD, FFTW_ESTIMATE);

虽然我的数据是真实数据，但我正在使用复杂到复杂的转换，因为我想 稍后将其替换为仅支持复杂到复杂转换的 opencl fft。

在 3D 傅立叶空间中，我做了一个非常简单的低通滤波器：

for all x, y, z:

// global position of the current bin
int gid = (y * w + x) + (z * w * h);

// position of the symmetric bin
vec3 conPos(M - x - 1, N - y - 1, L - z - 1);

// global position of the symmetric element
int conGid = (conPos.y * w + conPos.x) + (conPos.z * w * h);

if (sqrt(x * x + y * y + z * z) > 500)
{
    complex[gid].real = 0.0f;
    complex[gid].imag = 0.0f;
    complex[conGid].real = 0.0f;
    complex[conGid].imag = 0.0f;
}

最后是逆变换：

fftwf_plan_dft_3d(_dimensions[0], _dimensions[1], _dimensions[2], (fftwf_complex*)_inputBuffer, (fftwf_complex*)_outputBuffer, FFTW_BACKWARD, FFTW_ESTIMATE);
// normalization ...

结果并不如我所料。在逆变换之后，虚部并不像他们应该的那样全为零。

据我所知，在对真实数据进行正向变换后，只使用了总缓冲区大小的一半，而另一半没有共轭复数值。 （请参阅：c2c with real data）如果是这种情况，我必须在反向转换之前自行计算它们，但我在 fftw 文档中找不到提示，哪一半是计算的，哪一半不是。

我编写了一个非常简单的 2D-Test-Case 来查看傅立叶空间中的这种对称性：

int w = 4;
int h = 4;
int size  = w * h;

cl_float rawImage[16] = ...; // loading image

fftwf_complex *complexImage = (fftwf_complex*) fftwf_malloc(sizeof(fftwf_complex) * size);
fftwf_complex *freqBuffer = (fftwf_complex*) fftwf_malloc(sizeof(fftwf_complex) * size);

for (int i = 0; i < size; i++)
{
    complexImage[i][0] = rawImage[i]; complexImage[i][1] = 0.0f;
}

fftwf_plan forward = fftwf_plan_dft_2d(w, h, complexImage, freqBuffer, FFTW_FORWARD, FFTW_ESTIMATE);

fftwf_execute(forward);

for (int y = 0; y < h; y++)
{
    for (int x = 0; x < w; x++)
    {
        int gid = y * w + x;
        qDebug() << gid << "real:" << freqBuffer[gid][0] << "imag:" << freqBuffer[gid][1];
    }
}

这给了我以下输出：

gid
0    real 3060 imag 0 
1    real 510 imag 510 
2    real 0 imag 0 
3    real 510 imag -510 
4    real 510 imag 510 
5    real 0 imag -510 
6    real 0 imag 0 
7    real -510 imag 0 
8    real 0 imag 0 
9    real 0 imag 0 
10   real 0 imag 0 
11   real 0 imag 0 
12   real 510 imag -510 
13   real -510 imag 0 
14   real 0 imag 0 
15   real 0 imag 510 

据我所知，没有对称值。为什么？

如果有人能给我一个提示，那就太好了。

问候

狼

Answer 1

如果您想在逆 FFT 后得到严格的实数结果（减去使用有限大小算术的常见数值噪声），您必须确保提供完整 IFFT 的输入数据是完全共轭对称的（后半部分）向量是前半部分的镜像复共轭）。您似乎没有强迫您的数据如此。

Answer 2

该链接具有误导性。仅实数信号的 DFT（通常）不会导致一半的输出样本为零。它只是对它们施加（共轭）对称性。

因此，在您的过滤器代码中，您只操作了应有的输出值的一半。每次操作输出 bin n 时，还需要操作 bin N-n（其中 N 是 DFT 的长度），以便保持对称性，当您应用逆 DFT 时，这将为您提供仅实数的结果。

我的建议是首先解决一个更简单的问题——一维过滤器。一旦你得到了正确的结果，那么它应该很容易扩展到 3D。