CUDA_note_cuFFT

CUDA_note_cuFFT

Charles Lv7

cuFFT简介

FFT介绍

傅里叶变换是数字信号处理领域一个很重要的数学变换,它用来实现将信号从时域到频域的变换,在物理学、数论、组合数学、信号处理、概率、统计、密码学、声学、光学等领域有广泛的应用。离散傅里叶变换(Discrete Fourier Transform,DFT)是连续傅里叶变换在离散系统中的表示形式,由于DFT的计算量很大,因此在很长一段时间内其应用受到了很大的限制。20世纪60年代(1965年)由Cooley和Tukey提出了快速傅里叶变换(Fast Fourier Transform,FFT)算法,它是DFT的快速算法,使得离散傅里叶变换和卷积这类难度很大的计算工作的复杂度从$N^2$量级降到了$Nlog_2N$量级,大大提高了DFT的运算速度,从而使DFT在实际应用中得到了广泛的应用。

传统上,GPU只负责图形渲染,而大部分的处理都交给了CPU。自二十世纪九十年代开始,GPU的发展迅速。由于GPU具有强大的并行计算能力,加之其可编程能力的不断提高,GPU也用于通用计算,为科学计算的应用提供了新的选择。

2007年6月,NVIDIA公司推出了CUDA (Compute Unified Device Architecture),CUDA 不需要借助图形学API,而是采用了类C语言进行开发。同时,CUDA采用了统一处理架构,降低了编程的难度,同时,NVIDIA GPU引入了片内共享存储器,提高了效率。这两项改进使CUDA架构更加适合进行GPU通用计算。

快速傅里叶变换(FFT)

快速傅里叶变换(英语:Fast Fourier Transform, FFT),是快速计算序列的离散傅里叶变换(DFT)或其逆变换的方法。傅里叶分析将信号从原始域(通常是时间或空间)转换到頻域的表示或者逆过来转换。

FFT会通过把DFT矩阵分解为稀疏(大多为零)因子之积来快速计算此类变换。因此,它能够将计算DFT的复杂度从只用DFT定义计算需要的 $O(n^2)$,降低到$O(nlog n)$,其中n为数据大小。

快速傅里叶变换广泛的应用于工程、科学和数学领域。这里的基本思想在1965年才得到普及,但早在1805年就已推导出来。1994年美国数学家吉尔伯特·斯特朗把FFT描述为“我们一生中最重要的数值算法”,它还被IEEE科学与工程计算期刊列入20世纪十大算法。

FFT的CPU实现

一维FFT基2算法的实现

我们使用按频率抽取的方法实现了一维FFT基2算法。算法的关键代码如下:

(1) 声明双精度复数的结构体:

1
2
3
4
struct Complex {
double re; //复数的实部
double im; //复数的虚部
};

(2) 通过幂数获得快速傅里叶变换的长度,并初始化:

1
2
3
4
count = 1 << power; //power为幂数,count为快速傅里叶变换的长度
a = new Complex[count]; //a为初始化的数组,用来存放时域复数的值
b = new Complex[count]; //b为变换后存放结果的数组
memcpy(a, t, sizeof(Complex)*count); //初始化,将时域数据存放在a中,t为时域数据

(3) 计算旋转因子:

1
2
3
4
5
6
w = new Complex[count / 2];
for (i = 0; i<count / 2; i++) {
angle = -i*pi * 2 / count;
w[i].re = cos(angle);
w[i].im = sin(angle);
}

(4) 采用频率分解法进行蝶形运算:

1
2
3
4
5
6
7
8
9
10
11
12
13
for (k = 0; k<power; k++) {
for (j = 0; j<1 << k; j++) {
bfsize = 1 << (power - k);
for (i = 0; i<bfsize / 2; i++) {
p = j*bfsize;
b[i + p] = Add(a[i + p], a[i + p + bfsize / 2]); //Add()函数实现复数的加法
b[i + p + bfsize / 2] = Mul(Sub(a[i + p], a[i + p + bfsize / 2]), w[i*(1 << k)]);//Mul()函数实现复数的乘法,Sub()函数实现复数的减法
}
}
x = a;
a = b;
b = x;
}

(5) 蝶形运算全部完成后,对结果进行整序,从而得到正确的输出顺序:

1
2
3
4
5
6
7
8
for (j = 0; j<count; j++) {
p = 0;
for (i = 0; i<power; i++) {
if (j&(1 << i))
p += 1 << (power - i - 1);
}
f[j] = a[p];
}

