2012-04-29 46 views
0

我试图做一个3D FFT与FFTW库,但我有逆变换一些困难。3D C2C FFT与FFTW库

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

虽然我的数据是我使用的复杂复杂的转换,如想 真实数据后由OpenCL的FFT替换它仅支持复杂:

起初我通过做有序转变到复杂的转换。

在三维傅立叶空间我做一个很简单的低通滤波器:

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-测试用例查看此对称的傅立叶空间:据我看到

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 

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]; 
    } 
} 

这给了我下面的输出它没有对称值。为什么?

这将是很好,如果有人可以给我一个提示。

问候

回答

0

如果你想逆FFT后严格的实际结果(减去使用有限尺寸算术通常的数字噪声),你必须确保你喂一个完整的IFFT输入的数据完全是共轭对称(去年一半矢量是前半部分的镜像复共轭)。它似乎没有强迫你的数据是这样的。

+0

我已更新我的代码(请参阅上文)以保持对称性,但我没有按预期工作。 – DerHandwerk

+0

我还添加了一些代码来打印复杂的光谱。 – DerHandwerk

1

这种联系是一种误导。真实信号的DFT通常不会导致输出采样的一半为零。它只是施加(共轭)对称。

在筛选器代码

所以,你只操纵的输出值,你应该的一半。每当您操纵输出库n时,您还需要操作库Nn(其中N是DFT的长度),以便保持对称性,以便在应用时为您提供真实的结果逆DFT。

我的建议是先解决一个更简单的问题 - 一维滤波器。一旦你有正确的,那么它应该很容易扩展到3D。

+0

这听起来很合理。这是否意味着在1D FFT的情况下,我只需要从n = 0到N/2,因为在每一步中我都会过滤complex [n]和complex [N-n]? – DerHandwerk

+0

@DerHandwerk:不知道我明白你在问什么。长度为N的1D DFT将只有N/2 + 1个独立输出样本;关系是y [n] = y [N-n] *(其中“*”表示[复共轭](http://en.wikipedia.org/wiki/Complex_conjugate))。 –

+0

好吧,就我所见的一维FFT而言,我的低通必须执行以下操作:'如果振幅>值,则y [n] = 0; y [N - n - 1] = 0; end' for n = 0 to N - 1 – DerHandwerk