如何：CUDA IFFT
How to: CUDA IFFT
在 Matlab 中，当我输入一个一维复数数组时，我得到一个具有相同大小和相同维度的实数数组的输出。
尝试在 CUDA C 中重复此操作，但输出不同。
你能帮忙吗？在 Matlab 中，当我输入 ifft(array)
我的 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;
}
使用 MATLAB 计算 ifft 时，默认行为如下：
- 输入信号无零填充
- 输出信号无缩放
您的 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
在 Matlab 中，当我输入一个一维复数数组时，我得到一个具有相同大小和相同维度的实数数组的输出。 尝试在 CUDA C 中重复此操作，但输出不同。 你能帮忙吗？在 Matlab 中，当我输入 ifft(array)
我的 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;
}
使用 MATLAB 计算 ifft 时，默认行为如下：
- 输入信号无零填充
- 输出信号无缩放
您的 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