-2

我是新来的。我正在使用 FFT，我需要编写一个简单的代码，但它不起作用。我需要用 cufft a sin(x) 进行变换并返回，但在变换之间，我需要将结果乘以 2，这样，当我使用逆变换返回结果时，我将得到 2*sin( x) 例如。使用 fftw.h，我只需将我的 d_signal[i] 乘以 2，当我返回时，我有 2*sin(x)，但我曾经使用 complex.h。任何想法？谢谢。

# define SIGNAL_SIZE 64
# define PI acos(-1.0)
# define x 2*PI/SIGNAL_SIZE
void runTest(int argc, char **argv)
{
    printf("[simpleCUFFT] is starting...\n");

    findCudaDevice(argc, (const char **)argv);

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

    cufftComplex *h_reversed_signal = (Complex *)malloc(sizeof(Complex) * SIGNAL_SIZE);

// Initalize the memory for the signal
for (unsigned int i = 0; i < SIGNAL_SIZE; ++i)
{
    h_signal[i].x = sin(i*x);
    h_signal[i].y = 0;
}
cufftComplex *d_signal;
checkCudaErrors(cudaMalloc((void **)&d_signal, SIGNAL_SIZE*sizeof(cufftComplex)));
// Copy host memory to device
    checkCudaErrors(cudaMemcpy(d_signal, h_signal, SIGNAL_SIZE*sizeof(cufftComplex),
                           cudaMemcpyHostToDevice));
cufftHandle plan;
checkCudaErrors(cufftPlan1d(&plan, SIGNAL_SIZE, CUFFT_C2C, 1));

// Transform signal and kernel
printf("Transforming signal cufftExecC2C\n");
checkCudaErrors(cufftExecC2C(plan, (cufftComplex *)d_signal, (cufftComplex *)d_signal, CUFFT_FORWARD));
getLastCudaError("Kernel execution failed [ ComplexPointwiseMulAndScale ]");
for (unsigned int i = 0; i < SIGNAL_SIZE; ++i)
{
    d_signal[i].x = 2*d_signal[i].x;
    d_signal[i].y = 2*d_signal[i].y;
}
// Transform signal back
printf("Transforming signal back cufftExecC2C\n");
checkCudaErrors(cufftExecC2C(plan, (cufftComplex *)d_signal, (cufftComplex *)d_signal, CUFFT_INVERSE));

   // Copy device memory to host
checkCudaErrors(cudaMemcpy(h_reversed_signal, d_signal, SIGNAL_SIZE*sizeof(cufftComplex),
                           cudaMemcpyDeviceToHost));
// check result
    for (unsigned int i = 0; i < SIGNAL_SIZE; ++i)
    {
            h_reversed_signal[i].x = h_reversed_signal[i].x / (float)SIGNAL_SIZE;
            h_reversed_signal[i].y = h_reversed_signal[i].y/(float)SIGNAL_SIZE;
            printf("first : %f %f  after %f %f \n", h_signal[i].x, h_signal[i].y, h_reversed_signal[i].x, h_reversed_signal[i].y);
    }
bool bTestResult = sdkCompareL2fe((float *)h_reversed_signal, (float *)h_signal, 2 * SIGNAL_SIZE, 1e-5f);
//Destroy CUFFT context
checkCudaErrors(cufftDestroy(plan));
// cleanup memory
free(h_signal);
free(h_reversed_signal);
checkCudaErrors(cudaFree(d_signal));
cudaDeviceReset();
}

// Pad data
    int PadData(const Complex *signal, Complex **padded_signal, int signal_size,
        const Complex *filter_kernel, Complex **padded_filter_kernel, int filter_kernel_size)
{
int minRadius = filter_kernel_size / 2;
int maxRadius = filter_kernel_size - minRadius;
int new_size = signal_size + maxRadius;

// Pad signal
Complex *new_data = (Complex *)malloc(sizeof(Complex) * new_size);
memcpy(new_data +           0, signal,              signal_size * sizeof(Complex));
memset(new_data + signal_size,      0, (new_size - signal_size) * sizeof(Complex));
*padded_signal = new_data;

// Pad filter
new_data = (Complex *)malloc(sizeof(Complex) * new_size);
memcpy(new_data +                    0, filter_kernel + minRadius,                       maxRadius * sizeof(Complex));
memset(new_data +            maxRadius,                         0, (new_size - filter_kernel_size) * sizeof(Complex));
memcpy(new_data + new_size - minRadius,             filter_kernel,                       minRadius * sizeof(Complex));
*padded_filter_kernel = new_data;

return new_size;

}

