2012-04-29 105 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