基于CUDA的离散傅里叶变换（Discrete Fourier Transform，DFT）

最近在做地震勘探的全波形反演，用分频反演的方法，需要对地震波场按照特定的频段进行傅里叶变换，这要用到DFT。在CPU上，DFT的计算非常耗时，当处理三维数据时耗时更加严重，所以，本人用CUDA+SU（seismic Unix），在GPU上来做DFT。话不多说，直接上代码：

程序代码：

#include <stdio.h>
#include "time.h"
#include "par.h"
#include "su.h"
#include "segy.h"
#include <cuda.h>
#include <cuda_runtime.h>
#include <helper_cuda.h>
#include "common.h"
#include "cufft.h"
#include <mpi.h>/****************************** self documentation ****************************/
const char *sdoc[] = {
"                                                                             ",
"                                                                             ",
NULL};
/** Authors:CUP** References:**/
/****************************** end self doc **********************************/__global__ void P_fft(int nx,int ny,int nf,cufftComplex *d_cf,float *d_p,cufftComplex *d_caf)
{unsigned int ix = threadIdx.x + blockIdx.x * blockDim.x;unsigned int iy = threadIdx.y + blockIdx.y * blockDim.y;unsigned int idy = ix * ny + iy;if (ix<nx && iy<ny ){d_p[idy]=idy;for (int i = 0; i < nf; ++i){d_caf[idy].x += d_cf[i].x*d_p[idy];d_caf[idy].y += d_cf[i].y*d_p[idy];}}
}int main(int argc, char **argv)
{time_t start,finish;double duration;int it;int nx,ny,nf,nt;int dimx=8,dimy=8;float damp=1;float dt=0.002;float t=0;float *invf;float *d_p;/* Transform real field to complex field */cufftComplex *h_cf;cufftComplex *d_cf;cufftComplex *d_caf;cufftComplex *h_caf;/* hook up getpar to handle the parameters */initargs(argc,argv);requestdoc(0);/* get required parameters */if (!getparint("nx",&nx))                   err("must specify nx!");if (!getparint("ny",&ny))                   err("must specify ny!");if (!getparint("nt",&nt))                   err("must specify nt!");;/* get optional parameters */if (!getparint("nt",&nt))                   nt = 0;if (!getparfloat("dt",&dt))                 dt = 1.0;nf = countparval("invf");warn(" nx=%d,ny=%d,nt=%d",nx,ny,nt);warn(" nt=%d;dt=%f",nt,dt);warn(" nf=%d",nf);invf = alloc1float(nf);getparfloat("invf", invf);for (int iff=0;iff<nf;iff++) {        warn("invf[%d]=%f",iff,invf[iff]);}h_cf = (cufftComplex *)malloc(sizeof(cufftComplex) * nf);    cudaMalloc((void**)&d_cf, sizeof(cufftComplex)*nf);  cudaMalloc((void**)&d_caf, sizeof(cufftComplex)*nx*ny);    cudaMalloc((void**)&d_p, sizeof(float)*nx*ny);  h_caf = (cufftComplex *)malloc(sizeof(cufftComplex)*nx*ny);CHECK(cudaMemset (d_p,   1, nx*ny*sizeof(float))); dim3 block(dimx,dimy);dim3 grid( (nx + block.x - 1) / block.x,(ny + block.y - 1) / block.y);for ( it = 0,t=0.0; it < nt; ++it,t+=dt){CHECK(cudaMemset (d_caf,  0, nx*ny*sizeof(cufftComplex))); for (int iff=0;iff<nf;iff++) {h_cf[iff].x=exp(-damp*t)*cos(2.0*PI*invf[iff]*t);h_cf[iff].y=exp(-damp*t)*sin(2.0*PI*invf[iff]*t);}CHECK(cudaMemcpy(d_cf,h_cf,sizeof(cufftComplex)*nf,cudaMemcpyHostToDevice));P_fft<<<grid,block>>>(nx,ny,nf,d_cf,d_p,d_caf);CHECK(cudaMemcpy(h_caf,d_caf,sizeof(cufftComplex)*nx*ny,cudaMemcpyDeviceToHost));}for (int i = 0; i < nx*ny; ++i){warn("h_caf[%d].x=%f",i,h_caf[i].x);warn("h_caf[%d].y=%f",i,h_caf[i].y);}free(h_cf);free(h_caf);free1float(invf);CHECK(cudaFree(d_cf));CHECK(cudaFree(d_caf));CHECK(cudaFree(d_p));return(0);
}