简单的转换代码

// includes, system
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>

// includes, project
#include <cufft.h>
#include <helper_functions.h>
#include <helper_cuda.h>

// Complex data type
typedef float2 Complex;

// Filtering functions
void Convolve(const Complex *, int, const Complex *, int, Complex *);

// Padding functions
int PadData(const Complex *, Complex **, int,
const Complex *, Complex **, int);

////////////////////////////////////////////////////////////////////////////////
// declaration, forward
void runTest(int argc, char **argv);
// The filter size is assumed to be a number smaller than the signal size
#define SIGNAL_SIZE        32
#define FILTER_KERNEL_SIZE 11
#define PONTOS             32
#define PI           acos(-1)
#define dx        2*PI/PONTOS
        ////////////////////////////////////////////////////////////////////////////////
// Program main
////////////////////////////////////////////////////////////////////////////////
int main(int argc, char **argv)
{
runTest(argc, argv);
system("Pause");
}

////////////////////////////////////////////////////////////////////////////////
//! Run a simple test for CUDA
////////////////////////////////////////////////////////////////////////////////
void runTest(int argc, char **argv)
{
    printf("[simpleCUFFT] is starting...\n");

findCudaDevice(argc, (const char **)argv);

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

cufftComplex *h_reversed_signal = (Complex *)malloc(sizeof(Complex) * SIGNAL_SIZE);

// Initalize the memory for the signal
for (unsigned int i = 0; i < SIGNAL_SIZE; ++i)
{
    h_signal[i].x = rand() / (float)RAND_MAX;
    h_signal[i].y = sin(i*dx);;
}
cufftComplex *d_signal;
checkCudaErrors(cudaMalloc((void **)&d_signal, SIGNAL_SIZE*sizeof(cufftComplex)));
// Copy host memory to device
checkCudaErrors(cudaMemcpy(d_signal, h_signal, SIGNAL_SIZE*sizeof(cufftComplex),
    cudaMemcpyHostToDevice));
cufftHandle plan;
checkCudaErrors(cufftPlan1d(&plan, SIGNAL_SIZE, CUFFT_C2C, 1));

// Transform signal and kernel
printf("Transforming signal cufftExecC2C\n");
checkCudaErrors(cufftExecC2C(plan, (cufftComplex *)d_signal, (cufftComplex *)d_signal, CUFFT_FORWARD));
getLastCudaError("Kernel execution failed [ ComplexPointwiseMulAndScale ]");
// Transform signal back
printf("Transforming signal back cufftExecC2C\n");
checkCudaErrors(cufftExecC2C(plan, (cufftComplex *)d_signal, (cufftComplex *)d_signal, CUFFT_INVERSE));

// Copy device memory to host
checkCudaErrors(cudaMemcpy(h_reversed_signal, d_signal, SIGNAL_SIZE*sizeof(cufftComplex),
    cudaMemcpyDeviceToHost));
// check result
for (unsigned int i = 0; i < SIGNAL_SIZE; ++i)
{
    h_reversed_signal[i].x = h_reversed_signal[i].x / (float)SIGNAL_SIZE;
    h_reversed_signal[i].y /= (float)SIGNAL_SIZE;
    printf("first : %f %f  after %f %f \n", h_signal[i].x, h_signal[i].y, h_reversed_signal[i].x, h_reversed_signal[i].y);
    printf("1 Error %g %g \n", fabs(h_signal[i].x - h_reversed_signal[i].x), fabs(h_signal[i].y - h_reversed_signal[i].y));
}
bool bTestResult = sdkCompareL2fe((float *)h_reversed_signal, (float *)h_signal, 2 * SIGNAL_SIZE, 1e-5f);
//Destroy CUFFT context
checkCudaErrors(cufftDestroy(plan));
// cleanup memory
free(h_signal);
free(h_reversed_signal);
checkCudaErrors(cudaFree(d_signal));
cudaDeviceReset();
}

    // Pad data
    int PadData(const Complex *signal, Complex **padded_signal, int signal_size,
const Complex *filter_kernel, Complex **padded_filter_kernel, int filter_kernel_size)
{
int minRadius = filter_kernel_size / 2;
int maxRadius = filter_kernel_size - minRadius;
int new_size = signal_size + maxRadius;

// Pad signal
Complex *new_data = (Complex *)malloc(sizeof(Complex) * new_size);
memcpy(new_data + 0, signal, signal_size * sizeof(Complex));
memset(new_data + signal_size, 0, (new_size - signal_size) * sizeof(Complex));
*padded_signal = new_data;

// Pad filter
new_data = (Complex *)malloc(sizeof(Complex) * new_size);
memcpy(new_data + 0, filter_kernel + minRadius, maxRadius * sizeof(Complex));
memset(new_data + maxRadius, 0, (new_size - filter_kernel_size) * sizeof(Complex));
memcpy(new_data + new_size - minRadius, filter_kernel, minRadius * sizeof(Complex));
*padded_filter_kernel = new_data;

return new_size;
}

