用户
 找回密码
 立即注册
quanzhang100 该用户已被删除
发表于 2013-11-19 12:26:21
2689612
#define ROW  256
#define COL  256
int main(int argc, char* argv[])
{

/***********************************************************************
*****************  256*256 2D FFT             **************************
/***********************************************************************/
        double *idata_h = (double *)malloc(ROW*COL*sizeof(double));
        memset(idata_h,0,ROW*COL*sizeof(double));

        double2 *idataC_h = (double2 *)malloc(ROW*COL*sizeof(double2));
        memset(idataC_h,0,ROW*COL*sizeof(double2));

        double2 *odata_h = (double2 *)malloc(ROW*COL*sizeof(double2));
        memset(odata_h,0,ROW*COL*sizeof(double2));

        double2 *odataC_h = (double2 *)malloc(ROW*COL*sizeof(double2));
        memset(odataC_h,0,ROW*COL*sizeof(double2));
       
        for (int i = 0; i < ROW*COL; i++)                               
        {
                idataC_h[i].x=idata_h[i] = rand();
                idataC_h[i].y=0;
        }
       
        double *idata_d;
        CudaSafeCall(cudaMalloc((void **)&idata_d, ROW*COL*sizeof(double)));
        CudaSafeCall(cudaMemset((void *)idata_d, 0 ,ROW*COL * sizeof(double)));

        double2 *odata_d;
        CudaSafeCall(cudaMalloc((void **)&odata_d, ROW*COL*sizeof(double2)));
        CudaSafeCall(cudaMemset((void *)odata_d, 0 ,ROW*COL * sizeof(double2)));

        double2 *idataC_d;
        CudaSafeCall(cudaMalloc((void **)&idataC_d, ROW*COL*sizeof(double2)));
        CudaSafeCall(cudaMemset((void *)idataC_d, 0 ,ROW*COL * sizeof(double2)));

        double2 *odataC_d;
        CudaSafeCall(cudaMalloc((void **)&odataC_d, ROW*COL*sizeof(double2)));
        CudaSafeCall(cudaMemset((void *)odataC_d, 0 ,ROW*COL * sizeof(double2)));


        //拷贝内存信号到显存
        CudaSafeCall(cudaMemcpy(idata_d, idata_h, ROW*COL * sizeof(double), cudaMemcpyHostToDevice));
        CudaSafeCall(cudaMemcpy(idataC_d, idataC_h, ROW*COL * sizeof(double2), cudaMemcpyHostToDevice));

        printf("2D ROW = %d COL = %d points fft start\n",ROW,COL);
        cufftHandle planD2Z;                                //创建CUFFT D2Z句柄
        checkCudaErrors(cufftPlan2d(&planD2Z, ROW,COL, CUFFT_D2Z));

        cufftHandle planZ2Z;                                //创建CUFFT Z2Z句柄
        checkCudaErrors(cufftPlan2d(&planZ2Z, ROW,COL, CUFFT_Z2Z));

        checkCudaErrors(cufftExecD2Z(planD2Z, idata_d, odata_d));       
        checkCudaErrors(cufftExecZ2Z(planZ2Z, idataC_d, odataC_d,CUFFT_FORWARD));

        //拷贝显存到内存信号
        CudaSafeCall(cudaMemcpy(odata_h,odata_d,ROW*COL * sizeof(double2),cudaMemcpyDeviceToHost));
        CudaSafeCall(cudaMemcpy(odataC_h,odataC_d,ROW*COL * sizeof(double2),cudaMemcpyDeviceToHost));
        //验证结果的正确性
        for(int i=0;i<ROW*COL;i++)
        {
                if ((fabs(odata_h[i].x - odataC_h[i].x)>0.1)||(fabs(odata_h[i].y - odataC_h[i].y)>0.1))
                {
                        printf("fft restult is wrong\n");       
                }
        }

        printf("2D ROW = %d COL = %d points fft finishes\n",ROW,COL);

        cudaFree(idata_d);
        cudaFree(odata_d);
        cudaFree(idataC_d);
        cudaFree(odataC_d);
        free(idata_h);
        free(odata_h);
        free(idataC_h);
        free(odataC_h);

        cudaDeviceReset();
        return 0;
}

以上代码是想完成idata_h是double型的矩阵,idataC_h是复数数据,其实部与idata_h相等,虚部为0,对idata_h和idataC_h做FFT,计算的结果最终存入odata_h,odataC_h中进行比较发现前128项相等,从129项就不等了?帮我看那看是不是程序哪里出了问题?