二维FFT基2算法的实现

通过两次调用一维快速傅里叶变换即可实现二维快速傅里叶变换。关键代码如下:

(1) 计算进行傅里叶变换的宽度、高度,以及垂直方向上、水平方向上的迭代次数:

1
2
3
4
5
6
7
8
9
10
while (w * 2 <= width)  //计算进行傅立叶变换的宽度(2的整数次方)
{
w *= 2; //w为变换的宽度
wp++; //wp为垂直方向上的迭代次数
}
while (h * 2 <= height)//计算进行傅立叶变换的高度(2的整数次方)
{
h *= 2; //h为变换的高度
hp++; //hp为水平方向上的迭代次数
}

(2) 两次调用一维快速傅里叶变换:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
for (j = 0; j<h; j++)    //在垂直方向上进行快速傅立叶变换
{
FFT1D(&t[w*j], &f[w*j], wp);
}
for (j = 0; j<h; j++) //转换变换结果
{
for (i = 0; i<w; i++)
{
t[j + h*i] = f[i + w*j];
}
}
for (j = 0; j<w; j++) //水平方向进行快速傅立叶变换
{
FFT1D(&t[j*h], &f[j*h], hp);
}

FFT的GPU实现

对一维或多维信号进行离散傅里叶变换的FFT变换自身具有可“分治”实现的特点,因此能高效地在GPU平台上实现。CUDA提供了封装好的CUFFT库,它提供了与CPU上的FFTW库相似的接口,能够让使用者轻易地挖掘GPU的强大浮点处理能力,又不用自己去实现专门的FFT内核函数。使用者通过调用CUFFT库的API函数,即可完成FFT变换。

常见的FFT库在功能上有很多不同。有些库采用了基2变换,只能处理长度为2的指数的FFT,而其他的一些可以处理任意长度的FFT。CUFFT4.0提供了以下功能:

(1) 对实数或复数进行一维、二维或三维离散傅里叶变换;

(2) 可以同时并行处理一批一维离散傅里叶变换;

(3) 对任意维的离散傅里叶变换,单精度最大长度可达到6400万,双精度最大长度可达到12800万(实际长度受限于GPU存储器的可用大小);

(4) 对实数或复数进行的FFT,结果输出位置可以和输入位置相同(原地变换),也可以不同;

(5) 可在兼容硬件(GT200以及之后的GPU)上运行双精度的变换;

(6)支持流执行:数据传输过程中可以同时执行计算过程。

一维FFT算法的CUDA实现

关键代码如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
#include <cufft.h>          //CUFFT文件头
#define NX 1024
#define BATCH 1

cufftDoubleComplex *data; //显存数据指针

//在显存中分配空间
cudaMalloc((void**)&data, sizeof(cufftDoubleComplex)*NX*BATCH);

//创建CUFFT句柄
cufftHandle plan;
cufftPlan1d(&plan, NX, CUFFT_Z2Z, BATCH);

//执行CUFFT
cufftExecZ2Z(plan, data, data, CUFFT_FORWARD); //快速傅里叶正变换

//销毁句柄,并释放空间
cufftDestroy(plan);
cudaFree(data);

二维FFT算法的CUDA实现

二维FFT算法实现中,同一维FFT不同的是:

(1) 输入参数:没有BATCH,增加了NY。NX为行数,NY为列数;

(2) FFT的正变换的输入位置和输出位置不同。

关键代码如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
#include <cufft.h>          //CUFFT文件头

#define NX 1024
#define NY 1024

cufftDoubleComplex *idata, *odata; //显存数据指针

//在显存中分配空间
cudaMalloc((void**)&idata, sizeof(cufftDoubleComplex)*NX*NY);
cudaMalloc((void**)&odata, sizeof(cufftDoubleComplex)*NX*NY);

//创建CUFFT句柄
cufftHandle plan;
cufftPlan2d(&plan, NX, NY, CUFFT_Z2Z);

//执行CUFFT
cufftExecZ2Z(plan, idata, odata, CUFFT_FORWARD); //快速傅里叶正变换

//销毁句柄,并释放空间
cufftDestroy(plan);
cudaFree(idata);
cudaFree(odata);

三维FFT算法的CUDA实现

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
#define NX 64
#define NY 64
#define NZ 128

cufftHandle plan;
cufftComplex *data1, *data2;
cudaMalloc((void**)&data1, sizeof(cufftComplex)*NX*NY*NZ);
cudaMalloc((void**)&data2, sizeof(cufftComplex)*NX*NY*NZ);
/* Create a 3D FFT plan. */
cufftPlan3d(&plan, NX, NY, NZ, CUFFT_C2C);