结果

[simpleCUFFT] is starting...
GPU Device 0: "GeForce GTX 570" with compute capability 2.0

Transforming signal cufftExecC2C
Transforming signal back cufftExecC2C
first : 0.001251 0.000000  after 0.001251 0.000000
first : 0.563585 0.195090  after 0.563585 0.195090
first : 0.193304 0.382683  after 0.193304 0.382683
first : 0.808740 0.555570  after 0.808740 0.555570
first : 0.585009 0.707107  after 0.585009 0.707107
first : 0.479873 0.831470  after 0.479873 0.831470
first : 0.350291 0.923880  after 0.350291 0.923879
first : 0.895962 0.980785  after 0.895962 0.980785
first : 0.822840 1.000000  after 0.822840 1.000000
first : 0.746605 0.980785  after 0.746605 0.980785
first : 0.174108 0.923880  after 0.174108 0.923879
first : 0.858943 0.831470  after 0.858943 0.831470
first : 0.710501 0.707107  after 0.710501 0.707107
first : 0.513535 0.555570  after 0.513535 0.555570
first : 0.303995 0.382683  after 0.303995 0.382683
first : 0.014985 0.195090  after 0.014985 0.195090
first : 0.091403 0.000000  after 0.091403 0.000000
first : 0.364452 -0.195090  after 0.364452 -0.195090
first : 0.147313 -0.382683  after 0.147313 -0.382683
first : 0.165899 -0.555570  after 0.165899 -0.555570
first : 0.988525 -0.707107  after 0.988525 -0.707107
first : 0.445692 -0.831470  after 0.445692 -0.831470
first : 0.119083 -0.923880  after 0.119083 -0.923879
first : 0.004669 -0.980785  after 0.004669 -0.980785
first : 0.008911 -1.000000  after 0.008911 -1.000000
first : 0.377880 -0.980785  after 0.377880 -0.980785
first : 0.531663 -0.923880  after 0.531663 -0.923879
first : 0.571184 -0.831470  after 0.571184 -0.831470
first : 0.601764 -0.707107  after 0.601764 -0.707107
first : 0.607166 -0.555570  after 0.607166 -0.555570
first : 0.166234 -0.382683  after 0.166234 -0.382683
first : 0.663045 -0.195090  after 0.663045 -0.195090
L2 =  3.00296e-007
4

1 回答 1

1

此代码序列是非法的：

for (unsigned int i = 0; i < SIGNAL_SIZE; ++i)
{
    d_signal[i].x = 2*d_signal[i].x;
    d_signal[i].y = 2*d_signal[i].y;
}

d_signal之前已创建为设备指针（通过cudaMalloc）。在主机代码中直接使用即取消引用设备指针是非法的。

解决此问题的一种可能方法是通过调用执行类似功能的 CUDA 内核来替换将复杂中间结果乘以 2 的这部分代码。

