CUFFT输出不正确

Question

我有这个程序有问题：

#include <stdlib.h> 
#include <stdio.h> 
#include <string.h> 
#include <math.h> 
#include <cufft.h> 
#include <cuComplex.h> 

#define SIGNAL_SIZE  1024 

int main(int argc, char **argv) { 
    cudaEvent_t start, stop; 
    cudaEventCreate(&start); 
    cudaEventCreate(&stop); 

    // Allocate host memory for the signal 
    cuDoubleComplex *h_signal = (cuDoubleComplex *) malloc(sizeof(cuDoubleComplex) * SIGNAL_SIZE); 

    // Initalize the memory for the signal 
    for (unsigned int i = 0; i < SIGNAL_SIZE; ++i) { 
     if((double)i/SIGNAL_SIZE>=0 && (double)i/SIGNAL_SIZE<0.5) h_signal[i].x = (double)i/SIGNAL_SIZE; 
     else if((double)i/SIGNAL_SIZE>=0.5 && (double)i/SIGNAL_SIZE<1) h_signal[i].x = (double)i/SIGNAL_SIZE-1; 
     h_signal[i].y = 0; 
    } 

// Allocate device memory for signal 
    cuDoubleComplex *d_signal; 

    cudaMalloc((void **) &d_signal, SIGNAL_SIZE*sizeof(cuDoubleComplex)); 
    // Copy host memory to device 
    cudaMemcpy(d_signal, h_signal, SIGNAL_SIZE*sizeof(cuDoubleComplex), cudaMemcpyHostToDevice); 


cudaEventRecord(start, 0); 
    cufftHandle plan; 
    cufftPlan1d(&plan, SIGNAL_SIZE , CUFFT_C2C, 1); 

    // FFT computation 
    cufftExecC2C(plan, (cufftComplex *) d_signal, (cufftComplex *) d_signal, 
     CUFFT_FORWARD); 

    cufftExecC2C(plan, (cufftComplex *) d_signal, (cufftComplex *) d_signal, CUFFT_INVERSE); 

    cuDoubleComplex *h_signal_inv =(cuDoubleComplex *) malloc(sizeof(cuDoubleComplex) * SIGNAL_SIZE); 
    cudaMemcpy(h_signal_inv, d_signal, sizeof(cuDoubleComplex) * SIGNAL_SIZE, cudaMemcpyDeviceToHost); 


    cudaEventRecord(stop, 0); 
    cudaEventSynchronize(stop); 

    float elapsedTime; 
    cudaEventElapsedTime(&elapsedTime, start, stop); 
    printf("Elapsed Time: %3.1f ms\n", elapsedTime); 


    for(int i=0;i<SIGNAL_SIZE;i++) printf("\n%f %f", h_signal[i].x, h_signal_inv[i].x); 

    cufftDestroy(plan); 

    free(h_signal); 
    free(h_signal_inv); 

    cudaFree(d_signal); 

    cudaDeviceReset(); 
    return 0; 
} 

我想变换的信号，然后再回来与逆，但输出是错误的上半年。

你能帮我找到错误吗？

非常感谢！

Answer 1

你正在让你的数据类型感到困惑。

cufftDoubleComplex与cufftComplex不一样。当使用cufftDoubleComplex时，your transform type should be Z2Z, not C2C。

另外，为了看到数据奇偶性做时通过逆一个正向变换，然后使用变换CUFFT，它necessary to divide the result by the signal size：

CUFFT执行未归一化的FFT;也就是说，对输入数据组执行正向FFT，然后对所得到的组进行逆FFT，得到等于输入的数据，按照元素的数量缩放。通过数据集尺寸的倒数来缩放变换，留给用户以适合的方式执行。

下面的代码具有上述问题解决，应该给你更好的结果：

#include <stdlib.h> 
#include <stdio.h> 
#include <string.h> 
#include <math.h> 
#include <cufft.h> 
#include <cuComplex.h> 

#define SIGNAL_SIZE  1024 

int main(int argc, char **argv) { 
    cudaEvent_t start, stop; 
    cudaEventCreate(&start); 
    cudaEventCreate(&stop); 

    // Allocate host memory for the signal 
    cuDoubleComplex *h_signal = (cuDoubleComplex *) malloc(sizeof(cuDoubleComplex) * SIGNAL_SIZE); 

    // Initalize the memory for the signal 
    for (unsigned int i = 0; i < SIGNAL_SIZE; ++i) { 
     if((double)i/SIGNAL_SIZE>=0 && (double)i/SIGNAL_SIZE<0.5) h_signal[i].x = (double)i/SIGNAL_SIZE; 
     else if((double)i/SIGNAL_SIZE>=0.5 && (double)i/SIGNAL_SIZE<1) h_signal[i].x = (double)i/SIGNAL_SIZE-1; 
     h_signal[i].y = 0; 
    } 

// Allocate device memory for signal 
    cuDoubleComplex *d_signal; 

    cudaMalloc((void **) &d_signal, SIGNAL_SIZE*sizeof(cuDoubleComplex)); 
    // Copy host memory to device 
    cudaMemcpy(d_signal, h_signal, SIGNAL_SIZE*sizeof(cuDoubleComplex), cudaMemcpyHostToDevice); 


cudaEventRecord(start, 0); 
    cufftHandle plan; 
    cufftPlan1d(&plan, SIGNAL_SIZE , CUFFT_Z2Z, 1); 

    // FFT computation 
    cufftExecZ2Z(plan, d_signal, d_signal, CUFFT_FORWARD); 

    cufftExecZ2Z(plan, d_signal, d_signal, CUFFT_INVERSE); 

    cuDoubleComplex *h_signal_inv =(cuDoubleComplex *) malloc(sizeof(cuDoubleComplex) * SIGNAL_SIZE); 
    cudaMemcpy(h_signal_inv, d_signal, sizeof(cuDoubleComplex) * SIGNAL_SIZE, cudaMemcpyDeviceToHost); 


    cudaEventRecord(stop, 0); 
    cudaEventSynchronize(stop); 

    float elapsedTime; 
    cudaEventElapsedTime(&elapsedTime, start, stop); 
    printf("Elapsed Time: %3.1f ms\n", elapsedTime); 


    for(int i=0;i<SIGNAL_SIZE;i++) printf("\n%f %f", h_signal[i].x, h_signal_inv[i].x/SIGNAL_SIZE); 

    cufftDestroy(plan); 

    free(h_signal); 
    free(h_signal_inv); 

    cudaFree(d_signal); 

    cudaDeviceReset(); 
    return 0; 
}