→ loveme00835:顏色好豐富@@ 02/01 04:11
推 aada:謝謝大大, 假設我今天要對FFT後的數據處理的話,我是不是要再 02/01 11:12
→ aada:另外創造一個4x4的矩陣,將FFT後的值傳入到新的矩陣呢 02/01 11:13
→ lgen7604:其實結果用4x3的矩陣就存得下了 只是處理起來會比較麻煩 02/01 14:23
→ lgen7604:如果想要方便處理的話 就轉換到4x4的矩陣吧 02/01 14:24
其實你執行的結果並沒有錯
如果你仔細看document的話
會注意到在cufftExecR2C()的說明中有提到
This function stores the non-redundant Fourier coefficients in the odata array.
這句話指出執行cufftExecR2C()的結果是non-redundant
這就牽扯到實數DFT的特性問題
當DFT的input為實數數列的時候
第k項和第n-k項必定為共軛複數
假設n=8 共軛情形如下
a[0]a[1]a[2]a[3]a[4]a[5]a[6]a[7]
如果n=9 共軛情形如下
a[0]a[1]a[2]a[3]a[4]a[5]a[6]a[7]a[8]
所以實際上DFT的結果只需要記錄n/2+1項就可以了
你測試的範例為n=4
所以每次FFT的結果只要記錄3個複數就可以了
正確的範例我改成如下
#include <stdio.h>
#include <stdlib.h>
#include <cuda_runtime.h>
#include <cufft.h>
#include <math.h>
#define H 4
#define W 4
#define T H*W
int main(){
cufftReal *a;
a = (cufftReal*) malloc( sizeof(cufftReal)*(T) );
for(int i=0; i<T; i++)
a[i] = 2.0;
cufftHandle plan;
cufftComplex *odata;
cufftReal *idata;
cudaMalloc( (void**)&idata, sizeof(cufftReal)*(T) );
cudaMalloc( (void**)&odata, sizeof(cufftComplex)*(H*(W/2+1)) );
cudaMemcpy(idata, a, sizeof(cufftReal) * T,cudaMemcpyHostToDevice);
cufftPlan1d(&plan, W, CUFFT_R2C, H);
cufftExecR2C(plan, (cufftReal*)idata, odata);
cufftComplex *FFT_odata;
FFT_odata = (cufftComplex*) malloc( sizeof(cufftComplex)*H*(W/2+1) );
cudaMemcpy( FFT_odata, odata, sizeof(cufftComplex)*H*(W/2+1),
cudaMemcpyDeviceToHost );
for(int i=0; i <H; i++) {
for(int j=0; j <(W/2+1); j++) {
printf("(%d,%d)=%.3f%+.3fi ", i, j, FFT_odata[i*(W/2+1)+j].x,
FFT_odata[i*(W/2+1)+j].y );
}
for(int j=(W-1)/2; j >= 1; j--) {
printf("(%d,%d)=%.3f%+.3fi ", i, W-j, FFT_odata[i*(W/2+1)+j].x,
-FFT_odata[i*(W/2+1)+j].y );
}
printf("\n");
}
cufftDestroy(plan);
cudaFree(idata);
cudaFree(odata);
system("PAUSE");
return 0;
}
最後要注意的是MATLAB對2D的array做FFT的時候是對column執行
上述範例則是對row執行FFT
--
※ 發信站: 批踢踢實業坊(ptt.cc)
◆ From: 122.120.32.11
※ 編輯: lgen7604 來自: 122.120.32.11 (01/31 20:20)