-2

我对这个程序有一个问题：

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

我想转换一个信号，然后逆向返回，但是前半部分的输出是错误的。

你能帮我找出错误吗？

非常感谢！

4

1 回答 1

2

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

cufftDoubleComplex不一样cufftComplex。使用时cufftDoubleComplex，您的变换类型应为 Z2Z，而不是 C2C。

此外，为了在使用 CUFFT 进行正向变换和逆变换时查看数据奇偶性，有必要将结果除以信号大小：

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;
}
于 2014-10-05T19:29:36.297 回答