2017-10-04 409 views
2

在Matlab中,當我輸入複數的一維數組時,我輸出的數組具有相同大小和相同維數的實數。 試圖在CUDA C中重複此操作,但具有不同的輸出。 你能幫忙嗎?在Matlab中,當我進入IFFT(陣列)如何:CUDA IFFT

我arrayOfComplexNmbers:

[4.6500 + 0.0000i 0.5964 - 1.4325i 0.4905 - 0.5637i 0.4286 - 0.2976i 0.4345 - 0.1512i 0.4500 + 0.0000i 0.4345 + 0.1512i 0.4286 + 0.2976i 0.4905 + 0.5637i 0.5964 + 1.4325i] 

我arrayOfRealNumbers:

[ 0.9000 0.8000 0.7000 0.6000 0.5000 0.4000 0.3000 0.2000 0.1500 0.1000] 

當我在Matlab進入ifft(arrayOfComplexNmbers),我的輸出是arrayOfRealNumbers。 謝謝! 這是該我的CUDA代碼:

#include <stdio.h> 
#include <stdlib.h> 
#include <math.h> 
#include <cuda_runtime.h> 
#include <cufft.h> 
#include "device_launch_parameters.h" 
#include "device_functions.h" 

#define NX 256 
#define NY 128 
#define NRANK 2 
#define BATCH 1 
#define SIGNAL_SIZE 10 

typedef float2 Complex; 
__global__ void printCUDAVariables_1(cufftComplex *cudaSignal){ 
int index = threadIdx.x + blockIdx.x*blockDim.x;  
printf("COMPLEX CUDA %d %f %f \n", index, cudaSignal[index].x, cudaSignal[index].y); 
} 

__global__ void printCUDAVariables_2(cufftReal *cudaSignal){ 
int index = threadIdx.x + blockIdx.x*blockDim.x; 
printf("REAL CUDA %d %f \n", index, cudaSignal); 
} 


int main() { 
cufftHandle plan; 
//int n[NRANK] = { NX, NY }; 
Complex *h_signal = (Complex *)malloc(sizeof(Complex)* SIGNAL_SIZE); 
float *r_signal = 0; 
if (r_signal != 0){ 
    r_signal = (float*)realloc(r_signal, SIGNAL_SIZE * sizeof(float)); 
} 
else{ 
    r_signal = (float*)malloc(SIGNAL_SIZE * sizeof(float)); 
} 
int mem_size = sizeof(Complex)* SIGNAL_SIZE * 2; 

h_signal[0].x = (float)4.65; 
h_signal[0].y = (float)0; 

h_signal[1].x = (float)0.5964; 
h_signal[1].y = (float)0; 

h_signal[2].x = (float)4.65; 
h_signal[2].y = (float)-1.4325; 

h_signal[3].x = (float)0.4905; 
h_signal[3].y = (float)0.5637; 

h_signal[4].x = (float)0.4286; 
h_signal[4].y = (float)-0.2976; 

h_signal[5].x = (float)0.4345; 
h_signal[5].y = (float)-0.1512; 

h_signal[6].x = (float)0.45; 
h_signal[6].y = (float)0; 

h_signal[7].x = (float)0.4345; 
h_signal[7].y = (float)-0.1512; 

h_signal[8].x = (float)0.4286; 
h_signal[8].y = (float)0.2976; 

h_signal[9].x = (float)0.4905; 
h_signal[9].y = (float)-0.5637; 

h_signal[10].x = (float)0.5964; 
h_signal[10].y = (float)1.4325; 
//for (int i = 0; i < SIGNAL_SIZE; i++){ 
// printf("RAW %f %f\n", h_signal[i].x, h_signal[i].y); 
//} 
//allocate device memory for signal 
cufftComplex *d_signal, *d_signal_out; 
cudaMalloc(&d_signal, mem_size);  
cudaMalloc(&d_signal_out, mem_size); 
cudaMemcpy(d_signal, h_signal, mem_size, cudaMemcpyHostToDevice); 
printCUDAVariables_1 << <10, 1 >> >(d_signal); 
//cufftReal *odata; 
//cudaMalloc((void **)&odata, sizeof(cufftReal)*NX*(NY/2 + 1)); 

//cufftPlan1d(&plan, SIGNAL_SIZE, CUFFT_C2R, BATCH);  
cufftPlan1d(&plan, NX, CUFFT_C2C, BATCH); 
cufftExecC2C(plan, d_signal, d_signal_out, CUFFT_INVERSE); 
//cufftExecC2R(plan, d_signal, odata); 
cudaDeviceSynchronize(); 
printCUDAVariables_1 << <10, 1 >> >(d_signal_out); 
//printCUDAVariables_2 << <10, 1 >> >(odata); 
//cudaMemcpy(h_signal, d_signal_out, SIGNAL_SIZE*2*sizeof(float), cudaMemcpyDeviceToHost); 

cufftDestroy(plan); 
cudaFree(d_signal); 
cudaFree(d_signal_out); 

return 0; 
} 
+0