使用道具 举报 回复
发表于 2013-11-19 12:49:33
LZ您好:

FFT中,实值序列的FFT结果是共轭对称的,所以R2C或者D2Z的FFT函数利用了这个对称性,仅计算频域大概前一半点的数据(对应正频率的部分),所以您观察到的现象是正常的。

详情请参阅cuFFT 手册中的说明。

祝您编码顺利~
使用道具 举报 回复 支持 反对
发表于 2013-11-19 13:51:39
ice 发表于 2013-11-19 12:49
LZ您好:

FFT中,实值序列的FFT结果是共轭对称的,所以R2C或者D2Z的FFT函数利用了这个对称性,仅计算频域 ...

我的应用是对图像中的psf做FFT得到otf,本来是用cufftExecZ2Z,考虑到cufftExecD2Z耗时比cufftExecZ2Z少,想替换它,然而结果错了。那么请问一下如果是多实数做fft 用cufftExecD2Z怎么才能得到与cufftExecZ2Z一样的结果呢?
使用道具 举报 回复 支持 反对
发表于 2013-11-19 13:54:16
quanzhang100 发表于 2013-11-19 13:51
我的应用是对图像中的psf做FFT得到otf,本来是用cufftExecZ2Z,考虑到cufftExecD2Z耗时比cufftExecZ2Z少 ...

LZ您好:

请您复习一下FFT理论。

祝您好运~
使用道具 举报 回复 支持 反对
发表于 2013-11-19 16:46:32
ice 发表于 2013-11-19 13:54
LZ您好:

请您复习一下FFT理论。

cufftExecZ2Z变换之后的共轭对称的规律是什么呢?
如果要用cufftExecD2Z替换cufftExecZ2Z是不是还要把cufftExecD2Z计算的结果重新补齐(通过共轭对称)成和cufftExecZ2Z计算的结果呢?这样会不会得不偿失?
使用道具 举报 回复 支持 反对
发表于 2013-11-19 16:48:12
补充一下:上面说的是共轭对称的坐标的规律,我发现一维的挺好找规律的,对于二维好像不太好找。
使用道具 举报 回复 支持 反对
发表于 2013-11-19 16:58:33
quanzhang100 发表于 2013-11-19 16:46
cufftExecZ2Z变换之后的共轭对称的规律是什么呢?
如果要用cufftExecD2Z替换cufftExecZ2Z是不是还要把cuf ...

LZ您好:

1:共轭对称的具体形式请您查阅相关傅里叶变换理论的教本,解释这一点已经超出本版讨论范围。

2:R2C/D2Z的变换本来就是和C2R/Z2D的变换配合使用的,后者在计算中会自动从已有的复值数据中生成所需的另外一半共轭对称的数据,并予以计算。

3:如果您后续的自己的处理代码中需要Z2Z的格式的数据,那么您可以直接使用Z2Z变换,或用D2Z的数据手工补齐另外一半,或修改您后续的算法,从已有数据中就地生成所需的数据,并继续计算。您可以自行选择一种。

大致如此,请您定夺。

祝您编码顺利~
使用道具 举报 回复 支持 反对
发表于 2013-11-19 17:02:50
quanzhang100 发表于 2013-11-19 16:48
补充一下:上面说的是共轭对称的坐标的规律,我发现一维的挺好找规律的,对于二维好像不太好找。 ...

LZ您好:

1:根据FFTW的手册,在高维中,仅有一个维度的点数减半,cufft可能与之类似。具体规律请您查阅相关资料和手册。

2:您可以查阅手册看看有没有现成可用的函数,而无需自行推导。

祝您好运~
使用道具 举报 回复 支持 反对
发表于 2013-11-19 17:03:28
ice 发表于 2013-11-19 16:58
LZ您好:

1:共轭对称的具体形式请您查阅相关傅里叶变换理论的教本,解释这一点已经超出本版讨论范围。

好的 谢谢ice版主
使用道具 举报 回复 支持 反对
发表于 2013-11-19 17:05:02
quanzhang100 发表于 2013-11-19 17:03
好的 谢谢ice版主

不客气的,欢迎您常来
使用道具 举报 回复 支持 反对
12下一页
发新帖
您需要登录后才可以回帖 登录 | 立即注册