/* Transform the first signal in place. */
cufftExecC2C(plan, data1, data1, CUFFT_FORWARD);

/* Transform the second signal using the same plan. */
cufftExecC2C(plan, data2, data2, CUFFT_FORWARD);

/* Destroy the cuFFT plan. */
cufftDestroy(plan);
cudaFree(data1); cudaFree(data2);

cuFFT介绍

cuFFT库提供了一个优化的且基于CUDA实现的快速傅里叶变换(FFT)。FFT在信号处理中可以将信号从时域转换到频域,逆FFT过程则相反。换句话说,一个FFT以规则的时间间隔接受信号中的序列样本并作为输入。然后使用这些样本生成一组叠加的分量频率,按频率抽取作为输入样本的信号。如下图所示,两个信号叠加生成信号cos(x)+cos(2x),并通过FFT将两个信号的分量转为频率1.0和2.0。

使用cuFFT API

cuFFT通常指两个独立的库:核心高性能的cuFFT库和可移植的cuFFTW库。cuFFT库是在CUDA中能提供自身API的FFT实现。另一方面,cuFFTW与标准的FFTW(快速傅里叶变换的标准C语言程序集)主机端FFT库有相同的API。和cuBLAS与传统的BLAS库共享大部分的API的情况类似,cuFFW则是用来最大限度地提高使用FFTW现有代码的可移植性。FFTW库的很多函数在cuFFTW中同样适用。此外,cuFFTW库假设所有要传输的输入数据都存储在主机内存中,并为用户处理所有的内存分布(cudaMalloc)和内存拷贝(cudaMemcpy)。虽然这可能对性能有所影响,但它大大加快了程序移植的过程。至于cuFFTW和cuFFT支持的操作,请参阅cuFFT用户指南。
  cuFFT库的配置是用FFT plan完成的。即cuFFT是用来指代它的操作术语。一个plan定义了一个要进行的单一变换操作。cuFFT使用plan来获取内存分配、内存转移以及内存启动来执行变换请求。不同的plan创建函数可以用来生成增大复杂性和维数的plan:

1
2
3
cufftResult cufftPlan1d(cufftHandle *plan,int nx,cufftType type,int batch);
cufftResult cufftPlan2d(cufftHandle *plan,int nx,int ny,cufftType type);
cufftResult cufftPlan3d(cufftHandle *plan,int nx,int ny,int nz,cufftType type);

cuFFT还支持多种输入和输出数据类型,包括以下几种:复数到复数、实数到复数、复数到实数。当然,对于许多实际的应用程序,最实用的是实数到复数这种类型,它允许我们从实际系统输入实际的测量结果到cuFFT中。
  一旦配置好一个cuFFT plan,使用cufftExec* 函数来对它进行调用执行(例如,cufftExecC2C)。一般来说,无论该变换是一种正向FFT(时域到频域)还是怒向FFT(频域到时域),函数调用都可以将plan、输入数据的存储位置、输出数据的存放位置作为输入。

cuFFT功能示范

一个cuFFT应用程序的工作流的不同取决于变换的复杂性。一个cuFFT应用程序的工作流一般应包括:
  1.创建并配置一个cuFFT plan;
  2.用cudamalloc函数分配设备内存来存储输入样本和输出频率。注意,所分配的内存必须支持对应执行的变换类型(例如,复数到复数、实数到复数、复数到实数)我们可以使用相同的设备内存对输入和输出直接进行变换;
  3.用cudaMemcpy用输入样本传送设备内存;
  使用cudaExec* 函数执行plan;
  用cudaMemcpy取回设备内存中的结果;
  用cudaFree和cufftDestroy释放CUDA和cuFFT资源;
  接下来的示例代码从函数cos(x)中生成一个输入样本序列,将其转换为复数后传给GPU,在将结果拷贝回主机端之前执行复数到复数的一维plan。值得注意的是,对于输入和输出参数,因为可以将相同的内存位置dComplexSamples传给cufftExecC2C,所以这是一个就地FFT运算。可以预先分配一个独立的输出缓存区,并用来存储输出结果。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
// Setup the cuFFT plan
cufftPlan1d(&plan, N, CUFFT_C2C, 1);

// Allocate device memory
cudaMalloc((void **)&dComplexSamples,
sizeof(cufftComplex) * N);