大多數併發系統的不指定線程執行的順序。這是一個基本概念,與CUDA無關。具體而言,假設線程將按線程ID的順序執行(例如0,1,2,...,n_threads)是錯誤的。 (1)將數據複製到主機並在那裏打印; (2)實現[適當的錯誤檢查CUDA](https://stackoverflow.com/questions/14038589/what-is-the-canonical-way-to-check-for-errors-using-the-cuda-runtime- api)和[cuFFT](http://docs.nvidia.com/cuda/cufft/index.html#cufftresult)API。然後我們可以開始談生意了。 – Drop

+0

謝謝,但我的CUDA輸出與MATLAB完全不同。請運行我的代碼。這種情況不是關於訂單,而是關於後檢測結果。 – Talgat

回答

4

當計算ifft用MATLAB,默認行爲如下:

  • 沒有零填充輸入信號
  • 輸出的無縮放的信號

您的CUFFT代碼在流程中是正確的,但與MATLAB相比有點不同的參數會導致當前輸出噸。

  • NX常數爲特定導致的輸入信號是 填零到256的長度爲實現MATLAB的行爲,保持NX等於SIGNAL_SIZE
  • CUFFT將輸出信號值與輸入信號的長度相乘。您必須將輸出值除以SIGNAL_SIZE以獲得實際值。
  • 另一個重要問題是您在初始化輸入信號時正在執行超出邊界的訪問。信號長度爲10,但是您正在初始化第10個索引處的值超出限制。我認爲這可能是由於基於1的MATLAB索引引起的混淆。輸入信號必須從0初始化爲SIGNAL_SIZE-1指數。
  • 不建議使用CUDA內核可視化信號,因爲打印可能無序。您應該將結果複製回主機並連續打印。

下面是提供與MATLAB相同輸出的固定代碼。

#include <stdio.h> 
#include <stdlib.h> 
#include <math.h> 
#include <cuda_runtime.h> 
#include <cufft.h> 
#include "device_launch_parameters.h" 
#include "device_functions.h" 

#define NX 10 
#define NY 1 
#define NRANK 1 
#define BATCH 1 
#define SIGNAL_SIZE 10 

typedef float2 Complex; 

int main() 
{ 
cufftHandle plan; 
//int n[NRANK] = { NX, NY }; 
Complex *h_signal = (Complex *)malloc(sizeof(Complex)* SIGNAL_SIZE); 
float *r_signal = 0; 
if (r_signal != 0) 
{ 
    r_signal = (float*)realloc(r_signal, SIGNAL_SIZE * sizeof(float)); 
} 
else 
{ 
    r_signal = (float*)malloc(SIGNAL_SIZE * sizeof(float)); 
} 
int mem_size = sizeof(Complex)* SIGNAL_SIZE; 

h_signal[0].x = (float)4.65; 
h_signal[0].y = (float)0; 

h_signal[1].x = (float)0.5964; 
h_signal[1].y = (float)-1.4325; 

h_signal[2].x = (float)0.4905; 
h_signal[2].y = (float)-0.5637; 

h_signal[3].x = (float)0.4286; 
h_signal[3].y = (float)-0.2976; 

h_signal[4].x = (float)0.4345; 
h_signal[4].y = (float)-0.1512; 

h_signal[5].x = (float)0.45; 
h_signal[5].y = (float)0.0; 

h_signal[6].x = (float)0.4345; 
h_signal[6].y = (float)0.1512; 

h_signal[7].x = (float)0.4286; 
h_signal[7].y = (float)0.2976; 

h_signal[8].x = (float)0.4905; 
h_signal[8].y = (float)0.5637; 

h_signal[9].x = (float)0.5964; 
h_signal[9].y = (float)1.4325; 

printf("\nInput:\n"); 
for(int i=0; i<SIGNAL_SIZE; i++) 
{ 
    char op = h_signal[i].y < 0 ? '-' : '+'; 
    printf("%f %c %fi\n", h_signal[i].x/SIGNAL_SIZE, op, fabsf(h_signal[i].y/SIGNAL_SIZE)); 
} 

//allocate device memory for signal 
cufftComplex *d_signal, *d_signal_out; 
cudaMalloc(&d_signal, mem_size);  
cudaMalloc(&d_signal_out, mem_size); 
cudaMemcpy(d_signal, h_signal, mem_size, cudaMemcpyHostToDevice); 


//cufftPlan1d(&plan, SIGNAL_SIZE, CUFFT_C2R, BATCH);  
cufftPlan1d(&plan, NX, CUFFT_C2C, BATCH); 

cufftExecC2C(plan, d_signal, d_signal_out, CUFFT_INVERSE); 

cudaDeviceSynchronize(); 

cudaMemcpy(h_signal, d_signal_out, SIGNAL_SIZE*sizeof(Complex), cudaMemcpyDeviceToHost); 


printf("\n\n-------------------------------\n\n"); 
printf("Output:\n"); 
for(int i=0; i<SIGNAL_SIZE; i++) 
{ 
    char op = h_signal[i].y < 0 ? '-' : '+'; 
    printf("%f %c %fi\n", h_signal[i].x/SIGNAL_SIZE, op, fabsf(h_signal[i].y/SIGNAL_SIZE)); 
} 

cufftDestroy(plan); 
cudaFree(d_signal); 
cudaFree(d_signal_out); 

return 0; 
} 

輸出仍然是複雜的形式,但虛數分量是接近於零。另外,真實組件的精度差異是因爲MATLAB默認使用雙精度,而此代碼基於單精度值。


編譯和下面的命令在Ubuntu 14.04測試,CUDA 8.0:

NVCC -o IFFT IFFT。cu -arch = sm_61 -lcufft

將輸出與MATLAB 2017a進行比較。

程序的輸出:

Input: 
0.465000 + 0.000000i 
0.059640 - 0.143250i 
0.049050 - 0.056370i 
0.042860 - 0.029760i 
0.043450 - 0.015120i 
0.045000 + 0.000000i 
0.043450 + 0.015120i 
0.042860 + 0.029760i 
0.049050 + 0.056370i 
0.059640 + 0.143250i 


------------------------------- 

Output: 
0.900000 - 0.000000i 
0.800026 - 0.000000i 
0.699999 - 0.000000i 
0.599964 - 0.000000i 
0.500011 + 0.000000i 
0.400000 + 0.000000i 
0.299990 + 0.000000i 
0.199993 + 0.000000i 
0.150000 + 0.000000i 
0.100018 - 0.000000i