尽管它没有将中间结果精确地缩放 2，但有一个cuda 示例代码演示了 FFT->kernel->IFFT 的序列，可能有助于理解。

这是一个完整的示例，围绕您显示的代码构建：

$ cat t612.cu
# define SIGNAL_SIZE 64
# define PI acos(-1.0)
# define MY_X 2*PI/SIGNAL_SIZE

// includes, system
#include <stdlib.h>
#include <stdio.h>
#include <math.h>

// includes, project
#include <cufft.h>
#include <helper_cuda.h>


#define MPY_FACTOR 2.0
#define nTPB 256

__global__ void multiply_kernel(cufftComplex *data, int len){

  int idx = threadIdx.x+blockDim.x*blockIdx.x;
  if (idx < len){
    data[idx].x *= MPY_FACTOR;
    data[idx].y *= MPY_FACTOR;
  }
}

void runTest()
{
    printf("FFT multiply is starting...\n");

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

    cufftComplex *h_reversed_signal = (cufftComplex *)malloc(sizeof(cufftComplex) * SIGNAL_SIZE);

// Initalize the memory for the signal
  for (unsigned int i = 0; i < SIGNAL_SIZE; ++i)
  {
    h_signal[i].x = sin(i*MY_X);
    h_signal[i].y = 0;
  }
  cufftComplex *d_signal;
  checkCudaErrors(cudaMalloc((void **)&d_signal, SIGNAL_SIZE*sizeof(cufftComplex)));
// Copy host memory to device
  checkCudaErrors(cudaMemcpy(d_signal, h_signal, SIGNAL_SIZE*sizeof(cufftComplex), cudaMemcpyHostToDevice));
  cufftHandle plan;
  checkCudaErrors(cufftPlan1d(&plan, SIGNAL_SIZE, CUFFT_C2C, 1));

// Transform signal and kernel
  printf("Transforming signal cufftExecC2C\n");
  checkCudaErrors(cufftExecC2C(plan, (cufftComplex *)d_signal, (cufftComplex *)d_signal, CUFFT_FORWARD));
  getLastCudaError("Kernel execution failed [ ComplexPointwiseMulAndScale ]");
/*
  for (unsigned int i = 0; i < SIGNAL_SIZE; ++i)
  {
    d_signal[i].x = 2*d_signal[i].x;
    d_signal[i].y = 2*d_signal[i].y;
  }
*/
  multiply_kernel<<<(SIGNAL_SIZE+nTPB-1)/nTPB, nTPB>>>(d_signal, SIGNAL_SIZE);
  checkCudaErrors(cudaDeviceSynchronize());
  getLastCudaError("Kernel execution failed [ scaling kernel]");
// Transform signal back
  printf("Transforming signal back cufftExecC2C\n");
  checkCudaErrors(cufftExecC2C(plan, (cufftComplex *)d_signal, (cufftComplex *)d_signal, CUFFT_INVERSE));

   // Copy device memory to host
  checkCudaErrors(cudaMemcpy(h_reversed_signal, d_signal, SIGNAL_SIZE*sizeof(cufftComplex), cudaMemcpyDeviceToHost));
// check result
  for (unsigned int i = 0; i < SIGNAL_SIZE; ++i)
    {
            h_reversed_signal[i].x = h_reversed_signal[i].x / (float)SIGNAL_SIZE;
            h_reversed_signal[i].y = h_reversed_signal[i].y/(float)SIGNAL_SIZE;
            printf("first : %f %f  after %f %f \n", h_signal[i].x, h_signal[i].y, h_reversed_signal[i].x, h_reversed_signal[i].y);
    }
//Destroy CUFFT context
  checkCudaErrors(cufftDestroy(plan));
// cleanup memory
  free(h_signal);
  free(h_reversed_signal);
  checkCudaErrors(cudaFree(d_signal));
  cudaDeviceReset();
}