// Transfer inputs into device memory
cudaMemcpy(dComplexSamples, complexSamples,
sizeof(cufftComplex) * N, cudaMemcpyHostToDevice);

// Execute a complex-to-complex 1D FFT
cufftExecC2C(plan, dComplexSamples, dComplexSamples,
CUFFT_FORWARD);

// Retrieve the results into host memory
cudaMemcpy(complexFreq, dComplexSamples,
sizeof(cufftComplex) * N, cudaMemcpyDeviceToHost);

FFT详解&大数乘法

参考链接:FFT详解&大数乘法

首先是一道非常好的FFT解决大数乘法的练手题——P1919 【模板】A*B Problem 升级版(FFT 快速傅里叶变换)

提供一个讲解非常清楚的视频——超硬核FFT快速傅里叶变换讲解,高效进行高精度乘法运算!_哔哩哔哩_bilibili

底下提供一个我用C++写的题解,算是勉强对FFT使用入了门。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
#include<bits/stdc++.h>
using namespace std;
typedef complex<double> cp;
#define N 2097153
const double pie = acos(-1);
int n;
cp a[N], b[N];
int rev[N], ans[N];
char s1[N], s2[N];

//初始化每个位置最终到达的位置
void init(int k) {
int len = 1 << k;
for (int i = 0; i < len; i++)
rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (k - 1));
}

//a表示要操作的系数,n表示序列长度
//若flag为1,则表示FFT,为-1则为IFFT(需要求倒数)
void fft(cp *a, int n, int flag) {
for (int i = 0; i < n; i++) {
//i小于rev[i]时才交换,防止同一个元素交换两次,回到它原来的位置。
if (i < rev[i])swap(a[i], a[rev[i]]);
}
for (int h = 1; h < n; h *= 2)//h是准备合并序列的长度的二分之一
{
cp wn = exp(cp(0, flag * pie / h));//求单位根w_n^1
for (int j = 0; j < n; j += h * 2)//j表示合并到了哪一位
{
cp w(1, 0);
for (int k = j; k < j + h; k++)//只扫左半部分,得到右半部分的答案
{
cp x = a[k];
cp y = w * a[k + h];
a[k] = x + y; //这两步是蝴蝶变换
a[k + h] = x - y;
w *= wn; //求w_n^k
}
}
}
//判断是否是FFT还是IFFT
if (flag == -1)
for (int i = 0; i < n; i++)
a[i] /= n;
}

int main() {
scanf("%s%s", s1, s2);
int len1 = strlen(s1);
int len2 = strlen(s2);
n = len1 + len2 - 1;
//读入的数的每一位看成多项式的一项,保存在复数的实部
for (int i = 0; i < len1; i++)a[i] = (double) (s1[len1 - i - 1] - '0');
for (int i = 0; i < len2; i++)b[i] = (double) (s2[len2 - i - 1] - '0');
//k表示转化成二进制的位数
int k = 1, s = 2;
while ((1 << k) < n)k++, s <<= 1;
init(k);
//FFT 把a的系数表示转化为点值表示
fft(a, s, 1);
//FFT 把b的系数表示转化为点值表示
fft(b, s, 1);
//FFT 两个多项式的点值表示相乘
for (int i = 0; i < s; i++)
a[i] *= b[i];
//IFFT 把这个点值表示转化为系数表示
fft(a, s, -1);
//保存答案的每一位(注意进位)
for (int i = 0; i < s; i++) {
//取实数四舍五入,此时虚数部分应当为0或由于浮点误差接近0
ans[i] += (int) (a[i].real() + 0.5);
ans[i + 1] += ans[i] / 10;
ans[i] %= 10;
}
while (!ans[s] && s > -1)s--;
if (s == -1)printf("0");
else
for (int i = s; i >= 0; i--)
printf("%d", ans[i]);
return 0;
}

参考资料

[1]. CUDA学习笔记3:CUFFT

[2].cuFFT (nvidia.com)

[3].cuFFT库_伴君的博客-CSDN博客

[4].FFT详解&大数乘法

[5].P1919 【模板】A*B Problem 升级版(FFT 快速傅里叶变换)

  • Title: CUDA_note_cuFFT
  • Author: Charles
  • Created at : 2023-07-28 14:16:54
  • Updated at : 2024-05-03 09:19:01
  • Link: https://charles2530.github.io/2023/07/28/cuda-note-cufft/
  • License: This work is licensed under CC BY-NC-SA 4.0.
Comments