运行脚本：

#!/bin/sh
#
mpirun -np 1  SRC/complex_test nx=200 ny=100 nt=1000 dt=0.001 invf=0.5,1,1.5,2,2.5

此代码不限于地震反演，其他方向也可用。希望对大家有帮助，欢迎留言~~

基于CUDA的离散傅里叶变换（Discrete Fourier Transform，DFT）相关推荐

基于OpenCV完成离散傅里叶变换
基于OpenCV完成离散傅里叶变换目标学会使用函数: cv::copyMakeBorder() , cv::merge() , cv::dft() , cv::getOptimalDFTSize( ...
【OI向】快速傅里叶变换(Fast Fourier Transform)
[OI向]快速傅里叶变换(Fast Fourier Transform) 转存于本人博客园地址 FFT的作用在学习一项算法之前,我们总该关心这个算法究竟是为了干什么. (以下应用只针对OI ...
Discrete Fourier Transform离散傅里叶变换算法
DFT 公式: Ak=∑n=0N−1e−i2πNknan,k∈{0,1,...N−1}A_k = \sum_{n=0}^{N-1}e^{-i\frac{2\pi}{N}kn}a_n , k\in \{ ...
Discrete Fourier Transform（DCT）的理解
参考链接:简要理解DFT_我叫夏满满的博客-CSDN博客_dft原理 FFT(快速傅里叶变换)其本质就是DFT,只不过可以快速的计算出DFT结果.要理解FFT,必须先理解DFT(DiscreteFou ...
【matlab 图像处理】离散傅里叶变换离散余弦变换K-L变换小波变换
[matlab 图像处理]离散傅里叶变换&离散余弦变换&K-L变换&小波变换正交变换是信号处理的一种有效工具.图像信号不仅可以在空间域表示,也可以在频域表示,后者将有利于许多 ...
压缩感知的常见稀疏基名称及离散傅里叶变换基
题目:压缩感知的常见稀疏基名称及离散傅里叶变换基一.首先看九篇文献中有关稀疏基的描述: [1]喻玲娟,谢晓春.压缩感知介绍[J]. 电视技术,2008,32(12):16-18. 常用的稀疏基有:正 ...
【DSP】第4章离散傅里叶变换（NPU+BUAA）
离散傅立叶变换(DFT, Discrete Fourier Transform) 在学习之前先low一眼以前学过的傅里叶变换应用于连续周期信号--傅里叶技术展开 x a ( t ) = ∑ k = ...
dft计算傅里叶级数系数_DFT（离散傅里叶变换）与FFT（快速傅里叶变换）初识
一. 简介离散傅里叶变换(Discrete Fourier Transform, DFT)是数字信号处理最重要的基石之一,也是对信号进行分析和处理时最常用的工具之一.在200多年前法国数学家.物理学 ...
离散傅里叶变换（DFT）（一）
离散傅里叶变换(DFT) 傅里叶变换的优点就是能够将信号从时序空间转换到频域,从频率的角度去分析信号,能够容易发现一些时域内不容易发现的频率.这句话在N多博客里都有类似的描述.那么为什么呢?怎么转换到 ...

基于CUDA的离散傅里叶变换（Discrete Fourier Transform，DFT）

基于CUDA的离散傅里叶变换（Discrete Fourier Transform，DFT）相关推荐

最新文章

热门文章