int main(){

  runTest();

}
$ nvcc -arch=sm_20 -I/usr/local/cuda/samples/common/inc -o t612 t612.cu -lcufft
$ ./t612
[simpleCUFFT] is starting...
Transforming signal cufftExecC2C
Transforming signal back cufftExecC2C
first : 0.000000 0.000000  after 0.000000 0.000000
first : 0.098017 0.000000  after 0.196034 -0.000000
first : 0.195090 0.000000  after 0.390181 0.000000
first : 0.290285 0.000000  after 0.580569 0.000000
first : 0.382683 0.000000  after 0.765367 0.000000
first : 0.471397 0.000000  after 0.942793 0.000000
first : 0.555570 0.000000  after 1.111140 -0.000000
first : 0.634393 0.000000  after 1.268787 0.000000
first : 0.707107 0.000000  after 1.414214 0.000000
first : 0.773010 0.000000  after 1.546021 0.000000
first : 0.831470 0.000000  after 1.662939 -0.000000
first : 0.881921 0.000000  after 1.763842 0.000000
first : 0.923880 0.000000  after 1.847759 0.000000
first : 0.956940 0.000000  after 1.913881 -0.000000
first : 0.980785 0.000000  after 1.961570 0.000000
first : 0.995185 0.000000  after 1.990369 -0.000000
first : 1.000000 0.000000  after 2.000000 0.000000
first : 0.995185 0.000000  after 1.990370 0.000000
first : 0.980785 0.000000  after 1.961570 -0.000000
first : 0.956940 0.000000  after 1.913881 0.000000
first : 0.923880 0.000000  after 1.847759 0.000000
first : 0.881921 0.000000  after 1.763842 0.000000
first : 0.831470 0.000000  after 1.662939 0.000000
first : 0.773010 0.000000  after 1.546021 0.000000
first : 0.707107 0.000000  after 1.414214 0.000000
first : 0.634393 0.000000  after 1.268786 -0.000000
first : 0.555570 0.000000  after 1.111140 0.000000
first : 0.471397 0.000000  after 0.942793 0.000000
first : 0.382683 0.000000  after 0.765367 0.000000
first : 0.290285 0.000000  after 0.580569 0.000000
first : 0.195090 0.000000  after 0.390181 0.000000
first : 0.098017 0.000000  after 0.196034 0.000000
first : 0.000000 0.000000  after -0.000000 0.000000
first : -0.098017 0.000000  after -0.196034 0.000000
first : -0.195090 0.000000  after -0.390181 0.000000
first : -0.290285 0.000000  after -0.580569 0.000000
first : -0.382683 0.000000  after -0.765367 0.000000
first : -0.471397 0.000000  after -0.942793 0.000000
first : -0.555570 0.000000  after -1.111140 0.000000
first : -0.634393 0.000000  after -1.268787 0.000000
first : -0.707107 0.000000  after -1.414214 0.000000
first : -0.773010 0.000000  after -1.546021 -0.000000
first : -0.831470 0.000000  after -1.662939 0.000000
first : -0.881921 0.000000  after -1.763842 0.000000
first : -0.923880 0.000000  after -1.847759 0.000000
first : -0.956940 0.000000  after -1.913881 0.000000
first : -0.980785 0.000000  after -1.961570 -0.000000
first : -0.995185 0.000000  after -1.990369 0.000000
first : -1.000000 0.000000  after -2.000000 -0.000000
first : -0.995185 0.000000  after -1.990370 -0.000000
first : -0.980785 0.000000  after -1.961570 0.000000
first : -0.956940 0.000000  after -1.913881 -0.000000
first : -0.923880 0.000000  after -1.847759 0.000000
first : -0.881921 0.000000  after -1.763842 -0.000000
first : -0.831470 0.000000  after -1.662939 -0.000000
first : -0.773010 0.000000  after -1.546021 0.000000
first : -0.707107 0.000000  after -1.414214 0.000000
first : -0.634393 0.000000  after -1.268786 0.000000
first : -0.555570 0.000000  after -1.111140 -0.000000
first : -0.471397 0.000000  after -0.942793 0.000000
first : -0.382683 0.000000  after -0.765367 0.000000
first : -0.290285 0.000000  after -0.580569 0.000000
first : -0.195090 0.000000  after -0.390181 -0.000000
first : -0.098017 0.000000  after -0.196034 0.000000
$
于 2015-01-30T14:59:33.693 回答